The Study of Quasi Monte Carlo in the Parallel Computation of Invariant Measures

合集下载

1.1.2VectorMagnitude

two parallel vectors is zero - and
4. (signed) magnitude is equal to
, where is the angle between the two vectors,
measured from to .
The schoolbook formula is
where is the number of columns in matrix . The same inner dimension condition applies as noted above: the number of columns in must equal the number of rows in . Matrix multiplication is:
We are multiplying column vectors of by the scalar elements of .
1.2.3 Multiplying a Matrix by a Matrix
The multiplication , so that
is equivalent to a side-by-side arrangement of column vectors
application in solving systems of linear equations such as
1.2.8 Trace The trace of a matrix is simply the sum of the diagonals:
1.2.9 Eigenvalues and Eigenvectors
where is the determinant of the sub-matrix formed by neglecting the 'th row and the 'th column. The formula is symmetric, in the sense that one could also target the 'th column:

Study of B - rho pi decays at Belle

a r X i v :h e p -e x /0207007v 1 1 J u l 2002BELLEBelle Prerpint 2002-18KEK Preprint 2002-59Study of B →ρπdecays at BelleBelle Collaboration A.Gordon u ,Y.Chao z ,K.Abe h ,K.Abe aq ,N.Abe at ,R.Abe ac ,T.Abe ar ,Byoung Sup Ahn o ,H.Aihara as ,M.Akatsu v ,Y.Asano ay ,T.Aso aw ,V.Aulchenko b ,T.Aushev ℓ,A.M.Bakich an ,Y.Ban ag ,A.Bay r ,I.Bedny b ,P.K.Behera az ,jak m ,A.Bondar b ,A.Bozek aa ,M.Braˇc ko t ,m ,T.E.Browder g ,B.C.K.Casey g ,M.-C.Chang z ,P.Chang z ,B.G.Cheon am ,R.Chistov ℓ,Y.Choi am ,Y.K.Choi am ,M.Danilov ℓ,L.Y.Dong j ,J.Dragic u ,A.Drutskoy ℓ,S.Eidelman b ,V.Eiges ℓ,Y.Enari v ,C.W.Everton u ,F.Fang g ,H.Fujii h ,C.Fukunaga au ,N.Gabyshev h ,A.Garmash b ,h ,T.Gershon h ,B.Golob s ,m ,R.Guo x ,J.Haba h ,T.Hara ae ,Y.Harada ac ,N.C.Hastings u ,H.Hayashii w ,M.Hazumi h ,E.M.Heenan u ,I.Higuchi ar ,T.Higuchi as ,L.Hinz r ,T.Hokuue v ,Y.Hoshi aq ,S.R.Hou z ,W.-S.Hou z ,S.-C.Hsu z ,H.-C.Huang z ,T.Igaki v ,Y.Igarashi h ,T.Iijima v ,K.Inami v ,A.Ishikawa v ,H.Ishino at ,R.Itoh h ,H.Iwasaki h ,Y.Iwasaki h ,H.K.Jang a ℓ,J.H.Kang bc ,J.S.Kang o ,N.Katayama h ,Y.Kawakami v ,N.Kawamura a ,T.Kawasaki ac ,H.Kichimi h ,D.W.Kim am ,Heejong Kim bc ,H.J.Kim bc ,H.O.Kim am ,Hyunwoo Kim o ,S.K.Kim a ℓ,T.H.Kim bc ,K.Kinoshita e ,S.Korpar t ,m ,P.Krokovny b ,R.Kulasiri e ,S.Kumar af ,A.Kuzmin b ,Y.-J.Kwon bc ,nge f ,ai ,G.Leder k ,S.H.Lee a ℓ,J.Li ak ,A.Limosani u ,D.Liventsevℓ,R.-S.Lu z,J.MacNaughton k,G.Majumder ao, F.Mandl k,D.Marlow ah,S.Matsumoto d,T.Matsumoto au,W.Mitaroﬀk,K.Miyabayashi w,Y.Miyabayashi v,H.Miyake ae,H.Miyata ac,G.R.Moloney u,T.Mori d,T.Nagamine ar,Y.Nagasaka i,T.Nakadaira as,E.Nakano ad, M.Nakao h,J.W.Nam am,Z.Natkaniec aa,K.Neichi aq, S.Nishida p,O.Nitoh av,S.Noguchi w,T.Nozaki h,S.Ogawa ap, T.Ohshima v,T.Okabe v,S.Okuno n,S.L.Olsen g,Y.Onuki ac, W.Ostrowicz aa,H.Ozaki h,P.Pakhlovℓ,H.Palka aa,C.W.Park o,H.Park q,L.S.Peak an,J.-P.Perroud r, M.Peters g,L.E.Piilonen ba,J.L.Rodriguez g,F.J.Ronga r, N.Root b,M.Rozanska aa,K.Rybicki aa,H.Sagawa h,S.Saitoh h,Y.Sakai h,M.Satapathy az,A.Satpathy h,e,O.Schneider r,S.Schrenk e,C.Schwanda h,k,S.Semenovℓ,K.Senyo v,R.Seuster g,M.E.Sevior u,H.Shibuya ap,V.Sidorov b,J.B.Singh af,S.Staniˇc ay,1,M.Stariˇc m,A.Sugi v, A.Sugiyama v,K.Sumisawa h,T.Sumiyoshi au,K.Suzuki h,S.Suzuki bb,S.Y.Suzuki h,T.Takahashi ad,F.Takasaki h, K.Tamai h,N.Tamura ac,J.Tanaka as,M.Tanaka h,G.N.Taylor u,Y.Teramoto ad,S.Tokuda v,S.N.Tovey u,T.Tsuboyama h,T.Tsukamoto h,S.Uehara h,K.Ueno z, Y.Unno c,S.Uno h,hiroda h,G.Varner g,K.E.Varvell an,C.C.Wang z,C.H.Wang y,J.G.Wang ba,M.-Z.Wang z,Y.Watanabe at,E.Won o,B.D.Yabsley ba,Y.Yamada h, A.Yamaguchi ar,Y.Yamashita ab,M.Yamauchi h,H.Yanai ac,P.Yeh z,Y.Yuan j,Y.Yusa ar,J.Zhang ay,Z.P.Zhang ak,Y.Zheng g,and D.ˇZontar aya Aomori University,Aomori,Japanb Budker Institute of Nuclear Physics,Novosibirsk,Russiac Chiba University,Chiba,Japand Chuo University,Tokyo,Japane University of Cincinnati,Cincinnati,OH,USAf University of Frankfurt,Frankfurt,Germanyg University of Hawaii,Honolulu,HI,USAh High Energy Accelerator Research Organization(KEK),Tsukuba,Japani Hiroshima Institute of Technology,Hiroshima,Japanj Institute of High Energy Physics,Chinese Academy of Sciences,Beijing,PRChinak Institute of High Energy Physics,Vienna,Austria ℓInstitute for Theoretical and Experimental Physics,Moscow,Russiam J.Stefan Institute,Ljubljana,Slovenian Kanagawa University,Yokohama,Japano Korea University,Seoul,South Koreap Kyoto University,Kyoto,Japanq Kyungpook National University,Taegu,South Korear Institut de Physique des Hautes´Energies,Universit´e de Lausanne,Lausanne,Switzerlands University of Ljubljana,Ljubljana,Sloveniat University of Maribor,Maribor,Sloveniau University of Melbourne,Victoria,Australiav Nagoya University,Nagoya,Japanw Nara Women’s University,Nara,Japanx National Kaohsiung Normal University,Kaohsiung,Taiwany National Lien-Ho Institute of Technology,Miao Li,Taiwanz National Taiwan University,Taipei,Taiwanaa H.Niewodniczanski Institute of Nuclear Physics,Krakow,Polandab Nihon Dental College,Niigata,Japanac Niigata University,Niigata,Japanad Osaka City University,Osaka,Japanae Osaka University,Osaka,Japanaf Panjab University,Chandigarh,Indiaag Peking University,Beijing,PR Chinaah Princeton University,Princeton,NJ,USAai RIKEN BNL Research Center,Brookhaven,NY,USAaj Saga University,Saga,Japanak University of Science and Technology of China,Hefei,PR ChinaaℓSeoul National University,Seoul,South Koreaam Sungkyunkwan University,Suwon,South Koreaan University of Sydney,Sydney,NSW,Australiaao Tata Institute of Fundamental Research,Bombay,Indiaap Toho University,Funabashi,Japanaq Tohoku Gakuin University,Tagajo,Japanar Tohoku University,Sendai,Japanas University of Tokyo,Tokyo,Japanat Tokyo Institute of Technology,Tokyo,Japanau Tokyo Metropolitan University,Tokyo,Japanav Tokyo University of Agriculture and Technology,Tokyo,Japanaw Toyama National College of Maritime Technology,Toyama,Japanay University of Tsukuba,Tsukuba,Japanaz Utkal University,Bhubaneswer,Indiaba Virginia Polytechnic Institute and State University,Blacksburg,VA,USAbb Yokkaichi University,Yokkaichi,Japanbc Yonsei University,Seoul,South KoreaB events collected with the Belle detector at KEKB.Thebranching fractions B(B+→ρ0π+)=(8.0+2.3+0.7−2.0−0.7)×10−6and B(B0→ρ±π∓)=(20.8+6.0+2.8−6.3−3.1)×10−6are obtained.In addition,a90%conﬁdence level upper limitof B(B0→ρ0π0)<5.3×10−6is reported.Key words:ρπ,branching fractionPACS:13.25.hw,14.40.Nd1on leave from Nova Gorica Polytechnic,Nova Gorica,Sloveniamodes are examined.Here and throughout the text,inclusion of charge con-jugate modes is implied and for the neutral decay,B0→ρ±π∓,the notation implies a sum over both the modes.The data sample used in this analysis was taken by the Belle detector[9]at KEKB[10],an asymmetric storage ring that collides8GeV electrons against3.5GeV positrons.This produces Υ(4S)mesons that decay into B0B pairs.The Belle detector is a general purpose spectrometer based on a1.5T su-perconducting solenoid magnet.Charged particle tracking is achieved with a three-layer double-sided silicon vertex detector(SVD)surrounded by a central drift chamber(CDC)that consists of50layers segmented into6axial and5 stereo super-layers.The CDC covers the polar angle range between17◦and 150◦in the laboratory frame,which corresponds to92%of the full centre of mass(CM)frame solid angle.Together with the SVD,a transverse momen-tum resolution of(σp t/p t)2=(0.0019p t)2+(0.0030)2is achieved,where p t is in GeV/c.Charged hadron identiﬁcation is provided by a combination of three devices: a system of1188aerogelˇCerenkov counters(ACC)covering the momentum range1–3.5GeV/c,a time-of-ﬂight scintillation counter system(TOF)for track momenta below1.5GeV/c,and dE/dx information from the CDC for particles with very low or high rmation from these three devices is combined to give the likelihood of a particle being a kaon,L K,or pion, Lπ.Kaon-pion separation is then accomplished based on the likelihood ratio Lπ/(Lπ+L K).Particles with a likelihood ratio greater than0.6are identiﬁed as pions.The pion identiﬁcation eﬃciencies are measured using a high momentum D∗+data sample,where D∗+→D0π+and D0→K−π+.With this pion selection criterion,the typical eﬃciency for identifying pions in the momentum region0.5GeV/c<p<4GeV/c is(88.5±0.1)%.By comparing the D∗+data sample with a Monte Carlo(MC)sample,the systematic error in the particle identiﬁcation(PID)is estimated to be1.4%for the mode with three charged tracks and0.9%for the modes with two.Surrounding the charged PID devices is an electromagnetic calorimeter(ECL) consisting of8736CsI(Tl)crystals with a typical cross-section of5.5×5.5cm2 at the front surface and16.2X0in depth.The ECL provides a photon energy resolution of(σE/E)2=0.0132+(0.0007/E)2+(0.008/E1/4)2,where E is in GeV.Electron identiﬁcation is achieved by using a combination of dE/dx measure-ments in the CDC,the response of the ACC and the position and shape of the electromagnetic shower from the ECL.Further information is obtained from the ratio of the total energy registered in the calorimeter to the particle momentum,E/p lab.Charged tracks are required to come from the interaction point and have transverse momenta above100MeV/c.Tracks consistent with being an elec-tron are rejected and the remaining tracks must satisfy the pion identiﬁcation requirement.The performance of the charged track reconstruction is studied using high momentumη→γγandη→π+π−π0decays.Based on the relative yields between data and MC,we assign a systematic error of2%to the single track reconstruction eﬃciency.Neutral pion candidates are detected with the ECL via their decayπ0→γγ. Theπ0mass resolution,which is asymmetric and varies slowly with theπ0 energy,averages toσ=4.9MeV/c2.The neutral pion candidates are selected fromγγpairs by requiring that their invariant mass to be within3σof the nominalπ0mass.To reduce combinatorial background,a selection criteria is applied to the pho-ton energies and theπ0momenta.Photons in the barrel region are required to have energies over50MeV,while a100MeV requirement is made for photons in the end-cap region.Theπ0candidates are required to have a momentum greater than200MeV/c in the laboratory frame.Forπ0s from BE2beam−p2B and the energy diﬀerence∆E=E B−E beam.Here, p B and E B are the momentum and energy of a B candidate in the CM frame and E beam is the CM beam energy.An incorrect mass hypothesis for a pion or kaon produces a shift of about46MeV in∆E,providing extra discrimination between these particles.The width of the M bc distributions is primarily due to the beam energy spread and is well modelled with a Gaussian of width 3.3MeV/c2for the modes with a neutral pion and2.7MeV/c2for the mode without.The∆E distribution is found to be asymmetric with a small tail on the lower side for the modes with aπ0.This is due toγinteractions withmaterial in front of the calorimeter and shower leakage out of the calorimeter. The∆E distribution can be well modelled with a Gaussian when no neutral particles are present.Events with5.2GeV/c2<M bc<5.3GeV/c2and|∆E|< 0.3GeV are selected for theﬁnal analysis.The dominant background comes from continuum e+e−→qB events and jet-like qi,j|p i||p j|P l(cosθij)i,k|p i||p k|,r l=),where L s and L qqD0π+ decays.By comparing the yields in data and MC after the likelihood ratiorequirement,the systematic errors are determined to be4%for the modes with aπ0and6%for the mode without.Theﬁnal variable used for continuum suppression is theρhelicity angle,θh, deﬁned as the angle between the direction of the decay pion from theρin the ρrest frame and theρin the B rest frame.The requirement of|cosθh|>0.3 is made independently of the likelihood ratio as it is eﬀective in suppressing the background from B decays as well as the qB events is used[14].The largest component of this background is found to come from decays of the type B→Dπ;when the D meson decays via D→π+π−,events can directly reach the signal region while the decay D→K−π+can reach the signal region with the kaon misidentiﬁed as a pion.Decays with J/ψandψ(2S) mesons can also populate the signal region if both the daughter leptons are misidentiﬁed as pions.These events are excluded by making requirements on the invariant mass of the intermediate particles:|M(π+π−)−M D0|>0.14 GeV/c2,|M(π+π0)−M D+|>0.05GeV/c2,|M(π+π−)−M J/ψ|>0.07GeV/c2 and|M(π+π−)−Mψ(2S)|>0.05GeV/c2.The widest cut is made around the D0mass to account for the mass shift due to misidentifying the kaons in D0 decays as pions.Fig.1shows the∆E and M bc distributions for the three modes analysed after all the selection criteria have been applied.The∆E and M bc plots are shown for events that lie within3σof the nominal M bc and∆E values,respectively. The signal yields are obtained by performing maximum likelihoodﬁts,each using a single signal function and one or more background functions.The signal functions are obtained from the MC and adjusted based on comparisons of B+→B0are assumed to be equal.The M bc distribution for all modes isﬁtted with a single Gaussian and an ARGUS background function[15].The normalization of the ARGUS function is left toﬂoat and shape of the function isﬁxed from the∆E sideband:−0.25 GeV<∆E<−0.08GeV and5.2GeV/c2<M bc<5.3GeV/c2.For the mode with only charged pions in theﬁnal state,the∆E distribution isﬁtted with a single Gaussian for the signal and a linear function withﬁxed shape for the continuum background.The normalization of the linear function is left to ﬂoat and the slope isﬁxed from the M bc sideband,5.2GeV/c2<M bc<5.26GeV/c2,|∆E|<0.3GeV.There are also other rare B decays that are expected to contaminate the∆E distribution.For the mode without aπ0,these modes are of the type B0→h+h−(where h denotes aπor K),B→ρρ(including all combinations of charged and neutralρmesons,where the polarizations of theρmesons are assumed to be longitudinal)and B→Kππ(including the decays B+→ρ0K+,B+→K∗0π+,B+→K∗0(1430)0π+,B+→f0(980)K+ and B+→f0(1370)K+)[16].These background modes are accounted for by using smoothed histograms whose shapes have been determined by combining MC distributions.The three B→ρρmodes are combined into one histogram. The normalization of this component is allowed toﬂoat in theﬁt due to the uncertainty in the branching fractions of the B→ρρmodes.Likewise,the B→hh and all the B→Kππmodes are combined to form one hh and one Kππcomponent.The normalizations of these components areﬁxed to their expected yields,which are calculated using eﬃciencies determined by MC and branching fractions measured by previous Belle analyses[16,17].The∆Eﬁts for the modes with aπ0in theﬁnal state have the signal compo-nent modelled by a Crystal Ball function[18]to account for the asymmetry in the∆E distribution.As for the B+→ρ0π+mode,the continuum background is modelled by a linear function withﬁxed slope.Unlike the B+→ρ0π+mode, a component is included for the background from the b→c transition.The pa-rameterization for rare B decays includes one component for the B→Kππ0 modes(B0→ρ+K−and B0→K∗+π−)[19]and one for all the B→ρρmodes.The normalization of the B→ρρcomponent is left toﬂoat while the other components from B decays areﬁxed to their expected yields.Table1summarizes the results of the∆Eﬁts,showing the number of events, signal yields,reconstruction eﬃciencies,statistical signiﬁcance and branching fractions or upper limits for eachﬁt.The statistical signiﬁcance is deﬁned assystematic error in theﬁtted signal yield is estimated by independently varying eachﬁxed parameter in theﬁt by1σ.Theﬁnal results are B(B+→ρ0π+)=(8.0+2.3+0.7−2.0−0.7)×10−6and B(B0→ρ±π∓)=(20.8+6.0+2.8−6.3−3.1)×10−6where theﬁrst error is statistical and the second is systematic.For theρ0π0mode,one standard deviation of the systematic error is added to the statistical limit to obtain a conservative upper limit at90%conﬁdence of5.3×10−6.The possibility of a nonresonant B→πππbackground is also examined.To check for this type of background,the M bc and∆E yields are determined for diﬀerentππinvariant mass bins.Byﬁtting the M bc distribution inππinvariant mass bins with B→ρπand B→πππMC distributions,the nonresonant contribution is found to be below4%.To account for this possible background, errors3.7%and3.2%are added in quadrature to the systematic errors of the ρ+π−andρ0π+modes,respectively.Theππinvariant mass distributions are shown in Fig.2.Two plots are shown for theρ+π−andρ0π+modes,one with events from the M bc sideband superimposed over the events from the signal region(upper)and one with events from signal MC superimposed over events from the signal region with the sideband subtracted(lower).Fig.3 shows the distribution of the helicity variable,cosθh,for the two modes with all selection criteria applied except the helicity condition.Events fromρπdecays are expected to follow a cos2θdistribution while nonresonant and other background decays have an approximately uniform distribution.The helicity plots are obtained byﬁtting the M bc distribution in eight helicity bins ranging from−1to1.The M bc yield is then plotted against the helicity bin for each mode and the expected MC signal distributions are superimposed.Both the ππmass spectrum and the helicity distributions provide evidence that the signal events are consistent with being fromρπdecays.The results obtained here can be used to calculate the ratio of branching frac-tions R=B(B0→ρ±π∓)/B(B+→ρ0π+),which gives R=2.6±1.0±0.4, where theﬁrst error is statistical and second is systematic.This is consistent with values obtained by CLEO[20]and BaBar[21,22]as shown in Table2. Theoretical calculations done at tree level assuming the factorization approx-imation for the hadronic matrix elements give R∼6[3].Calculations that include penguin contributions,oﬀ-shell B∗excited states or additionalππres-onances[4–8]might yield better agreement with the the measured value of R.In conclusion,statistically signiﬁcant signals have been observed in the B→ρπmodes using a31.9×106BWe wish to thank the KEKB accelerator group for the excellent operation of the KEKB accelerator.We acknowledge support from the Ministry of Ed-ucation,Culture,Sports,Science,and Technology of Japan and the Japan Society for the Promotion of Science;the Australian Research Council and the Australian Department of Industry,Science and Resources;the National Science Foundation of China under contract No.10175071;the Department of Science and Technology of India;the BK21program of the Ministry of Education of Korea and the CHEP SRC program of the Korea Science and Engineering Foundation;the Polish State Committee for Scientiﬁc Research under contract No.2P03B17017;the Ministry of Science and Technology of the Russian Federation;the Ministry of Education,Science and Sport of the Republic of Slovenia;the National Science Council and the Ministry of Education of Taiwan;and the U.S.Department of Energy.References[1] A.E.Snyder and H.R.Quinn,Phys.Rev.D48,2139(1993).[2]I.Bediaga,R.E.Blanco,C.G¨o bel,and R.M´e ndez-Galain,Phys.Rev.Lett.81,4067(1998).[3]M.Bauer,B.Stech,and M.Wirbel,Z.Phys.C34,103(1987).[4] A.Deandrea et al.,Phys.Rev.D62,036001(2000).[5]Y.H.Chen,H.Y.Cheng,B.Tseng and K.C.Yang,Phys.Rev.D60,094014(1999).[6] C.D.Lu and M.Z.Yang,Eur.Phys.J C23,275(2002).[7]J.Tandean and S.Gardner,SLAC-PUB-9199;hep-ph/0204147.[8]S.Gardner and Ulf-G.Meißner,Phys.Rev.D65,094004(2002).[9]Belle Collaboration,A.Abashian et al.,Nucl.Instr.and Meth.A479,117(2002).[10]E.Kikutani ed.,KEK Preprint2001-157(2001),to appear in Nucl.Instr.andMeth.A.[11]G.C.Fox and S.Wolfram,Phys.Rev.Lett.41,1581(1978).[12]This modiﬁcation of the Fox-Wolfram moments wasﬁrst proposed in a seriesof lectures on continuum suppression at KEK by Dr.R.Enomoto in May and June of1999.For a more detailed description see Belle Collaboration,K.Abe et al.,Phys.Lett.B511,151(2001).[13]CLEO Collaboration,D.M.Asner et al.,Phys.Rev.D53,1039(1996).[14]These MC events are generated with the CLEO group’s QQ program,see/public/CLEO/soft/QQ.The detector response is simulated using GEANT,R.Brun et al.,GEANT 3.21,CERN Report DD/EE/84-1,1984.[15]The ARGUS Collaboration,H.Albrecht et al.,Phys.Lett.B241,278(1990).[16]Belle Collaboration,A.Garmash et al.,Phys.Rev.D65,092005(2002).[17]Belle Collaboration,K.Abe et al.,Phys.Rev.Lett.87,101801(2001).[18]J.E.Gaiser et al.,Phys.Rev.D34,711(1986).[19]Belle Collaboration,K.Abe et al.,BELLE-CONF-0115,submitted as acontribution paper to the2001International Europhysics Conference on High Energy Physics(EPS-HEP2001).[20]CLEO Collaboration,C.P.Jessop et al.,Phys.Rev.Lett.85,2881(2000).[21]Babar Collaboration,B.Aubert et al.,submitted as a contribution paper tothe20th International Symposium on Lepton and Photon Interactions at High Energy(LP01);hep-ex/0107058.[22]BaBar Collaboration,B.Aubert et al.,submitted as a contribution paper tothe XXXth International Conference on High Energy Physics(ICHEP2000);hep-ex/0008058.Table1∆Eﬁt results.Shown for each mode are the number of events in theﬁt,the signal yield,the reconstruction eﬃciency,the branching fraction(B)or90%conﬁdence level upper limit(UL)and the statistical signiﬁcance of theﬁt.Theﬁrst error in the branching fraction is statistical,the second is systematic.ρ0π+15424.3+6.9−6.29.68.0+2.3+0.7−2.0−0.74.4σρ+π−30144.6+12.8−13.46.820.8+6.0+2.8−6.3−3.13.7σρ0π0116−4.4±8.58.5<5.3-Experiment B(B0→ρ±π∓)×10−6B(B+→ρ0π+)×10−6RE v e n t s /16 M e VE v e n t s /3 M e V /c2(b) ρ0π+Signal backgrd02.557.51012.51517.52022.55.25.225 5.25 5.2755.3E v e n t s /18 M e VE v e n t s /2 M e V /c2(d) ρ+π-Signal backgrd051015202530355.25.225 5.25 5.2755.3∆E(GeV)E v e n t s /18 M e V(e) ρ0π024681012-0.2-0.10.10.2(GeV/c 2)E v e n t s /2 M e V /c2M bc (f) ρ0πSignal backgrd02468101214165.25.225 5.25 5.2755.3Fig.1.The ∆E (left)and M bc (right)ﬁts to the three B →ρπmodes:ρ0π+,ρ+π−and ρ0π0.The histograms show the data,the solid lines show the total ﬁt and the dashed lines show the continuum component.In (a)the contribution from the B →ρρand B →hh modes is shown by the cross hatched component.In (c)the cross hatched component shows the contribution from the b →c transition and B →ρρmodes.102030405060+0(GeV/c 2)E v e n t s /0.1 G e V /c2M(π+π0)(GeV/c 2)E v e n t s /0.1 G e V /c2(GeV/c 2)E v e n t s /0.1 G e V /c2+-(GeV/c 2)E v e n t s /0.1 G e V /c2M(π+ π-)510152025Fig.2.The M (ππ)distributions for B 0→ρ±π∓(left)and B +→ρ0π+(right)events in the signal region.Plots (a)and (b)show sideband events superimposed;plots (c)and (d)show the sideband subtracted plots with signal MC superimposed.-1-0.500.51M b c y i e l d (E v e n t s )cos θh-1-0.500.51M b c y i e l d (E v e n t s )cos θhFig.3.The ρmeson helicity distributions for B 0→ρ±π∓(a)and B +→ρ0π+(b).Signal MC distributions are shown superimposed.。

