A Synchronous Bisimulation Based Approach for Information Flow Analysis

中科院计算机类SCI期刊及分区（201
期刊 IEEE TRANSACTIONS ON FUZZY SYSTEMS International Journal of Neural Systems IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION INTEGRATED COMPUTER-AIDED ENGINEERING IEEE Transactions on Cybernetics IEEE Transactions on Neural Networks and Learning Systems MEDICAL IMAGE ANALYSIS Information Fusion INTERNATIONAL JOURNAL OF COMPUTER VISION IEEE TRANSACTIONS ON IMAGE PROCESSING IEEE Computational Intelligence Magazine EVOLUTIONARY COMPUTATION IEEE INTELLIGENT SYSTEMS PATTERN RECOGNITION ARTIFICIAL INTELLIGENCE KNOWLEDGE-BASED SYSTEMS NEURAL NETWORKS EXPERT SYSTEMS WITH APPLICATIONS Swarm and Evolutionary Computation APPLIED SOFT COMPUTING DATA MINING AND KNOWLEDGE DISCOVERY INTERNATIONAL JOURNAL OF APPROXIMATE REASONING SIAM Journal on Imaging Sciences DECISION SUPPORT SYSTEMS Swarm Intelligence Fuzzy Optimization and Decision Making IEEE TRANSACTIONS ON KNOWLEDGE AND DATA ENGINEERING JOURNAL OF MACHINE LEARNING RESEARCH ACM Transactions on Intelligent Systems and Technology NEUROCOMPUTING ENGINEERING APPLICATIONS OF ARTIFICIAL INTELLIGENCE CHEMOMETRICS AND INTELLIGENT LABORATORY SYSTEMS ARTIFICIAL INTELLIGENCE IN MEDICINE COMPUTER VISION AND IMAGE UNDERSTANDING JOURNAL OF AUTOMATED REASONING INTERNATIONAL JOURNAL OF INTELLIGENT SYSTEMS COMPUTATIONAL LINGUISTICS ADVANCED ENGINEERING INFORMATICS JOURNAL OF INTELLIGENT MANUFACTURING Cognitive Computation IEEE Transactions on Affective Computing JOURNAL OF CHEMOMETRICS MECHATRONICS IEEE Transactions on Human-Machine Systems Semantic Web IMAGE AND VISION COMPUTING Wiley Interdisciplinary Reviews-Data Mining and Knowledge Discovery issn 1063-6706 0129-0657 0162-8828 1089-778X 1069-2509 2168-2267 2162-237X 1361-8415 1566-2535 0920-5691 1057-7149 1556-603X 1063-6560 1541-1672 0031-3203 0004-3702 0950-7051 0893-6080 0957-4174 2210-6502 1568-4946 1384-5810 0888-613X 1936-4954 0167-9236 1935-3812 1568-4539 1041-4347 1532-4435 2157-6904 0925-2312 0952-1976 0169-7439 0933-3657 1077-3142 0168-7433 0884-8173 0891-2017 1474-0346 0956-5515 1866-9956 1949-3045 0886-9383 0957-4158 2168-2291 1570-0844 0262-8856 1942-4787

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.。