油藏数值模拟基础

?活性剂驱?聚合物驱?混相驱思考题水平井模型第二章数学模型第一节数学模型的分类和推导原则一数学模型的分类按空间维数来分零维物质平衡方程一维岩心水驱油注采井间动态二维三维按流体相数来分单相气藏油藏弹性开发两相气藏水驱油藏水驱三相按流体组分来分单组分两组分按岩石类型来分单重介质砂岩双重介质碳酸盐岩低渗透油田按模型功能来分黑油模型凝析气藏模型双重介质模型热采模型基本模型聚合物驱模型三元复合驱模型水平井模型守恒关系运动方程状态方程辅助方程物质平衡关系能量平衡关系渗流darcy定律扩散fick定律导热fourier定律流体状态方程岩石状态方程流动辅助方程参数辅助方程化学辅助方程物理辅助方程偏微分方程组质量守恒方程组能量守恒方程质量守恒方程单位时间内流入单元体的流体质量单位时间内流出单元体的流体质量单位时间内从单元体注入或采出的流体质量单位时间内单元体中的质量增加量积分法2运动方程多相gradp为流体的屈服应力4含有启动压力低渗透岩石岩石压缩系数孔隙压缩系数气体理想气体实际气体nrtpv常数actualideal实际理想4
从离散的程度看，精度和速度是矛盾的。
三、用途
• 油藏描述
• 油藏动态预测
• 驱油机理研究
1. 油藏描述
油藏描述是油田开发的基础，是一项系统工程，由多学科各种方
法联合研究的结果。油藏数值模拟作为一种方法，在油藏描述中起了一
定的作用。－不同的方法研究的尺度不同
1) 孔隙结构研究 ~10μm级
CT 、核磁共振、图象分析仪、微观驱油机理、毛管压力实验
油藏数值模拟基础
华北油田培训班课程
中国石油大学石油工程学院 2008年9月
第一章油藏数值模拟进展
• 油藏数值模拟的基本概念
• 80年代的油藏数值模拟进展 • 90年代的油藏数值模拟进展

动态对等的最佳策略

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 simpliﬁcation 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 signiﬁcant loss of accuracy is to reduce the size of the system model by means of dynamic equiva-lents.The dynamic equivalent is a simpliﬁed 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 deﬁned 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)identiﬁcation 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 identiﬁcation approach is based on either parametric or non-parametric identiﬁcation [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 identiﬁcation 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 identiﬁcation 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].Artiﬁcial 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,ﬁctitious gen-erators are allocated.The external boundary is deﬁned 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 reductionAﬁrst 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 ﬂowchart.voltage calculation.In this paper,all nodes of the external system are eliminated,except the frontier nodes.By a load ﬂow study,the complex power that should inject some ﬁctitious 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 ﬂowing 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 ﬁve 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 ﬁve ﬁctitious 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 deﬁned,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 ﬁtness assignment strategies.To obtain diversiﬁcation 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 efﬁcient 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 ﬁtness 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 ﬁve 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 ﬂowchart 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 ﬁve 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 areﬁctitious 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 identiﬁcation of the equivalent gener-ators’parameters based on a speciﬁc 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 speciﬁc 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 signalﬁdelity measureis to compare two signals by providing a quantitative score thatdescribes the degree of similarity/ﬁdelity 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}aretwoﬁnite-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 deﬁned 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 speciﬁc 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 bemodiﬁed.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 deﬁned,as well as additional objective functions,as line’s power ﬂows,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,Identiﬁcation 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,Identiﬁcation 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,Identiﬁcation 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 lineﬂows 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,Identiﬁcationof dynamic voltage–current power system equivalents through artiﬁcial neural networks,Bulk power system dynamics and control–VI,2004.[19]A.M.Azmy,I.Erlich,Identiﬁcation 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 artiﬁcial neural networks,Electr.Power Energy Syst.24(2002)113–120.[22]A.M.Stankovic, A.T.Saric,osevic,Identiﬁcation of nonparametricdynamic power system equivalents with artiﬁcial 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 signalﬁdelity measures,IEEE Signal Process.Mag.26(1)(2009)98–117.。