中国科学:信息科学2014年第44卷第?期:1–? 中国高速铁路列控系统的形式化分析与验证郭丹青x,吕继东y,王淑灵x,唐涛y,詹乃军x∗,周达天y,邹亮xx中国科学院软件研究所计算机科学国家重点实验室,北京100190y北京交通大学,北京100044*詹乃军.E-mail:znj@收稿日期:2014–01–19;接受日期:2014–05–08国家重点基础研究发展计划(973计划)(批准号:2014CB340700),国家高技术研究发展计划(863计划)(批准号:2012AA112801)和国家自然科学基金(批准号:91118007,6110006,61304185)项目资助摘要高速铁路列控系统的安全与否直接涉及人民的生命财产安全,对高速铁路列控系统进行严格的形式化验证具有重要意义.但是随着高速铁路列控系统软件以及硬件规模的不断增大,系统的复杂性有了很大的提高,直接对高速铁路列控系统进行形式化验证已经变得越来越困难.另一方面,由于图形化建模和仿真表现方式直观,易于理解,在工程实践中已经得到了广泛的应用.因此,为了更好的保证铁路系统的安全,对系统进行仿真排除部分安全隐患显得尤为重要.本文通过使用Simulink/Stateflow建模工具对高速铁路列控系统的行车许可,等级升级及部分模式转换场景进行了建模.该模型具有普适性,通过修改参数信息,可以对不同的等级转换以及模式转换的组合情况进行仿真.通过使用该模型我们对10个组合情况进行了仿真,最后发现在某些情况下可能会出现不正常停车或者等级转换失败的现象.类似于测试,仿真仅仅能够发现错误,不能证明系统没有错误.因为仿真的这种不完备性,对仿真辅助形式验证在安全攸关系统设计中是非常必要的.为此,我们取其中一个不正常停车的场景进行了形式验证,验证结果证明了在任何情况下都会不正常停车.关键词中国高速铁路列控系统Simulink/Stateflow仿真模式转换等级转换形式化验证1引言高速铁路具有节能,环保,大运量,安全舒适等明显优势,是交通运输体系中最具可持续性和环境友好性的运输模式,客运高速化已成为世界潮流.我国高速铁路建设列为拉动内需的“火车头”,重点发展高速铁路,将建成以北京为中心,2小时大中城市圈为主节点,8小时快速铁路交通为主干的高速铁路客运网络.可以预见,以高速铁路为核心的快速铁路交通网的建成,将使我国经济社会发展步入高速铁路时代.高速铁路是庞大复杂的系统工程,集成了多学科,多领域的高新技术.其中中国高速铁路列控系统是保证列车高速,高效,安全运行的坚实基础,所以确保高速铁路列控系统的安全性有非常重要的意义.列车运行控制系统是高速铁路的核心技术之一,是一个深度融合了计算,通信和控制的系统,同时又是安全攸关的系统,即系统的任何错误都可能导致灾难性后果,是一种典型的信息物理融合系统(Cyber-Physical Systems,CPS).例如,2011年7月23日发生在温州事故导致40多位乘客失去生命.它郭丹青,吕继东,唐涛等:中国高速铁路列控系统的形式化分析与验证通过3C(Computation,Communication,Control)技术的有机融合与深度协作,实现列车运行过程与信息交互系统的实时感知,动态控制和信息服务.随着列车运行速度的不断提高(运行时速将超过500公里/小时),地面,轨道网络与列车之间,列车与列车之间的交互作用极其复杂,使得列车运行控制系统建模,性质刻画和性能抽取变得复杂,最终必然导致对系统安全性分析和验证的复杂程度急剧增加.列车运行安全控制标准是列车运行控制系统的技术标准,是实现列车安全运行的技术保障之一.因而必须确保制定的列车运行安全控制标准是安全可靠的,即首先保证内部不存在互相矛盾的地方;其次,要保证根据标准列车能够安全运行在给定路段上;再次,一旦发生事故,根据标准能够采取相应应对措施不至于发生灾难性后果.1.1本文主要贡献本文主要研究如何利用信息物理融合系统方面最新结果,特别是我们最近几年在混成系统建模,分析和验证方面的结果,提出一套中国高速铁路列控系统图形建模,仿真,形式建模和形式验证于一体的方法,从而提高中国高速铁路列控系统的可靠性.具体讲：Simulink/Stateflow[1,2]是Matlab中一个重要的商业组件,拥有Matlab的强大计算能力支持,在工业界使用广泛,拥有深厚的用户基础.它是一种图形化的建模工具,对模型的描述比较直观,并且支持仿真,可以检查模型是否与预计行为相符,以便排除部分早期的设计失误,已经成为一种事实上的工业标准.同时,它还可以将建立的模型转化为对应的C语言,有助降低开发成本.而过去,很少有使用Simulink/Stateflow对列控系统进行建模和仿真的工作.基于此,我们首先针对目前中国高速铁路列控系统的规范,提出一种一般性Simulink/Stateflow图形模型框架.对高速铁路列控系统的行车许可,等级升级及部分模式转换场景建立了Simlink/Stateflow模型,并且使用该模型对二级到三级与模式转换相结合的部分情况进行了仿真,发现在部分情况下会出现不正常停车以及等级转换失败的现象.另一方面,因为仿真仅仅能够根据环境的有穷输入在有限时间内观察系统行为是否满足要求,而环境的输入通常有无穷多种,并且这些系统的运行可能没有时间限制,所以类似于测试,仿真仅仅能够发现系统错误,不能够证明系统没有错误.特别对于安全攸关系统,仅仅仿真无法保证系统的正确性.因而对Simulink/Stateflow图形模型进行形式验证,作为仿真的一种补充,是非常必要的.形式验证Simulink/Stateflow图形模型的前提是需要提供一个形式语义以及针对该形式语义的验证技术.在[3,4]中,我们考虑如何将Simulink/Stateflow图形模型转换成HCSP[5,6]形式模型,并使用HHL[7]及其定理证明器[8]来形式验证转化后的形式模型.本文的另一贡献是：我们以其中某个场景为例,说明如何使用上述方法将其转换为形式模型并严格验证,验证结果表明在任何输入下均不能正常停车.1.2相关工作为了支持混成系统建模,人们提出很多建模语言,例如:Modelica[9],HybridUML[10],时间自动机[11],混成自动机[12],Esterel[13],等等.Modelica[9]与Simulink同是图形化的建模语言.它是一种开源语言,支持用户自定义修改.它同时支持图形建模和代码生成,表达能力很强.但是它没有类似于Stateflow这样的控制流图建模工具,所以对控制逻辑建模支持比较薄弱.虽然Modelica使用算法和函数对控制逻辑进行建模,但是由于这种建模方式是基于文本的,表现不够直观,不利于对高层模型进行建模.HybridUML[10]是一种对混成系统进行图形化建模的语言,它是基于传统UML的混成扩展. UML在工业界有广泛的使用,用户基础较好.但是HybridUML不能进行仿真,也不能进行验证,使用上有很大的限制.混成自动机[12]是自动机在混成系统方面的拓展,基于可达集计算的各种验证工作在学2中国科学:信息科学第44卷第?期术界研究深入,应用广泛.基于自动机的建模技术直观,易于理解,但是它没有系统结构方面的信息,组合性差,不适合大型系统的建模.特别是,没有图形化工具支持,不能进行仿真,形式验证技术也不适用于复杂大型系统.Esterel[13]是一种基于时间同步模型的建模语言.与常见的同步异步通讯方式不同,当Esterel中信号被发送之后,该时间点上其他并发进程都可以将该信号接收任意有限次.Esterel拥有成熟的商业工具Scade[14]对其进行支持,该工具生成的代码可适用于实际系统,减少了后期开发的工作量.总之,Simulink/Stateflow是更实用的混成系统建模,分析和仿真工具.国内外有许多关于高速铁路列控系统形式建模和验证的工作,例如:Platzer等在[15]利用一种基于微分不变量的微分动态逻辑等技术对欧洲高速铁路列控系统进行建模,分析和验证,从欧洲高速铁路列控系统的非形式描述中发现了许多不严格的地方,之后他们还进一步通过设计一个可以严格证明安全性质的形式模型,以改进欧洲高速铁路列控系统的设计.国内也有很多科研工作者尝试对中国高速铁路列控系统进行形式建模和验证,例如:在[16]中,作者结合UML在业界的广泛应用及SMV模型检验工具的优点,提出了一套列控系统规范的建模与验证方法.而在[17]中,作者使用Petri网对高速铁路列控系统信道模型和数据传输的时间特性进行建模和分析,为CTCS-3级列控系统规范中的相关参数设计提供了前提和基础.综上分析,目前还没有一套针对高速铁路列控系统的集图形建模,仿真,形式建模和验证于一体的方法.本文填补了这方面的空白.2背景知识在这一节中,我们将介绍高速铁路列控系统,Simulink/Stateflow和HCSP以及HHL的一些背景知识,因为我们关心等级转换和模式转换相结合的场景,所以我们使用第一小节和第二小节分别介绍高速铁路列控系统等级和列车在不同等级下的不同模式,第三小节介绍Simulink/Stateflow,第四小节介绍HCSP和HHL.2.1高速铁路列控系统等级概述中国列车控制系统CTCS(Chinese Train Control System),根据总体原则,从国情,路情实际出发,共划分为5级[18],其中我们关心如下两个等级:1.CTCS-2(C2)级基于轨道传输信息的列车运行控制系统,面向提速干线和高速新线,采用车,地一体化设计.适用于各种限速区段,地面可不设通过信号机,机车凭车载信号行车.地面子系统中增加列控中心,根据列车占用情况及进路状态,计算行车许可及静态列车速度曲线并传送给列车.点式信息设备用于向车载设备传输定位信息,进路参数,线路参数,限速和停车信息等.2.CTCS-3(C3)级基于无线传输信息,采用轨道电路等方式检查列车占用的列车运行控制系统.面向提速干线,高速新线或特殊线路,基于无线通信的固定闭塞或虚拟自动闭塞.CTCS-3在CTCS-2的基础上有了很多改进.主要体现在控制中心,和信息交互方式的改变上.3郭丹青,吕继东,唐涛等:中国高速铁路列控系统的形式化分析与验证2.2CTCS-2级和CTCS-3级下模式概述每个等级下列车都有不同的模式,其中我们只关心C2下的完全监控模式（Full Supervision mode, FS）,目视行车模式(OnSight mode,OS),部分监控模式(Partial Supervision mode,PS)模式,以及C3下的FS,OS,引导模式(Calling on mode,CO),冒进模式(Trip mode,TR).1.CTCS-2下模式介绍完全监控模式(FS):在完全监控模式下,列控车载设备应能判断列车位置和停车位置,在保证列车速度满足线路固定限速,车辆构造速度,停车位置,临时限速等条件的前提下,生成目标距离连续速度控制模式曲线,并连续监控列车速度,与模式速度比较,自动输出紧急制动或常用制动命令,同时,应能通过DMI显示列车实际速度,允许速度,目标速度和目标距离等信息.目视行车模式(OS):车载设备显示停车信号或位置不确定时,在停车状态下司机按压专用按钮可使车载设备转入目视行车模式.在该模式下,列控车载设备生成NBP(normal brake profile)∗为20km/h的模式曲线.部分监控模式(PS):由于应答器信息接收异常导致线路数据缺失,或者由于其它原因列控车载设备无线路数据,以及侧线接车和在车站办理引导接车时,列控车载设备的工作模式都定义为部分监控模式.2.CTCS-3下模式介绍完全监控模式(FS):当车载设备具备列车控制所需的全部基本数据(包括列车数据,行车许可和线路数据等)时,车载设备生成目标距离连续速度控制模式曲线,并通过DMI显示列车运行速度,允许速度,目标速度和目标距离等信息,监控列车安全运行.目视行车模式(OS):车载设备显示停车信号或位置不确定时,在停车状态下司机按压专用按钮可使车载设备转入目视行车模式.目视行车模式下,车载设备按固定限制速度40km/h监控列车运行,列车每运行一定距离司机需确认一次.引导模式(CO):当开放引导信号进行接发车时,车载设备生成目标距离连续速度控制模式曲线,并通过DMI显示列车运行速度,允许速度,目标速度和目标距离等,车载设备按固定限制速度40km/h监控列车运行,司机负责在列车运行时检查轨道占用情况.冒进模式(TR):列车执行冒进防护时,车载设备进入冒进模式.这个模式下列车会紧急制动,直至完全停止.2.3Simulink/Stateflow基本概念2.3.1Simulink介绍Simulink[1]是MATLAB最重要的组件之一,它提供一个图形化的动态系统建模,仿真和综合分析的集成环境.在该环境中,无需大量书写程序,而只需要通过简单直观的鼠标操作,就可构造出复杂的系统.Simulink是用于动态系统和嵌入式系统的多领域仿真和基于模型的设计工具,拥有丰富的可扩充预定义模块库.对各种实时系统,包括通讯,控制,信号处理,视频处理和图像处理系统,Simulink提供了交互式图形化环境和可定制模块库来对其进行设计,仿真,执行和测试.∗常用制动介入曲线4中国科学:信息科学第44卷第?期Simulink模型是由一系列的模块和连接这些模块的边组成的.每个模块可以是一个Simulink模块库中的基本模块,表达了输入输出之间的一个数学关系,也可以是由若干模块组成的子系统.连接边反应了两个被连接模块之间的关系.Simulink模型通常由三部分组成:输入信号源(Source),系统(System)以及接收模块(Sink).1.输入信号源模块库(Source)输入信号源模块库用来向模型提供输入信号,没有输入口,但是至少有一个输出口.主要有: Constant,Step,Ramp,Sine Wave,Signal Generator,From File,From Workspace,Clock,In.2.系统模块库(System)系统模块主要包括连续模块,离散模块,数学模块,逻辑模块,信号处理模块,系统模块.连续模块主要有:Integrator,Derivative,State-Space,Transfer Fcn,Zero-Pole,Transport Delay.离散模块主要有:Discrete-time Integrator,Discrete Filter,Discrete State-Space,Discrete Transfer-Fcn,Discrete,Zero-Pole,First-Order Hold,Zero-Order Hold,Unit Delay.数学模块主要有:Sum, Product,Dot Product,Gain,Math Function,MinMax,Abs,Sign.逻辑模块主要有:Compare To Zero,Compare to Constant,Inteval Test,Relational Operator,Logical Operator,Bit Set,Bitwise Operator,Bit Clear.信号处理模块主要有:IC,Bus Creator,Bus Selector,Mux,Demux,Switch, Multiport Switch,Manual Switch,Selector.系统模块主要有:SubSystem,Triggered Subsystem, Enabled Subsystem,Function-Call Subsystem,If Action Subsystem.3.接收模块库(Sink)接收模块是用来接收模块信号的,没有输出口,但是至少有一个输入口.主要有:Scope,Dis-play,XY Graph,To File,To Workspace,Stop Simulation,Out.Simulink模块从其行为模式上分可以分为连续模块和离散模块两种.Simulink模型中的每个模块都有一个样本时间,它的取值范围是-1或者非负浮点数.如果这个模块的样本时间是-1,表示其样本时间由其来源模块的样本时间计算而来;如果这个模块的样本时间是0,表示这个模块是连续模块,其计算行为具有连续性;如果是一个正值x,表示该模块为离散模块,其每隔x时间长度重新计算一次.大多数Simulink模块都包含配置参数,用户可以通过修改这些参数来获得预想的功能.比如说常数模块可以设置输出常数值,积分模块可以设置初始值,Switch模块可以设置操作符号类型和阈值.图1描述了一个Simulink模型,其中的Add模块从In1,In2接收了两个输入,通过加法运算产生结果,并通过Out1输出.该模型中In1,In2,Out1和Add分别对应于输入信号源,输出信号源与系统.2.3.2Stateflow介绍Stateflow[2]是Simulink里面的一个工具箱,它是一个基于状态机和流程图来构建组合和时序逻辑决策模型并进行仿真的环境.Stateflow可以将图形表示和表格表示(包括状态转换图,流程图,状态转换表和真值表)结合在一起,针对系统对事件,基于时间的条件以及外部输入信号的反应方式进行建模.用户可以在使用Simulink仿真时,使用这种图形化的工具实现各个状态之间的转换,对复杂的监控逻辑进行建模.Stateflow模型主要由事件变量集合,状态,结点和状态迁移组成:5郭丹青,吕继东,唐涛等:中国高速铁路列控系统的形式化分析与验证图1Simulink 图表Figure 1A Simulink diagram 图2Stateflow 图表Figure 2A Stateflow diagram1.事件变量集合事件变量集合是记录Stateflow 模型中使用的广播事件和变量的集合.广播事件和变量均分为输入输出和局部三种,分别表示从外部环境输入的广播事件或变量,向外部输出的广播事件或变量和内部使用的广播事件或变量.2.状态状态代表了系统现在处在的情况.在Stateflow 下,状态有两种行为:活动的(active)和非活动的(inactive).可以对一个状态进行标记,包括给这个状态规定状态名以及这个状态在进入,退出和处于激活状态下接收到一个事件所应执行的动作,一般表示如下:nameentry:entry actionsduring:during actionsexit:exit actionsbind:data and eventson event_name:on event_name actions入口动作(entry:entry actions)是表示发生状态迁移,激活了该状态时需要执行的动作.中间动作(during:during actions)表示原处于激活的状态受到一个事件的触发,不存在从这个状态发出的状态迁移时,此状态仍处于激活状态需要执行的动作.出口动作(exit:exit actions)表示存在由此状态发出的有效状态迁移时,该状态退出时执行的动作.数据事件绑定动作(bind:data and events)将数据和事件绑定在此状态上.绑定的数据只能在此状态或其子状态内被改写,其他状态只能读取此数据.绑定的事件由此状态或其子状态广播.特定事件发生动作(on event_name:on event_name actions).event_name 规定一个特定的事件;on event_name actions 表示当该状态是激活状态且event_name 规定的事件发生时需要执行的动作.3.结点结点用于描述状态迁移过程中的迁移信号的分离和汇合.通过使用结点,状态和状态之间的迁移不仅仅是边的简单连接,而是一个复杂的迁移网络.6中国科学:信息科学第44卷第?期4.状态迁移状态迁移†一般用于连接两个状态或者结点,它可以被标记,该标记的一般形式为如下四元组: Event Trigger[Condition]Condition Action/Transition Action,触发事件(Evnet trigger)是事件集合中的一个元素,作为迁移条件的一部分.条件(Condition)是一个布尔表达式,与触发事件共同组成状态迁移的迁移条件.触发事件(Evnet trigger)为空或者与当前广播事件一致,同时条件(Condition)为真,则该状态迁移被触发.条件动作(Condition Action)和迁移动作(Transition Action)是使用Matlab的action language书写的代码片段.条件动作在迁移被触发时被立即执行,而迁移动作仅当迁移终点为状态时被立即执行,否则被寄存在一个队列中,若该迁移最终达到一个状态,队列中的迁移动作将被顺序执行.Stateflow模型具有层次化的特性,任意Stateflow状态均可嵌入一个Stateflow子模型来丰富该状态的行为.同一层次中,所有的状态是互斥(OR)或者并行(AND)的,互斥状态之间可以使用迁移进行连接,而并行状态之间不能使用迁移连接.Stateflow通讯采用广播机制,当一个事件被广播时,并行状态按照其预先定义的顺序先后接收到该广播事件,互斥状态仅有激活状态收到该广播事件.当状态接收到一个广播事件之后,进行如下操作:按照预先定义的顺序,对状态的所有外部出口迁移进行检查,如果该迁移网络能成功到达某个状态则执行对应的迁移动作,并完成该轮事件广播;否则,按照预定的顺序对内部出口迁移进行检查,如果该迁移网络能成功到达某个状态则执行对应的迁移动作,并完成该轮事件广播;如果上述两项均不成功,则执行中间动作.然后把广播事件对其包含的子图进行广播.当状态迁移成功时,需要找到源状态和目的状态之间的共同路径,从内至外,执行源状态出口动作至共同路径,然后从外向内,执行目标状态的入口动作.图2描述了一个Stateflow的模型,该图中有两个并行的状态A与B,同时处于激活状态.状态A和状态B中分别包含六个和两个互斥状态,它们层次结构如图2所示.图2中的迁移边包含三种特殊情况:以实心黑点开头的迁移边为默认迁移边,表明该箭头所指状态为默认激活状态;状态S11到S2的迁移为层间迁移;状态S21到S22的迁移为迁移网络,该迁移网络包含一个结点.2.4混成CSP简介2.4.1混成CSPHCSP(Hybrid Communicating Sequential Process)[5,6]在CSP(Communicating Sequential Process)的基础上引入了微分方程以描述在混成系统中的连续动态行为,并同时引入各种中断机制(条件中断和通讯中断)来中断一个连续的微分行为,进而使连续微分行为和离散控制行为交错运行.HCSP能同时刻画系统的连续行为和控制逻辑,并具有CSP高可组合性的特性,是一种很好的用来描述大型控制系统的形式化建模语言.HCSP的语法如下所示:P::=skip|x:=e|ch?x|ch!e|P;Q|B→P|P⊔Q|P∗|⟨F(˙s,s)=0&B⟩|⟨F(˙s,s)=0&B⟩¤ i∈I(io i→Q i)S::=P|S∥S.†默认迁移是一种特殊的迁移,它仅有一个目标状态,表示该状态为默认激活状态.7郭丹青,吕继东,唐涛等:中国高速铁路列控系统的形式化分析与验证上式中P,Q,Q i,S均为HCSP进程,x和s为进程变量,ch为通道名,io i为通信事件(输入事件ch?x或输出事件ch!e),B和e为布尔表达式和算术表达式.以上语法结构的非形式化语义表示如下:skip不执行任何操作,立即终止;x:=e将表达式e的值赋给x;ch?x从通道ch中得到一个值,并将其赋给x;ch!e将e的值发送到通道ch;P;Q先执行进程P,当进程P终止时执行进程Q;B→P在B为真时执行P,否则立刻终止;P⊔Q由系统随机选择是执行进程P还是Q;P∗表示有限次重复执行进程P;⟨F(˙s,s)=0&B⟩代表微分动态过程,该微分过程相关的变量取值必须使得布尔表达式B始终取真值,否则该语句立即终止;⟨F(˙s,s)=0&B⟩¤ i∈I(io i→Q i)与⟨F(˙s,s)=0&B⟩执行过程基本一致,但是当通信事件io i发生时会立刻执行对应进程Q i,否则它会一直执行到微分过程终止为止;S1∥S2中S1和S2为并发进程,在不需要通讯时他们各自独立运行,当需要通讯时通讯双方需要同步执行该通讯过程,并发进程之间不能存在共享变量或是共享输入或者输出通道.2.4.2混成Hoare逻辑在[7]中,我们通过在经典的Hoare逻辑中引入历史公式构建了HCSP的证明系统.之后,[8]将该逻辑在Isabelle/HOL中实现,并利用该工具验证了中国铁路控制系统中的一个实际案例.历史公式由一个时段验算公式[19,20]表达,用于描述系统在执行时间区段内满足的性质.在混成Hoare逻辑(HHL)中,顺序进程P的性质描述由三元组{pre}P{post;HF}来表达,其中pre,post,和HF分别表示前置条件,后置条件,和历史公式.前置条件和后置条件由一阶逻辑公式表达,历史公式由时段验算公式表达.对于并发进程其性质描述如下表示.{pre1,...,pre n}P1∥···∥P n{post1,...,post n;HF1,...,HF n},其中pre,post i和HF i分别表示第i个进程的前置条件,后置条件和历史公式.i3场景描述及其建模过程我们将高速铁路列控系统在等级转换和模式转换的两个场景中涉及的主体抽象为由列控系统,列车与司机组成.列控系统分为多个等级,我们只考虑列车在C2级和C3级下面的行为,所以只涉及C2级列控系统与C3级列控系统,分别是TCC(Train Control Center),RBC(Radio Block Center),其主要作用是根据车载子系统,地面子系统等提供的列车状态,轨道占用,临时限速命令,联锁进路状态,灾害防护等信息产生其控制范围内各个列车的行车许可等控制信息,并且传输给车载设备以控制列车的运行.列车运行的过程中,司机要监控列车的运行,需要经常对不同的情况作出反应.图3是表示整个列控系统的Stateflow模型,其中C2级控制器,C3级控制器,列车和司机四个部分分别由由TCC,RBC,Train和Driver四个状态表示,这四个部分需要并行执行,所以这四个状态之间的关系为并行.列车在运行过程中需要一个动力系统来计算速度,由于动力系统是连续的,因此我们在Stateflow外使用一个子系统来表示,如图4所示,这个子系统的输入是加速度,输出是速度和距离,使用了两个积分模块,分别代表加速度的积分为速度,速度的积分为距离.运营场景是对运营中系统工作方式的简要描述,在这个模型中,我们考虑三个运营场景:行车许可（Movement Authority,MA）场景,控制器会对列车的运行过程进行监控,并且通过MA授权控制列车的行为;等级转换场景,列车在运行中可以从C2级转换到C3级,图5是等级升级场景的一个示意图;模式转换场景,列车在不同的模式下面也会有不同的行为模式,不同等级下面不同的模式之间的转换有不同的行为模式.下面我们会详细介绍这三个场景.因为从文档中推测等级转换和模式转换相结合8。