数值积分和MonteCarlo方法数值积分令则零阶近似

第二章数值积分和Monte Carlo 方法第一节数值积分 ()ba S f x dx =⎰ 令 10,,k k n h x x x a x b +=-==，则()()()110(),''()k kn k k k k k k k k x x S f x dxf x f x x f f f x f f x +-=='=+-+==∑⎰()x fa k xb x零阶近似()()h f x f k O +=()()()∑∑-=-=O +=O +=110n k k n k k h f h h f h S一阶近似()()()21h hf f x x f x f kk k k O +--+=+ ∵()()⎰+=---=-++1212212122k kx x k k k kk k h x x x x x dx x x∴()()∑-=+O +⎪⎭⎫⎝⎛-+⋅=102121n k k k k h h f f h f S()()∑-=+O ++=102121n k k k h f f h 从直观看，用()112k k f f ++近似()f x 比只用k f 或1k f +好。

这方法也称Trapezoid 方法。

这样的数值积分方法的优点：● 简单直观，误差可以控制缺点：● “平均主义”，在()0≈x f 的区域，()k f x x ∆对S 贡献很小，但消耗同等的机时。

在多自由度系统这弱点尤为特出。

问题：直观地看，零级近似和一级近似的差别在哪？习题：编程序数值计算高斯积分。

第二节 Monte Carlo 方法如何用随机方法求积分？例如，可用‘抛石子’方法。

但这方法不比简单的数值积分有效。