文章编号：1000-4750(2021)02-0187-11基于元胞自动机微观模拟的随机车流与桥梁耦合振动数值研究周军勇，苏建旭，齐 飒(广州大学土木工程学院，广东，广州 510006)摘 要：将经典车桥耦合振动理论与最新提出的多轴单元胞自动机(MSCA)微观车流荷载模拟方法进行融合，形成了一种精细化的随机车流与桥梁耦合振动数值分析方法。

介绍了该研究所采用的车桥耦合振动理论及模型；提出了MSCA 实现车桥动力分析的思路和方法，并进行了程序开发；通过具有实测时程动态挠度的工程算例，验证MSCA 实现车桥耦合动力分析的准确性；将MSCA 用于随机车流激励下某斜拉桥的动力效应分析中，论证基于MSCA 的随机车流与桥梁耦合振动分析程序的可靠性。

研究结果表明：工程算例很好地证明了该文所提方法和模型在进行车桥耦合分析的准确性，最大误差仅为11.6%；斜拉桥在随机车流作用下的静力与动力时程挠度分析显示，两者具有很好的一致性，随着路面粗糙度等级提升两者差异更加显著，说明了该模型和方法在开展随机车流与桥梁耦合振动分析的可靠性。

该研究进一步拓展了MSCA 在随机车流激励下分析桥梁各类动态响应的能力，为该方法程序在实桥监测与评估的应用提供了基础。