1．简单抽样的Monte Carlo 方法均匀地随机地选取[b a ,]中{}k M x 个点，显然，(11()Mkk S f x M==+O ∑当M 足够大，当然可以得到足够好的积分值。

问题：为什么误差是(1/O ？答：不妨把这看成一个M 次测量的实验，假设每次测量都是独立的，由涨落理论，误差应为(1/O 。

一类动力学方程及流体力学方程解的Gevrey类正则性

Boltzmann 方程 . . . . . . . . . . . . . . . . . . . . . . . . 碰撞算子 Q(f, f ) 的基本性质 . . . . . . . . . . . . . . . . . Fokker-Planck 方程、Landau 方程以及 Boltzmann 方程线性化模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Navier-Stokes 方程 . . . . . . . . . . . . . . . . . . . . . . . Gevrey 函数空间 . . . . . . . . . . . . . . . . . . . . . . . .
研究现状及本文主要结果 . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 1.2.2 1.2.3 1.2.4 存在性及唯一性 . . . . . . . . . . . . . . . . . . . . . . . . . 动力学方程的正则性理论: 空间齐次情形 . . . . . . . . . . . 动力学方程的正则性理论: 空间非齐次情形 . . . . . . . . . . Navier-Stokes 方程的正则性理论 . . . . . . . . . . . . . . .
第二章预备知识 2.1 2.2 2.3 基本记号
Fourier 变换 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 基本函数空间及常用不等式 . . . . . . . . . . . . . . . . . . . . . . 2.3.1 2.3.2 Lp 空间及其性质 . . . . . . . . . . . . . . . . . . . . . . . . Sobolev 空间及其性质 . . . . . . . . . . . . . . . . . . . . .

propane

Quasiclassical trajectory calculations of collisional energy transfer in propane systemsApichart Linhananta and Kieran F.Lim*¤Centre for Chiral and Molecular T echnologies,School of Biological and Chemical Sciences,Deakin University,Geelong,V ictoria3217,Australia.E-mail:lim=.auRecei v ed6th December1999,Accepted27th January2000Published on the Web9th March2000Quasiclassical trajectory calculations of collisional energy transfer(CET)and rotational energy transfer from highly vibrationally excited propane to rare bath gases are reported.The calculations employed atomÈatom pairwise-additive Lennard-Jones,Buckingham exponential and hard-sphere intermolecular potentials to examine the dependence of CET on the intermolecular potential and to establish a protocol for future work on larger alkane systems.The role of the torsional(internal)and molecular(external)rotors in the energy-transfer mechanism were parison of the results with our earlier work on ethane]neon systems(A. Linhananta and K.F.Lim,Phys.Chem.Chem.Phys.,1999,1,3467)suggests that the internal and external rotors play a signiÐcant role in the deactivation mechanism for highly vibrationally excited alkanes.I.IntroductionGas-phase chemical reaction rates are strongly dependent on intermolecular collisional energy transfer(CET).CET is a vital component in any combustion-model and atmospheric-model systems.The only experimental CET quantities for hydrocarbon fuel molecules have been inferred““indirectlyÏÏfrom measurements of pressure-dependent reaction rates.1h3 Despite this,there have been no systematic theoretical dynamics studies of CET of hydrocarbon and halogenated hydrocarbons.In fact,most theoretical studies have been on small molecules.3h14The exceptions are the quasiclassical tra-jectory(QCT)calculations of azulene,toluene,benzene and hexaÑuorobenzene systems.15h23We have recently reported QCT calculations for highly vibrationally excited ethane in neon bath gas.24This and the recent work by Svedung et al. are,to our knowledge,theÐrst theoretical CET studies of an alkane with internal rotors.24,25Comparisons of theoretical and experimental studies of CET show that many of the dominant energy transfer mecha-nisms in small molecules are also present in large mol-ecules.3h6However,there are several di†erences between large-substrate and small-substrate behaviours.A notable example is that in the CET from a““smallÏÏsubstrate to a rare gas collider the trend He[Ne[Ar is observed,26h28 whereas the opposite trend of He\Ne\Ar is observed for ““largeÏÏsubstrates.29h39Theoretical studies of CET on small molecules employing various techniquesÈquantum,semi-classical and classical dynamicsÈall have correctly predicted the small-substrate behaviour.40,41This is not the case for large-substrate systems where QCT simulations incorrectly found the same smallsubstrate trend.15h18The discrepancy may be due to the lack of reliable data on the intermolecular potential surface involving large molecules and is most likely to be manifested in systems with the small collider helium bath gas.42,43¤Lim Pak Kwan.The aforementioned QCT calculations of large-substrate molecules have been on aromatic hydrocarbons because they have been most amenable to experimental studies using spec-troscopic probes.There have been fewer studies28,44h47on alkanes and branched-alkanes,which are the main com-ponents of common combustion fuels,and their halogenated analogues,which are important in ozone and““greenhouseÏÏchemistry.TheÐrst most obvious di†erences between alkanes and aromatics are their shapes,which are expected to a†ect the rotational energy transfer(RET).Since rotation to trans-lation(R]T)energy transfer and vibration to rotation (V]R)energy transfer are often more efficient than vibration to translation(V]T)energy transfer,this can have a strong inÑuence on the overall CET.Another crucial aspect is theÑexibility of alkanes.QCT simulations of alkanes would require the development of an efficient algorithm for sampling conformer space.Related to theÑexibility,as well as to RET,is the role of internal rotors in the CET mechanism.QCT calculations of highly vibra-tionally excited ethane in neon bath gas show that there is an interrelationship between the internal methyl rotors and the external rotation giving rise to V]torsion,R energyÑow in theÐrst collision,resulting in an““enhancedÏÏCET in sub-sequent collisions.24The same e†ect is also observed in experiments on the deactivation of highly vibrationally excited benzene and toluene,where toluene has much larger CET values than benzene.37This e†ect suggests that the torsional rotors in alkanes are important.Since the use of intramolecular torsional potential terms (nine such terms for each additional methylene unit)24plus a sampling of the conformational space may prove to be cost-prohibitive for large(r)alkane systems,there is a need to establish an e†ective protocol for QCT alkane simulations. Use of a hard-sphere potential will reduce a substrateÈcollider collision into a sequence of““sudden-impactÏÏatomÈatom encounters.Furthermore,there is no need to calculate molec-ular interactions at medium-to-large atomic separations.This paper““benchmarksÏÏCET using a hard-sphere potential against the more-commonly used Lennard-Jones andDOI:10.1039/a909614k Phys.Chem.Chem.Phys.,2000,2,1385È13921385This journal is The Owner Societies2000(Buckingham-type exponential-6models,by performing QCT calculations on the propane ]monatomic collider systems.The role of the torsional (internal)and molecular (external)rotors in the energy-transfer mechanism are reported.II.Quasiclassical trajectory calculationsA.Intermolecular potentialThe lack of knowledge of the detailed form of intermolecular potentials has always been a hindrance to quasiclassical mod-elling of CET.This is especially true for large-substrate systems,where there is a paucity of reliable theoretical and experimental data.Previous trajectory calculations of large molecules usually modelled the intermolecular potential by pairwise-additive atom Èatom potentials:7h 24,48h 50the inter-action parameters were usually obtained by semiempirical methods.Collins and coworkers have ““builtÏÏintermolecular potentials by interpolation of ab initio data:51h 53thus far they have only applied their method to relatively small polyato-mics whereas we wish to use a protocol that can be consistent-ly and easily ““scaled upÏÏfor larger alkane systems.Hence in this work three pairwise-additive atom Èatom intermolecular potentials were employed.The Ðrst intermolecular potential was the pairwise-additive Lennard-Jones (LJ)potential with atom Èatom terms given by V ij \4e ijCA p ij r ij B 12[A p ij r ijB 6D,(1)(i \C,H;j \rare gas),where is the atom Èatom centre-of-r ijmass separation,and and are the Lennard-Jones radiusp ij e ijand well depth,respectively.The LJ parameters were chosen by the method of Lim to match empirical values.16,29,54The second intermolecular potential was the pairwise-additive Buckingham exponential (exp-6)potential with atom Èatom terms given byV ij \A ij exp([c ij r ij )[C ij r ij~6,(2)where the parameter determines the repulsive steepness ofc ijthe potential.55The parameters and were chosen toA ij C ijmatch empirical values.16,29,54The last intermolecular potential was a pairwise-additive hard-sphere (HS)potentialV ij \GO ,0,r ij O r ij vdW ,r ij [r ijvdW ,(3)where is the van der Waals radius 56between atoms i andr ijvdW j .This potential is in the spirit of the e†ective mass theory.57The HS potential is tested here to determine if it can be used to derive useful qualitative information:if so then it would be a useful model for simulations of larger alkanes.The intermolecular parameters for propane ]Rg (Rg \rare gases He,Ne and Ar)potentials are given in Table 1.B.Intramolecular potentialA simple harmonic valence force Ðeld,consisting of harmonic stretches,bends and torsions,was used to describe the propane substrate:V intra\;i V stretch,i ];i V bend,i ];iV torsion,i .(4)The Ðrst two terms have been deÐned previously.15,58,59The harmonic stretching and bending force constants were obtained by the empirical prescription of Lindner:60k str,CC\4.705]102J m ~2,J m ~2,k str,CH \4.702]102k bend,CCH\6.67]10~17J rad ~2,and J rad ~2.k bend,HCH\5.61]10~17The Ðnal term in eqn.(4)is a 3-fold methyl torsional potential,which was assumed to be:V torsion,i \V 0n ;j /1n cos 2A 3qij 2B.(5)The torsional angles are the nine H ÈC ÈC ÈH or H ÈC ÈC ÈCq ijdihedral angles for each of the i th C ÈC bonds.Each carbon centre was assumed to have perfect tetrahedral geometry with C ÈC and C ÈH bond lengths of 0.1543nm and 0.1093nm,respectively.To study the e†ect of the torsion,the torsional barrier parameter was taken to be 0(free rotors)and 13.8V 0kJ mol ~1(experimental barriers).61The direction of the bond vectors was deÐned so that the staggered conformer has the lowest-energy geometry.The free-rotor model has apparent harmonic torsional ““vibrationalÏÏfrequencies of 9.2and 9.3cm ~1while the hindered-rotor model has apparent harmonic torsional ““vibrationalÏÏfrequencies of 167.4and 186.3cm ~1.These fre-quencies arise from the numerical normal mode analysis and are used in the selection of initial conditions.58,59,62The other 25vibrational frequencies compare favourably with experi-mental group frequencies of putational detailsTrajectory calculations were performed using program MARINER 58which is a customised version of VENUS96.59The LJ and exp-6potential models,selection of initial condi-tions,and general methodology are standard options in program MARINER/VENUS96.58,59,62The initial impact energy was chosen from a 300K thermal distribution.InE transthe majority of cases,the initial rotational angular momentum of propane was chosen from a thermal distribution at 300K.The rotational temperature was varied from 100to 1500K to investigate the RET of propane ]argon by the HS model.The initial vibrational phases and displacements were chosen from microcanonical ensembles at E @\41000,30000or 15000cm ~1,where E @is the rovibrational energy above the zero-point energy.These initial conditions are appropriate for comparison with the Ðrst few collisions in time-resolved infra-red Ñuorescence and ultraviolet absorption experi-ments.3h 6,29h 38,64Note that experiments measure the CET values of a cascade of collisions.The rovibrational energy dis-tribution of subsequent collisions will not be microcanonical,Table 1Intermolecular potential parametersLJ model exp-6model HS model p (e /k B )Aij Cij c ij r vdW /nm/K /kJ mol ~1/10~6kJ mol ~1nm 6/nm 1/nm H ÉÉÉHe 0.28258.0882294712479.945.50.325C ÉÉÉHe 0.291517.6931859254179745.60.345H ÉÉÉNe 0.293817.00103476168.5345.70.305C ÉÉÉNe 0.302034.156********.6945.90.325H ÉÉÉAr 0.306628.87140033519.240.80.335C ÉÉÉAr0.321658.025809650187641.00.3551386Phys .Chem .Chem .Phys .,2000,2,1385È1392but the CET behaviour of these subsequent collisions can be inferred 18,65,66from the microcanonical values.For the models employing the LJ and exp-6intermolecular potentials,trajectories were initialised with a centre-of-mass separation of 1.2nm and the classical equations of motion were integrated by the Adams ÈMoulton algorithm 58,59,62until the distance between the monatomic collider and the closest hydrogen exceeded a critical value of 1.0nm,at which point the trajectory was terminated.The initial impact param-eter b was chosen with importance sampling 16,17,58,59,62for values between 0nm and nm (He and Ne)or 0.9nmb m\0.8(Ar).These initial and Ðnal conditions were chosen by per-forming preliminary runs which showed that an insigniÐcant amount of energy was transferred at larger distances.For the HS interaction model,there is no intermolecular interaction until the point of impact,when the propane sub-strate is still described by a (near)microcanonical putationally,this is achieved by initialising trajectories as above,but translating the colliders to the point of initial contact without altering the rovibrational phases and orienta-tion.The translation was performed using an algorithm devel-oped by Alder and Wainwright 67,68to model hard-sphere Ñuid systems.After this initial point of contact,the trajectory was propagated normally.At each time step,the interatomic distances between the rare gas collider and every propane atom were checked for overlap.If an atom Èatom encounter occurred,the trajectory was projected back to the point of impact and the impulsive momentum transfer was calcu-lated.68The process was repeated until another encounter occurred or until the distance between the monatomic collider and the closest hydrogen exceeded a critical value,at which point the trajectory was terminated.Program MARINER 58was altered to implement the HS potential and trajectory-propagation algorithms.The short-ranged HS interaction per-mitted critical values as low as 0.4nm.Since the equations of motion are integrated for a comparatively short period,the HS model required much less computing time than the LJ and exp-6models.For E @\15000and 30000cm ~1,the integration time step was chosen to be 0.085fs,which is sufficient to conserve total energy to within 0.5cm ~1.This is approximately four times larger than the time step used in our previous ethane trajec-tory calculations.24Propane has less excitation per vibra-tional mode and hence energy can be conserved by larger time steps.For E @\41000cm ~1,it was necessary to employ a time step of 0.075fs to conserve energy.The numerical insta-bilities associated with the inversion of the methyl group(s)previously observed in simulations of ethane 24and toluene 16,17were not observed here.The calculations were performed on a DEC Alpha 3000/300LX workstation and an SGI Power Challenge Super-computer.In calculations that employed the LJ or exp-6intermolecular potentials,batches of 3000trajectories required approximately 60CPU hours for He collider and 100CPU hours for Ar on the workstation.The HS model decreased the required CPU time by a factor of 10:this reduction will be very signiÐcant in the study of larger alkanes.CPU time was reduced by a factor of about 4on the supercomputer.D.Rotation energy and torsional angular momentum It is well documented that rotational energy transfer is an effi-cient pathway for CET.3,24,65,66,69However,while angular momenta are well-deÐned,rovibrational coupling gives rise to an ambiguity in the deÐnition of rotational energy.Previous quasiclassical simulations employed several di†erent methods to decouple the rotational and vibrational energies.One method 11deÐnes the rotational energy asE rot \1(JI ~1J ),(6)where I and J are,respectively,the instantaneous moment of inertia and angular momentum.In a second method isE rotapproximated by the instantaneous angular momentum,but the moment of inertia is taken to be the equilibrium geometry value.11Both deÐnitions give rotational energies that oscillate with time.There is an alternative deÐnition that is valid for symmetrical top rotors:65E rot \1B effJ 2,(7)where J is the magnitude of the rotational angular momentum and is an e†ective rotational constant.This deÐnitionB effdecoupled the rovibrational energy so that the rotational energy includes only the ““adiabatic partÏÏ,whereas the ““activeÏÏpart is included with the vibrational energyE V \E [E rot,(8)where and E are,respectively,the vibrational and totalE Vinternal energies.Eqn.(7)is a valid approximation for sym-metrical top molecules.70The main advantage of this deÐni-tion is that,classically,it is a conserved quantity.The equilibrium Cartesian principal moments of inertia of propane are kg m 2,kgI xx \1.11]10~45I yy\9.7]10~46m 2and kg m 2.Hence,propane is a goodI zz\2.97]10~46approximation of a symmetrical top and it is possible to deÐne the rotational energy by eqn.(7),with the approx-imationB eff \12hc (I xx I yy I zz)~1@3.(9)It was shown in our previous work on ethane 24that the coupling between external and internal rotors enhances the overall CET.Hence the torsional angular momentum of propane was also monitored in this work.Whereas ethane has only one torsional rotor which lies along its molecular axis,propane has two distinct and unparallel torsional rotors.The deÐnition of the torsional angular momentum introduced for ethane is generalised by calculating the rotational angular momentum of the methyl group and the associated ethyl groupJ methyl \;i /H,H,Hr i ]p iJ ethyl\;i /C,H,Hr i ]p i,(10)where is the angular momentum of the methyl groupJ methyland is the angular momentum of the associated ethylJ ethylrotor.Note that for consistency with eqn.(5),only the six atoms directly bonded to each torsional C ÈC bond have been included in the summation in eqn.(10).The torsional angular momentum is then deÐned asJ tor \o (J methyl [J ethyl)Éa o ,(11)where is a unit vector parallel to the CC torsional axis.The a CET to/from the torsional rotors was monitored by calcu-lating the average torsional angular momentum change*J tor \J tor (Ðnal)[J tor(initial).(12)E.Data analysisTrajectory data were analysed by a bootstrap algorithm 71,72in program PEERAN.16,73Some 3000È5000trajectories were performed for each potential model.This was sufficient to obtain average energy-transfer quantities with statistical uncertainties of about 10%.However,the uncertainties for the average rotational energy transfer were about 20%,due to the initial rotational-energy Boltzmann distribution (rather than an initial microcanonical distribution).Trajectory averagesPhys .Chem .Chem .Phys .,2000,2,1385È13921387deÐned by (for both overall CET and RET)S*E n TtrajS*E n T traj \1N ;i /1N bi bm(*E i )n(13)are related to experimentally obtained quantities S*E n T by ratio of collision cross-sectionsS*E n T \p b m 2p p LJ2X (2,2)RS*E n T traj (14)where is the LJ collision cross-section and is thep LJ 2X (2,2)R b mmaximum impact parameter in the trajectory simulation.This normalisation removes the ambiguity related to the elastic scattering at high impact parameter.74The input LJ param-eters were obtained from ref.29.At 300K,the LJ collision cross-section values of nm 2,0.4834nm 2p LJ2X (2,2)R \0.3976and 0.6945nm 2for propane ]He,propane ]Ne and propane ]Ar,respectively,were obtained using the program COLRATE.75This corresponds to the LJ collision frequencies of m 3s ~1,328.58]10~18m 3s ~1Z LJ,coll\523.29]10~18and 382.37]10~18m 3s ~1,respectively.In this paper,we have reported both the Ðrst and second moments of the trajectory data since the Ðrst moment is usually more useful for comparison with experiment,but the QCT second moment is statistically more reliable.74Some experiments can determine both the Ðrst and second moments of the CET probability.3,5III.Results and discussionA.The e†ect of the torsional barrierFigs.1and 2show the CET values,[S*E T and S*E 2T 1@2,and the RET values,as functions of energy E @aboveS*E RT ,zero-point energy for propane ]neon.One set of results areFig.1Dependence of energy-transfer quantities on torsional barrier for deactivation of vibrationally excited propane by neon bath gas:)Hindered-rotor (LJ);Free-rotor (LJ);Hindered-rotor (exp-6);L +…Free-rotor (exp-6).Fig.2Dependence of rotational energy transfer on torsional barrier for deactivation of vibrationally excited propane by neon bath gas:)Hindered-rotor (LJ);Free-rotor (LJ);Hindered-rotor (exp-6);L +…Free-rotor (exp-6).for the free-rotor model the other for the hindered-(V 0\0),rotor model kJ mol ~1).These results are for the LJ(V 0\13.8and exp-6intermolecular potentials.The overall deactivation,[S*E T and S*E 2T 1@2,is larger for the hindered-rotor model,similar to results for ethane ]neon.24The torsional angular momentum transfer is shownS*J torT in Fig.3.Note that for the hindered-rotor modelsS*J torT with both LJ and exp-6intermolecular potentials are virtually identical:the reason for this is unclear.Overall,S*J torTdecreases,but remains positive,with the presence of a barrier In contrast,for ethane ]neon changes from posi-V 0.S*J torT tive to negative over a similar range of values.24This di†er-V 0ence is probably due to the higher torsional moment of inertia for propane torsion compared to ethane(CH 3ÈCCH 2)This means that propane torsion has higher(CH 3ÈCH 3).density of states and can more readily gain torsional excita-tion than ethane torsion,explaining why is positiveS*J torT for propane,but negative for ethane.In ethane,the torsion acts like a vibration providing an efficient torsion ]T pathway.24The increase in [S*E T and S*E 2T 1@2(Fig.1)for the hindered-rotor model suggests that propane torsions play the same role in the CET mechanism.The RET is smaller for the propane free-rotor modelS*E RT than the hindered-rotor model (Fig.2),contrary to the ethane results.24For ethane,the torsion is aligned along the molecu-lar axis,hence any increase in methyl-rotor angular momen-tum contributes to both (internal)torsional excitation S*J torTand (external)rotational excitation The propane free-S*E RT .rotor model has Ðve (three external and two internal)indepen-Fig.3Dependence of torsional angular momentum change on tor-sional barrier for deactivation of vibrationally excited propane by neon bath gas:Hindered-rotor (LJ);Free-rotor (LJ);)L +Hindered-rotor (exp-6);Free-rotor (exp-6).Note that the two sets …of hindered-rotor results are almost identical.1388Phys .Chem .Chem .Phys .,2000,2,1385È1392dent rotors,none of which have coincident axes.The extra rotors mean that there is less energy available to the external rotors in any V ]torsion,R energy redistribution.Noteworthy is the fact that the di†erences between the free-rotor and hindered-rotor models persist up to E @\41000cm ~1.For ethane ]neon,there is an onset of near-free-rotor behaviour at E @\30000cm ~1:at E @\41000cm ~1there is no signiÐcant di†erence between the free-and hindered-rotor models.However,the larger number of vibrational modes in propane,which decreases the excitation per torsional mode,ensures that the di†erences remain even at very high excita-tion.Hence correct theoretical treatments of internal rotors become even more essential for larger molecules.B.Trajectory results for LJ and exp-6modelsThe CET results for the deactivation of highly excited propane by helium,neon and argon are shown in Fig.4,where the intermolecular interactions have been modelled by the LJ and exp-6potentials.Three important features are:(1)Energy transfer increases with increasing E @and is in accord with theoretical and experimental studies on the deac-tivation of highly vibrationally excited molecules.(2)The LJ potential results in larger CET values than the exp-6model,since the LJ potential has a much harder repul-sive part than the exp-6potential.There are numerous works which concluded that CET depends mainly on the repulsive part of the intermolecular potential and that,in general,a harder repulsive part results in larger energy transfers.9,16,17(3)The deactivator efficiency shows the trend He [Ne [Ar which,unfortunately,is in discord with experi-mental trends for Ñuorinated alkane systems.28To our knowledge,there has been no experimental study of CET in propane ]rare gas systems.““IndirectÏÏstudies of related systems include 2-bromopropane ]Ne ([S*E T \130cm ~1for E @\17000È21000cm ~1)76andFig.4Energy-transfer quantities for deactivation of vibrationally excited propane by rare gases:Helium (LJ);Neon (LJ);)K |Argon (LJ);Helium (exp-6);Neon (exp-6);Argon (exp-6).+=>isotopically-substituted cyclopropane ]He (S*E 2T 1@2\200È400cm ~1for E @D 22000cm ~1).2These CET quantities were not directly measured,but were inferred from pressure-dependent thermal reaction rates at elevated temperatures.Some more recent studies using time-resolved optoacoustic spectroscopy include ([S*E T \114cm ~1atC 3F 8]Ar E @\15000cm ~1and [S*E T \300cm ~1at E @\40000cm ~1).46These studies reveal no information about RET nor the role of torsional modes.These experimental CET quan-tities correlate well with our present calculations (Fig.4)but also indicate a need for fresh experimental studies.The decreasing trend with collider He [Ne [Ar has been observed in many other QCT studies.9,15,18,77Although the lack of qualitative agreement with experiment is disappoint-ing,these studies and the present work have used very crude intermolecular potential models.Given the lack of detailed information about polyatomic intermolecular potential sur-faces,the intention in the present and other studies has been to use a set of consistent and transferable potentials,16much in the spirit of molecular mechanics force Ðelds.Experience with simulations on other systems would suggest that the exp-6model potentials predict ““betterÏÏCET values than the LJ potentials.17Fig.5plots the RET of propane ]rare gas systems.For Ne and Ar,monotonically increases with E @,whereas forS*E RT He,it initially increases but decreases at higher excitation energy.In all cases,RET is larger for the LJ model which is in accord with the CET behaviour.Clary and Kroes 78and others 16,17,40have observed that RET is larger for heavier colliders because the collision duration is closer to the rota-tional period of the molecular substrate.Fig.6plots the torsional angular momentum transfer as a function of E @.is largest for He and smal-S*J tor T S*J torT lest for Ar,which is the same trend as for CET.This implies that,in addition to the external rotor gateway,the torsional rotor is a gateway for facile CET via V,torsion ]torsion,T.24An interesting feature of Fig.6is that seems to beS*J torT Fig.5Rotational energy transfer for deactivation of vibrationally excited propane by rare gases:Helium (LJ);Neon (LJ);)K |Argon (LJ);Helium (exp-6);Neon (exp-6);Argon (exp-6).+=>Phys .Chem .Chem .Phys .,2000,2,1385È13921389Fig.6Torsional angular momentum change for deactivation ofvibrationally excited propane by rare gases:Helium(LJ);Neon)K(LJ);Argon(LJ);Helium(exp-6);Neon(exp-6);Argon|+=>(exp-6).insensitive to the intermolecular potential.However,the factthat it depends on the type of bath gas indicates a dependenceon the mass of the deactivator.This suggests that isS*JtorTinsensitive to theÐne details of the intermolecular potentialand can be modelled by either LJ or exp-6potentials.C.Trajectory results for hard-sphere modelLJ and exp-6potentials have long-range attractive terms andare computationally expensive.Since HS is a short-rangepotential,it is computationally cheaper in terms of computertime than other potential models by an order of magnitude.Inthis section we compare the results of the short-range HS withthe longer-range potentials.Fig.7shows S*E T and S*E2T1@2for the HS model.Fig.8shows the RET for the HS model.The qualitative behavioursare the same as for the LJ and exp-6models but the energy-transfer values are several times larger than for the LJ andexp-6model.This is not surprising in view of the““hardnessÏÏof the HS potential.9,16,17Another important feature is thatS*E T and S*E2T1@2for He are several times larger than forNe and Ar.This is also true for the LJ model(Fig.3)whichindicates that the HS and LJ models tend to give CET valuesthat are much too high for helium colliders.Table2lists the average number of encounters per collision,for He,Ne and Ar colliders.This average includes onlyNC,trajectories in which collisions have occurred.As expected NCTable2Average number of atomÈatom encounters NcE@/cm~1150003000041000Propane]He 1.967 1.884 1.847Propane]Ne 3.145 2.952 2.852Propane]Ar 3.753 3.501 3.400Fig.7Energy-transfer quantities for deactivation of vibrationallyexcited propane by rare gases for the HS model:Helium;Neon;+=Argon.>is largest for Ar and smallest for He due to their reducedmasses.also decreases with increasing E@which suggestsNCthat a more highly excited substrate imparts more energy perencounter to the deactivator,reducing the collision duration.Fig.9shows S*E T,and for propane]argonS*EVT S*ERTsystems at rotational temperatures300,1000andTROT\100,1500K.In these simulations,initial excitation wasÐxed atE@\15000cm~1and the initial translational temperaturewas K.It can also be seen that RETTtrans\300S*ERTdecreases with increasing the magnitude of the vibra-TROT;tional energy transfer also decreases with increasingS*EVTThis implies that rotationally cold systems exhibitTROT.V]R,T energy transfer,whereas rotationally hot systemsexhibit V,R]R,T.It can be seen that the overall[S*E T islarger for larger which agrees with the hypothesis thatTROTthe external rotation is a facile CET path.This behaviour hasFig.8Rotational energy transfer for deactivation of vibrationallyexcited propane by rare gases for the HS model.Helium;Neon;+=Argon.>1390Phys.Chem.Chem.Phys.,2000,2,1385È1392。

系统辨识_6_多新息辨识理论与方法_丁锋

的最小二乘辨识算法或随机梯度等辨识算法有下列形式： ^ （ t） = θ ^ （ t － 1 ） + L （ t ） e（ t ）， θ e（ t）： = 其中 L（ t） ∈R 为算法增益向量（ gain vector）， T ^ （ t － 1 ） ∈R 为标量新息（ scalar innovay （ t）－ φ （ t ） θ tion），即单新息（ single innovation）．这个算法可以这样描述： t 时刻的参数估计向量 ^ （ t）是用增益向量 L （ t）与标量新息 e （ t ）的乘积， θ 对 ^ （ t － 1 ）进行修正， ^ （ t） t － 1 时刻参数估计向量 θ 即θ ^ （ t － 1 ）的基础上加上增益向量 L （ t ）与新息是在 θ e（ t）的乘积．这种方法也称为新息修正辨识方法或新息辨识方法．上述算法中新息 e （ t ）是标量，我们把这个标量 in新息加以推广，就导出了多新息辨识方法（ multinovation identification method ）［24］．多新息辨识理论（ multiinnovation identification theory ）就是将单新息从新息修正角度提出多新息修修正技术加以推广，正技术辨识的概念，建立多新息修正辨识方法，简称多新息辨识方法．顾名思义，多新息算法就是将新息加以推广．对将算法中的标量新息 e （ t ） ∈ R 推广标量系统而言， t ） ∈ Rp ， innova为新息向量 E （ p，即多新息（ multin tion），为使矩阵乘法维数兼容，增益向量 L （ t ） ∈ R t ） ∈R n × p ，须推广为增益矩阵（ gain matrix） Γ（ p，那么 n

张平文

张平文北京大学数学科学学院科学与工程计算系邮编：100871 电话：86‐10‐6275‐9851 传真：86‐10‐6275‐1801电子邮件：pzhang@ 个人网页：/pzhang教育背景1988 ‐ 1992 博士研究生，北京大学数学科学学院导师：应隆安教授1984 ‐ 1988 学士，北京大学数学系工作经历2008‐ 副院长，北京大学数学科学学院2001‐ 常务副主任，北京大学科学与工程计算中心1999‐ 系主任，北京大学数学科学学院科学与工程计算系1996‐ 教授，北京大学数学科学学院1994 ‐ 1996 副教授，北京大学数学科学学院1992 ‐ 1994 讲师，北京大学数学科学学院研究领域复杂流体的多尺度建模与计算多尺度分析与计算移动网格方法及其应用所获荣誉与奖励2007 教育部创新团队“偏微分方程数值计算及其应用”带头人2007 与汤涛、李若和汤华中合作的项目“基于调和映射的移动网格方法及其应用”获得教育部高校科学技术奖自然科学一等奖2002 国家教育部全国普通高等院校优秀教材二等奖2002 国家杰出青年基金2002 长江学者2002 北京市第十六届“五四青年奖章”2000 教育部首届高等学校优秀青年教师教学科研奖励计划‘青年教师奖’1999 霍英东教育基金会第七届高等院校青年教师奖（研究类）一等奖1999 冯康科学计算奖1995 与应隆安合著“涡度法”获全国优秀科技图书二等奖1994 北京计算数学学会优秀论文奖1992 光华奖学金1990 九章数学奖学金学术兼职2006 ‐ 学术委员会副主任，北京应用物理与计算数学研究所计算物理实验室2004 ‐ 副理事长，中国工业与应用数学学会(CSIAM)2004 ‐ 学术委员会主任，中国工业与应用数学学会(CSIAM)2006 ‐ 学术委员会委员，“大规模科学与工程计算”国家重点实验室2002 ‐ 2006 副理事长，中国计算数学学会2001 ‐ 2006 学术委员会副主任，“大规模科学与工程计算”国家重点实验室2005 ‐ 兼职教授，吉林大学2004 ‐ 兼职教授，湘潭大学2004 ‐ 兼职教授，苏州大学1999 ‐ 2001 兼职教授，清华大学学术交流2005‐2006 访问学者，Department of Mathematics, Hong Kong University of Science and Technology (3个月)2005 访问学者，Applied and Computational Mathematics, Princeton University (2个月)2004 访问学者，Applied and Computational Mathematics, Princeton University (3个月)2002 访问学者，Applied and Computational Mathematics, Princeton University (2 个月)2001 访问学者，Department of Mathematics, Hong Kong University of Science and Technology (2个月) 1999 访问学者，Department of Applied Mathematics, California Institute of Technology (2个月) 1999 访问学者，Department of Mathematics, Hong Kong Baptist University (2个月)1998‐1999 访问学者，Department of Mathematics, Hong Kong Baptist University (6个月)1997 访问学者，Department of Applied Mathematics, California Institute of Technology (10个月) 1995 访问学者，Department of Applied Mathematics, California Institute of Technology (11个月)杂志编委2007 Journal of Computational Mathematics (JCM)2006 Communications in Computational Physics (CiCP)2006 International Journal of Nonlinear Science(IJNS)2005 Communications in Mathematical Sciences (CMS)2005 SIAM Journal on Numerical Analysis (SINUM)2005 Journal of Information and Computational Science(JOICS)2002 Applied Mathematical Research Express (AMRX);2007 《工程数学学报》2006 《数学杂志》2004 《计算数学》2004 《计算物理》2004 《东北数学》学术会议会议名称时间地点第六届流体力学计算方法 1993.8 济南第三次中日计算数学会议 1996.8 大连第一界中意应用数学与计算数学会议 1996.8 北京力学中偏微分方程与数值计算 1999.5 香港科学计算会议 1999.6 香港海内外青年计算数学工作者会议 1999.7 北京第九届流体力学计算方法 1999.9 温州第六届全国计算数学年会 1999.10 上海中瑞偏微分方程数值解研讨会 2000.1 香港第二界中意应用数学与计算数学会议 2000.6 意大利第一届中韩数值分析及其应用研讨会 2001.2 汉城科学与工程计算国际研讨会 2001.3 北京计算与应用偏微分方程国际讨论会 2001.7 张家界2002.6 瑞典 Third China‐Sweden Workshop on ComputationalMathematics， GotebergWorkshop on Multiscale Analysis and Computation 2002.6 台湾2002.8 上海 11th International Conference of Fluid Dynamics and SoftCondensed Matter2002.8 西安 ICM2002‐Beijing Satelite Conference on ScientificComputing科学计算与应用国际研讨会 2003.1 香港中韩数值分析及其应用国际研讨会 2003.2 北京2003.11 意大利 Third China‐Italy Joint Conference on Computational andApplied Mathematics中国数学会第九次年会 2003.11 武汉 International Conference of Scientific Computing 2003.12 北京非线性动力系统和随机偏微分方程国际研讨会 2004.5 北京 International Conference on Superconvergence and A 2004.5 湖南长沙Posteriori Estimates in FEM2004.5 上海 Symposium on Multiscale Computation in Materials andBiologyInternational Workshop on Wave Propagations 2004.6 北京 The 2nd International Conference on Inverse Problem 2004.6 上海2004.6 北京 International Conference on Frontiers of AppliedMathematicsSummer School on Multiscale Modelings and Simulations2004.6 北京 International Conference on Numerical and Applied PDEs2004.6 吉林长春 Workshop on Multiscale Rheological Models for Fluids 2004.11 Montreal， Canada2005.1 Singapore Nanoscale Material Interfaces: Experiment, Theory andSimulationThe 3rd joint Chinese‐Korean Workshop on Recent2005.2 South Korea Progresses on Numerical Analysis and Its Applications2005.6 BeijingThe Fourth China‐Italy Conference: Mathematical Modelsin Life Science‐Theory and Simulation.The 1st China‐Germany Workshop on Computational and2005.9 Berlin,German Applied MathematicsThe second International Conference on Scientific2005.12 香港 Computing and Partial Differential Equations2006.5 北京 International Conference on Calculus of Variations, PDEsand Nonlinear AnalysisInterfacial Dynamics in Complex Fluids2006.5 Banff,Canada International Symposium on Polymer Physics2006,6 Suzhou2006 International Conferences on Applied Mathematics2006.6 天津 and Interdisciplinary Research2006.6 Beijing International Conference on Recent Advances in ScientificComputationsWorkshop on Multiscale Modeling of Complex Fluids2006.6 Beijing International Conference on PDE and Numerical Analysis2006.6 长沙2006.10 西安The Symposium on Multi‐physics and Muti‐ScaleComputation of Materials‐20062006.12 Singapore International Workshop on Multiscale Analysis andApplicationsMultiscale Modeling of Complex Fluids 2007.4 Maryland Canada‐China workshop on industrial mathematics 2007.8 Banff, Canada2008.6 Busan, Korea6th International Conference on Scientific Computing andapplicationsWorkshop on Nanoscale Interfacial Phenomena in Complex2008.6 BeijingFluids另外，1999年6月作为政府代表团成员参加由联合国教科文组织在匈牙利举办的世界科学家大会及青年科学家论坛。

【5A文】关于序列二次规划(SQP)算法求解非线性规划问题研究

关于序列二次规划（SQP）算法求解非线性规划问题研究兰州大学硕士学位论文关于序列二次规划（SQP）算法求解非线性规划问题的研究姓名：石国春申请学位级别：硕士专业：数学、运筹学与控制论指导教师：王海明20090602兰州大学2009届硕士学位论文摘要非线性约束优化问题是最一般形式的非线性规划NLP问题，近年来，人们通过对它的研究，提出了解决此类问题的许多方法，如罚函数法，可行方向法，Quadratic及序列二次规划SequentialProgramming简写为SOP方法。

本文主要研究用序列二次规划SOP算法求解不等式约束的非线性规划问题。

SOP算法求解非线性约束优化问题主要通过求解一系列二次规划子问题来实现。

本文基于对大规模约束优化问题的讨论，研究了积极约束集上的SOP 算法。

我们在约束优化问题的s一积极约束集上构造一个二次规划子问题，通过对该二次规划子问题求解，获得一个搜索方向。

利用一般的价值罚函数进行线搜索，得到改进的迭代点。

本文证明了这个算法在一定的条件下是全局收敛的。

关键字：非线性规划，序列二次规划，积极约束集Hl兰州人学2009届硕二t学位论文AbstractNonlinearconstrainedarethemostinoptimizationproblemsgenericsubjectsmathematicalnewmethodsareachievedtosolveprogramming．Recently,Manyasdirectionit，suchfunction，feasiblemethod，sequentialquadraticpenaltyprogramming??forconstrainedInthisthemethodspaper,westudysolvinginequalityabyprogrammingalgorithm．optimizationproblemssequentialquadraticmethodaofSQPgeneratesquadraticprogrammingQPsequencemotivationforthisworkisfromtheofsubproblems．OuroriginatedapplicationsinanactivesetSQPandSQPsolvinglarge-scaleproblems．wepresentstudyforconstrainedestablishontheQPalgorithminequalityoptimization．wesubproblemsactivesetofthesearchdirectionisachievedQPoriginalproblem．AbysolvingandExactfunctionsaslinesearchfunctionsubproblems．wepresentgeneralpenaltyunderobtainabetteriterate．theofourisestablishedglobalconvergencealgorithmsuitableconditions．Keywords：nonlinearprogramming，sequentialquadraticprogrammingalgorithm，activesetlv兰州大学2009届硕士学位论文原创性声明本人郑重声明：本人所呈交的学位论文，是在导师的指导下独立进行研究所取得的成果。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。