关键词：桥梁工程；车桥耦合；随机车流模拟；多轴单元胞自动机；数值分析中图分类号：U441+.2 文献标志码：A doi: 10.6052/j.issn.1000-4750.2020.04.0239NUMERICAL INVESTIGATION ON RANDOM TRAFFIC-BRIDGE COUPLED VIBRATION USING CELLULAR AUTOMATON-BASED MICROSCOPIC SIMULATIONZHOU Jun-yong , SU Jian-xu , QI Sa(College of Civil Engineering, Guangzhou University, Guangdong, Guangzhou 510006, China)Abstract: A numerical delicacy method for random traffic-bridge coupled vibration analysis is proposed.Incorporating the classical vehicle-bridge interaction theory, it is a newly established multi-axle single-cell cellular automaton (MSCA)-based microscopic traffic load simulation approach. The utilized equations and models in the classical vehicle-bridge interaction theory are introduced. The concepts and routes of the realization of MSCA for vehicle-bridge coupled dynamic analysis are proposed, and the relevant code program is developed.An engineering example with measured time-history dynamic deflections is utilized to verify the accuracy of the vehicle-bridge interaction analysis by MSCA. MSCA is used to analyze the dynamic load effects of a cable-stayed bridge under the excitation of random traffic loads, to demonstrate the reliability of the proposed approach. The results indicate that MSCA has good accuracy in vehicle-bridge coupling analysis. The maximum error in the engineering example is 11.6%. The static and dynamic time-history deflections of the cable-stayed bridge under random traffic loads show that they have good consistency, and the difference between them becomes more significant along with the increase in the pavement roughness grade. These prove the reliability of the proposed model and method in the random traffic-bridge coupled vibration analysis. This study forwards MSCA's ability to收稿日期：2020-04-19；修改日期：2020-07-29基金项目：国家自然科学基金项目(51808148)；广东省自然科学基金项目(2019A1515010701)；广州市科技计划项目(201904010188)通讯作者：周军勇(1990−)，男，江西人，讲师，博士，主要从事桥梁工程研究(E-mail: ***************.cn ).作者简介：苏建旭(1994−)，男，广东人，硕士生，主要从事桥梁工程研究(E-mail: ****************);齐　飒 (1994−)，女，河南人，硕士生，主要从事桥梁工程研究(E-mail: ***************).第 38 卷第 2 期Vol.38 No.2工 程 力 学2021年2 月Feb.2021ENGINEERING MECHANICS187analyze various types of dynamic load effects of bridges under the excitation of random traffic flow, which provides more applications of MSCA in monitoring and evaluation of real bridges.Key words: bridge engineering; vehicle-bridge interaction; random traffic simulation; multi-axle single-cell cellular automaton; numerical investigation车桥耦合振动特性是桥梁在移动车辆荷载作用下结构响应行为的重要表征，不仅可以揭示桥梁结构参数、力学行为和损伤特性[1]，还能反演移动车辆荷载特性[2]，是桥梁工程领域一直以来的研究热点[3 − 4]。