The Boson Normal Ordering Problem and Generalized Bell Numbers

a rXiv:h ep-ph/999571v13Se p1999DOE/ER/40762-170U.ofMd.PP#99-054NT@UW-99-53High Energy Theorems at Large-N Silas R.Beane Department of Physics,University of Maryland College Park,MD 20742-4111and Department of Physics,University of Washington Seattle,WA 98195-1560sbeane@ ABSTRACT Sum rules for products of two,three and four QCD currents are derived using chiral symmetry at infinite momentum in the large-N limit.These exact relations among meson decay constants,axialvector couplings and masses determine the asymptotic behavior of an infinite number of QCD correlators.The familiar spectral function sum rules for products of two QCD currents are among the relations derived.With this precise knowledge of asymptotic behavior,an infinite number of large-N QCD correlators can be constructed using dispersion relations.A detailed derivation is given ofthe exact large-N pion vector form factor and forward pion-pion scattering amplitudes.PACS:11.30.Rd;12.38.Aw;11.55.Jy;11.30.Er1.IntroductionIn the large-N limit all QCD mesonic correlation functions are sums of tree graphs[1,2]. The machinery of dispersion relations,which relates correlators at different energies,is transparent in this context since all dispersion relations are saturated by single-particle states.Of course dispersion relations require knowledge of asymptotic behavior and tra-ditionally this knowledge has been assumed or taken from phenomenologically inspired models like Regge pole theory.What is lacking is a fundamental principle based in QCD which determines the asymptotic behavior of correlation functions.The dominance of single-particle states at large-N has a related advantage:this limit is particularly suited to a symmetry point of view since it is easy to place single-particle states in representations of a symmetry group and extract consequences for matrix elements.In fact,there is a deep relation between asymptotic behavior and symmetries.For instance,requiring asymptotic behavior consistent with unitarity implies a spin-flavor algebra which constrains the spec-trum of the large-N baryons[3].Another classic example is string theory,where both the spectrum of states and the asymptotic behavior of S-matrix elements(which are QCD-like)are determined by symmetries on the string world sheet.In this paper I show that in large-N QCD,chiral symmetry determines the asymptotic behavior of an infinite number of mesonic correlation functions.Armed with this asymptotic information,exact large-N correlators can be constructed using dispersion relations in a facile way.There are well known QCD predictions of asymptotic behavior for products of two currents with nontrivial chiral quantum numbers[4].These asymptotic constraints—known as spectral function sum rules(SFSR’s)—are naturally derived in the context of the operator product expansion(OPE)[5].The fundamental ingredient of their derivation is one of the remarkable properties of the OPE:coefficient functions in the OPE are insensitive to vacuum properties and therefore transform with respect to the full unbroken global symmetry group of the underlying theory,regardless of whether this symmetry is spontaneously broken in the low-energy theory[5].Hence the SFSR’s exist because the full global chiral symmetry group of QCD can be used for classification purposes in the region of large Euclidean momenta.One might then wonder whether the SFSR’s can be derived directly using chiral symmetry without recourse to the OPE.In full QCD a direct symmetry approach is hindered by complicated analyticities arising from various continuum contributions to the two-point functions.On the other hand,in the large-N limit,the two-point functions are determined by the exchange of an infinite number of single-particle states.I will show that the SFSR’s at large-N are purely consequences of chiral symmetry and are easily derived directly from the chiral representation theory and/or the chiral algebra without recourse to the OPE.Relevant in this respect is a special classof Lorentz frames in which a component of momentum is taken to be very large compared to typical hadronic scales(e.g.the infinite momentum frame).The purpose of this work is not tofind an amusing way of reproducing the SFSR’s. Rather it is to see if the direct symmetry approach leads to other sum rules which have been overlooked or simply obscured by the complexities of the OPE.I will derive an infinite number of(new)sum rules in the large-N limit.These sum rules rely on precisely the same theoretical assumptions as the SFSR’s and in fact depend on the SFSR’s for consistency.Just as the SFSR’s determine the asymptotic behavior of two-point functions which in turn can be constructed using dispersion relations,all of the new chiral sum rules determine the asymptotic behavior of specific three-or four-point functions in large-N QCD.Dispersion relations can then be used to construct the exact correlators.I will consider several examples in detail.In section2I review some basic relevant facts about large-N QCD.The technology of the operator product expansion(OPE)is then used to review the derivation of the spectral function sum rules.In section3I derive the SFSR’s and an infinite number of new sum rules directly using chiral symmetry and its representation theory in the infinite momentum frame.In section4I identify some of these sum rules with asymptotic constraints and explicitly construct exact large-N correlators using dispersion relations.I conclude in section5.In an appendix I derive a set of sum rules directly from the chiral algebra.2.Spectral Function Sum Rules at Large-NThis section reviews known material which is essential for what follows.In this paper I consider QCD with two masslessflavors.This theory has an SU(2)L×SU(2)R chiral sym-metry.In the large-N limit,assuming confinement,this symmetry breaks spontaneously to the SU(2)V isospin subgroup[6].(Strictly speaking,large-N QCD has a U(2)L×U(2)R chiral symmetry.In this paper the additional Goldstone mode and its associated U(1) charge are ignored.)In the large-N limit the mesons have the most general quantum num-bers of the quark bilinear¯QΓQ whereΓis some arbitrary spin structure.This means that all mesons have zero or unit isospin.The corresponding chiral transformation properties of the mesons will be discussed below.The SU(2)×SU(2)algebra is expressed via the commutation relations[Q5a,Q5b]=iǫabc T c[T a,Q5b]=iǫabc Q5c[T a,T b]=iǫabc T c,(1)where T a are SU(2)V generators normalized so that T a=τa/2and Q5a are the SU(2)A generators that are broken by the vacuum,at rest.Note that V=L+R and A=L−R. The conserved QCD currents associated with these charges are:A aµ=¯Qγµγ5T a QV aµ=¯QγµT a Q.(2) These currents satisfy the commutation relations[Q a5,V bµ]=iǫabc A cµ[Q a5,A bµ]=iǫabc V cµ.(3) Therefore,V aµand A aµform a complete(six dimensional)chiral multiplet;they transform as(3,1)⊕(1,3)with respect to SU(2)×SU(2).In the low-energy theory,the axialvector current A aµhas a nonvanishing matrix element between the vacuum and the Goldstone boson(pion)states:0|A aµ|πb =δab Fπpµ.(4) Both conserved currents have nonvanishing matrix elements between the vacuum and vec-tor meson states:0|A aµ|A b (λ)=δab F A M Aǫ(λ)µ0|V aµ|V b (λ)=δab F V M Vǫ(λ)µ(5) whereǫ(λ)µis the vector meson polarization vector andλis a helicity label.Consider the time-ordered product of currents:ΠµνLR(q)δab=4i d4xe iqx 0|T Lµa(x)Rνb(0)|0 ,(6) where the SU(2)L,R QCD currents L and R transform as(3,1)and(1,3),respectively, with respect to SU(2)×SU(2).It follows thatΠLR transforms as(3,3).Lorentz invariance and current conservation allow the decompositionΠµνLR(q)=(qµqν−gµνq2)ΠLR(Q2),(7) where Q2=−q2.Because the functionΠLR(Q2)carries nontrivial chiral quantum numbers, it vanishes to all orders in QCD perturbation theory.Hence,asymptotic freedom implies ΠLR(Q2)→0as Q2→∞andΠLR satisfies the unsubtracted dispersion relation:ΠLR(Q2)=1t+Q2.(8)In the OPE we have the expansion,ΠLR(Q2)−−→Q2→∞∞n=1O d=2n(3,3)where the coefficients O d=2n(3,3),of mass dimension d=2n transform as(3,3)with respect to SU(2)×SU(2).This is the case regardless of whether SU(2)×SU(2)is spontaneously broken in the low-energy theory.The coefficient functions always transform with respect to the full unbroken symmetry group of the underlying theory[5].In the large-N limitΠLR is given by a sum of single-particle states[7].Using Eq.(4)-Eq.(7)one obtains−1Q2+M2A− V Q2F2VQ2+∞m=1 V F2V M2m V− A F2A M2m A (−1)m1Some interesting consequences of these sum rules when truncated to a small number of states are discussed in Ref.[7].Σ= Σ+ VA Σ = ΣV A Figure 1:Spectral function sum rules for products of two QCD currents.The top diagram is SFSR1and the bottom diagram is SFSR2.The crossrepresents a mass-squared insertion.ΠLR (q 2)=V F 2V M 4V q 4(M 2A −q 2),(15)which manifests the correct asymptotic behavior.Since narrow single-particle states are all that survive in the large-N limit,one may wonder whether it is possible to obtain the SFSR’s directly using SU (2)×SU (2)symmetry,with no reference to the OPE.One might expect that the pion together with the infinite number of vector and axialvector states in some sense fill out a representation of the full unbroken SU (2)×SU (2)group.We will see that this is the case.I reiterate that the main motivation of the symmetry approach is to find sum rules for other n -point functions in a simple way.For instance,consider adding an external pion line to one of the diagrams in Fig.1so that a vector current turns into an axial current (see Fig.4).Intuitively is would seem that there should be a sum rule analogous to the SFSR’s for such a graph.This is in fact the case,as we will see.3.Chiral Symmetry at Infinite Momentum 3.1The Infinite Momentum FrameThe fact that OPE coefficients do not feel vacuum properties enables the derivation of the SFSR’s by making use of the full unbroken chiral symmetry group [5].In order to obtain the SFSR’s using symmetry arguments it is necessary to work directly with the matrix elements of the currents between the states and the vacuum,rather than with the correlation function ΠLR .But clearly the manner in which the meson decay constants are defined through the matrix elements of Eq.(4)and Eq.(5)implies an asymmetry between states of different spin.The symmetric appearance of the decay constants in the SFSR’s (for instance,SFSR1is invariant with respect to F A ↔F πfor each A )is a consequence of taking Q 2→∞in the OPE.In studying matrix elements it is therefore convenient to work in a Lorentz frame in which the pions and the vector mesons appear the same.This is easy to achieve.At rest we define the vector meson polarization vectors:Figure2:Collinear kinematics of the infinite momentum frame wherep3→∞.Conservation of helicity,λ,follows from invariance with respect torotations about the3-axis.Parity takes A to B.A rotation byπabout the2-axis restores the original configuration with the sign of the spin(helicity)reversed.ǫ(+)µ=(0,1,0,0)ǫ(−)µ=(0,0,1,0)ǫ(0)µ=(0,0,0,1).(16) Consider boosting all particles along the3-axis(or the observer along the negative3-axis) to pµ=(p0,0,0,p3).We then haveǫ(±)µunchanged andǫ(0)µ=(p3M).(17)Now as we take p3→∞,ǫ(0)µ=pµp3)(18)and we have0|A aµ|πb =δab Fπpµ0|A aµ|A b (0)=δab F A pµ0|V aµ|V b (0)=δab F V pµ.(19)Therefore,as the momentum is taken large compared to the mass scales in the problem, theλ=0vector mesons act like Goldstone bosons.This is familiar from the Goldstone boson equivalence theorem in electroweak physics[8].The kinematical conditions described above are known as the infinite momentum frame.One advantage of this frame is that boost invariance is preserved along the3-direction.We can use boost invariance to define the amplitudes:0|A a0−A a3|α ≡(p0−p3) 0|A a|α0|V a0−V a3|α ≡(p0−p3) 0|V a|α ,(20) where the matrix elements of V a and A a are constants andαrepresents a physical meson state.It then follows from Eq.(19)and Eq.(20)that0|A a|πb =δab Fπ 0|A a|A b (0)=δab F A 0|V a|V b (0)=δab F V.(21)To reiterate,the matrix elements of the pions and theλ=0components of the vector mesons look the same in the infinite momentum frame.This is ideal when using symmetry arguments to relate matrix elements.3.2Selection Rules at Infinite MomentumWe have the following selection rules at infinite momentum2:•Invariance with respect to rotations about the3-axis implies helicity conservation.We can look at sectors of states labelled by helicity.This explains why theλ=0components of the vector mesons are on par with the pions at infinite momentum;allλ=0states look like scalars.We will be interested in the helicity conserving pion transition matrix elements between physical meson statesβandαwhich are related to the axial charges in the infinite momentum frame by:M a(p′λ′β,pλα)=(Fπ)−1(m2α−m2β)[(λ) β|Q5a|α (λ)]δλ′λ.(22)•The axial charges annihilate the vacuum in the broken phase:Q a5|0 =0.(23) This means that the chiral algebra is useful for classifying hadron states3.We will dis-cuss the manner in which symmetry breaking effects appear in this frame below.Since meson states are quark bilinears in the large-N limit and quarks transform as(2,1)and (1,2),mesons transform as combinations of(2,2),(3,1),(1,3)and(1,1)irreducible rep-resentations of SU(2)×SU(2).Charge conjugation leaves(2,2)and(1,1)unchangedand interchanges(1,3)and(3,1).Physical meson states have definite charge conjuga-tion,C,and isospin,I,and therefore are linear combinations of the isovectors|(2,2)a , {|(1,3)a −|(3,1)a }/√2≡|A a and the isoscalars |(2,2)0 and|(1,1) .Roman subscripts are isospin indices.Only|V a changes sign under charge conjugation.The action of the generators on the states of definite chirality can be obtained using tensor analysis:i A a|Q5b|V c j=iǫabcδij i (2,2)a|Q5b|(2,2)0 j=iδabδij(24a)i (2,2)a|T b|(2,2)c j=i A a|T b|A c j=i V a|T b|V c j=iǫabcδij.(24b)•Invariance with respect to a combined space inversion and rotation throughπabout an axis perpendicular to the3-axis implies:(λ) α|Q5a|β (λ)=(−λ) α|Q5a|β (−λ)PαPβ(−1)Jβ−Jα+1(25) where P and J are parity and spin,respectively.Forλ=0this implies the selection rule:ηβ=−ηα(26) else α|Q5a|β =0,whereη≡P(−1)J is normality.From Eq.(22)it then follows that only (λ=0)states of opposite normality communicate by pion exchange.•Conservation of G-parity implies the selection ruleGβ=−Gα(27) else α|Q5a|β =0.The product of normality and G-parity is a symmetry.Hence it follows that physical meson states fall into representations labelled by Gη.A geometric picture of several of the selection rules is given in Fig.2.3.3The Representation TheoryUsing the infinite momentum frame selection rules it is straightforward to identify the most general pion SU(2)×SU(2)representation in the large-N limit.Since the pion is a Lorentz scalar all states in its chiral representation haveλ=0.This representation can contain meson states of any spin.Of course all mesons haveλ=0components.Too,the pion has Gη=+1as do all states in its representation.Since the pion has I=1and C=+1it is a linear combination of|(2,2)a and|A a states.The other physical states which are linear combinations of|(2,2)a and|A a states necessarily haveη=−1andG Pηπ−−−P−−−A−+−V+−+ S+++|V a (0)1=m j=1w1j|V a j...|V a (0)m=m j=1w mj|V a j(28d)|S (0)1=n+m i=m+1v(m+1)i|(2,2)0 i+m+n+o i=m+n+1v(m+1)i|(1,1) i...|S (0)n+o=m+n i=m+1v(m+n+o)i|(2,2)0 i+m+n+o i=m+n+1v(m+n+o)i|(1,1) i(28e) whereˆu,ˆw andˆv are real mixing matrices,subject to the orthonormality relationsk=1u ik u jk=m k=1w ik w jk=m+n+o k=m+1v ik v jk=δij.(29)n+mThe dimension of the most general representation is Dim=2(2n+3m)+o where n is the number of(2,2)representations,m is the number of(1,3)⊕(3,1)representations,and o is the number of(1,1)representations.We work with Dim=finite but it should be understood that strictly speaking Dim=∞in large-N QCD.It is convenient to group the pion,the axialvector states and the pseudoscalar states into the n+m-component vectorΠi=(π,A2...A l,P l+1...P n+m).We can then express the pion multiplet in the compact form:|Πa i=m j=1u ij|A a j+n+m j=m+1u ij|(2,2)a j i=1...(m+n)(30a)|V a i=m j=1w ij|V a j i=1...m(30b)|S i=n+m j=m+1v ij|(2,2)0 j+m+n+o j=n+m+1v ij|(1,1) j i=1...(n+o).(30c)The helicity label has been supressed for ease of notation.It will be reintroduced below. The pion and vector meson decay constants generalized from Eq.(21)are given by0|A a|Πb i=δab FΠi(31a)0|V a|V b i=δab F V i(31b)and we now define the matrix elements of definite chirality:0|V a|V b j= 0|A a|A b j≡δab F j(32a)0|V a|A b j= 0|A a|V b j=0.(32b) Now using Eq.(30),Eq.(31)and Eq.(32)givesFΠi =mj=1u ij F j(33a)F Vi =mj=1w ij F j.(33b)Of course these decay constants are nonvanishing only for pion and vector meson states while the sums over states span all spins according to Table1.The axialvector couplings are defined byi Πb|Q5a|S j=−iδab G S jΠi/Fπ(34a)i Πb|Q5a|V c j=−iǫabc G V jΠi/Fπ.(34b) From Eq.(24),Eq.(30)and Eq.(34)it follows thatG SjΠi /Fπ=−m+nk=m+1u ik v jk(35a)G VjΠi /Fπ=mk=1u ik w jk.(35b)Note that the axial couplings are completely determined by the mixing matrices.It is now easy tofind relations that are independent of mixing matrices by contracting Eq.(33)and Eq.(35).3.4Chiral Symmetry ConstraintsThere is one relation involving decay constants only:m+ni=1FΠi FΠi=m i=1F V i F V i.(36) There are three relations which contract one decay constant with one axial coupling:mi=1F V i G V iΠk=FπFΠk(37a) m+nj=1FΠjG ViΠj=FπF Vi(37b)m+nj=1FΠj G S iΠj=0.(37c) There are three relations which contract two axial couplings:m+n+oj=m+1G S jΠi G S jΠk+m j=1G V jΠi G V jΠk=F2πδik(38a)m+ni=1G VjΠiG VkΠi=F2πδjk(38b)m+ni=1G S jΠi G V kΠi=0.(38c)These are the basic sum rules.(In an appendix these sum rules are derived directly from thechiral algebra.)They can be further contracted with axial couplings and decay constants to give other sum rules.In compact notation the complete set of sum rules which follow from Eq.(36)-Eq.(38)is:Two-point functions(SFSR1):V F V G Vπ=F2π(40a)V F V G(0)V A i=FπF A i(40b)V F V G(0)V P i=0(40c)A F A G(0)V i A+FπG V iπ=FπF V i(40d) A F A G(0)S i A+FπG S iπ=0(40e)ππΣV=πΣV=A iiπVΣ=iAπV i+iπΣ =iAπ+ 0iSπΣV=iPFigure3:Sum rules for three-point functions with one external current(Eq.(40)).The uppermost sum rule constrains the pion vector form factor.Three-point functions(two external currents):S G2Sπ+ V G2Vπ=F2π(42a)S G SπG(0)SA i+ V G VπG(0)V A i=0(42b)S G SπG(0)SP i+ V G VπG(0)V P i=0(42c)SG(0)SP iG(0)SA k+ V G(0)V P i G(0)V A k=0(42d) S G(0)SA i G(0)SA k+ V G(0)V A i G(0)V A k=F2πδik(42e)S G(0)SP i G(0)SP k+ V G(0)V P i G(0)V P k=F2πδik(42f)Σ =ΣV, AAFigure 4:Sum rule for three-point functions with two external currents(Eq.(41)).G V i πG V k π+A G(0)V i AG(0)V k A+ PG (0)V i P G (0)V k P =F 2πδik (42g )G S i πG V k π+AG (0)S i A G (0)V k A +P G (0)S i P G (0)V k P =0(42h )Four-point functions (one external current):A,S,A ′F AG (0)SAG (0)SA ′F A ′+A,V,A ′F A G(0)V AG (0)V A ′F A ′=F 2πAF 2A(44a )V,V ′F VG V πG V ′πF V ′+V,A,V ′F V G(0)V AG(0)V ′AF V ′+V,P,V ′F VG (0)V PG(0)V ′PF V ′=F 2πVF 2V .(44b )We have reintroduced the helicity label to denote those couplings which generally havenonzero helicity components which are not constrained by the pion representation.The first sum rule is the first spectral function sum rule,SFSR1.Here we see that there are an infinite number of additional sum rules which are at precisely the same level of theoretical rigor.The new sum rules are illustrated diagramatically in Fig.3(Eq.(40)),Fig.4(Eq.(41)),Fig.5(Eq.(42)),Fig.6(Eq.(43))and Fig.7(Eq.(44)).=S+1V=S+V=+=+kkδikδikV V kkS S kk=S+V 0ii=S+V 0i=S+VkkiP 0=S+VkkδikP P Figure 5:Sum rules for four-point functions with no external currents(Eq.(42)).Dotted lines are pions.3.5Symmetry Breaking ConstraintsGiven that the chiral charges annihilate the vacuum at infinite momentum,one may wonder how symmetry breaking manifests itself.Of course it must be present in order to split states within chiral multiplets.At infinite momentum spontaneous symmetry breaking implies[Q 5a ,ˆH ∞QCD ]=0(45)where ˆH ∞QCD is the QCD Hamiltonian at infinite momentum.This is simply the state-ment —familiar to practitioners of light-front field theory—that at infinite momentum,the effects of spontaneous symmetry breaking appear as explicit breaking terms in the Hamiltonian.At infinite momentum the Hamiltonian can be formally expanded as:ˆH ∞QCD=2|ˆP|+O (|ˆP|−3).(46)Therefore,spontaneous chiral symmetry breaking implies[Q 5a ,ˆM 2]=0(47)where ˆM2is the hadronic mass-squared matrix.The general solution to Eq.(47)is ˆM 2=ˆM 2(1,1)+RˆM 2R(48)ΣA, SΣA,V= +ΣA, SΣA,V=+A iV =iΣ VV i Σ+= Σ VΣ +S kS kiiΣA, SΣA,V=+iiFigure 6:Sum rules for four-point functions with one external current(Eq.(43)).Dotted lines are pions.where (1,1)is of course the singlet representation and R is a nontrivial SU (2)×SU (2)representation.There is no sense in which ˆM 2Ris small.From the allowed representations for the states it is straightforward to show that —without loss of generality—the most general symmetry breaking mass-squared matrix can be written asRˆM 2R =ˆM 2(2,2)+ˆM 2(3,3).(49)The matrix elements of ˆM 2(1,1)between the states of definite chirality are defined as i V a |ˆM 2(1,1)|V b j=i A a |ˆM 2(1,1)|A b j ≡δab m 2ij(50a )i (2,2)α|ˆM 2(1,1)|(2,2)β j ≡δαβn 2ij(50b )i (1,1)|ˆM 2(1,1)|(1,1) j≡o 2ij .(50c )The matrix elements of ˆM 2(2,2)are defined as i A a |ˆM 2(2,2)|(2,2)b j=i (2,2)a |ˆM 2(2,2)|A b j ≡δab ¯n 2ij(51a )ΣΣΣ =+Σ+Σ V=ΣA V, A, P, V’Figure 7:Sum rules for four-point functions with two external currents(Eq.(44)).Dotted lines are pions.i (1,1)|ˆM 2(2,2)|(2,2)0 j=i (2,2)0|ˆM 2(2,2)|(1,1) j ≡¯o 2ij .(51b )The matrix elements of ˆM 2(3,3)are defined as −i V a |ˆM 2(3,3)|V b j =i A a |ˆM 2(3,3)|A b j ≡δab ¯m 2ij .(52)The physical masses of the pion multiplet states arei Πa |ˆM2|Πb j =δij δab M 2Πi (53a )i V a |ˆM 2|V b j =δij δab M 2V i(53b )i S |ˆM 2|S j =δij M 2S i,(53c )from which we obtain using Eq.(30),Eq.(50),Eq.(51),and Eq.(52):δij M 2Πi=m k =1ml =1u ik u jl {m 2kl +¯m 2kl }+m +n k =m +1m +nl =m +1u ik u jl n 2kl +2m k =1m +nl =m +1u ik u jl ¯n 2kl(54a )δij M 2V i =m k =1ml =1w ik w jl {m 2kl −¯m 2kl }(54b )δij M 2S i =m +nk =m +1m +nl =m +1v ik v jl n 2kl+m +n +ok =n +m +1m +n +ol =n +m +1v ik v jl o 2kl+2m +nk =m +1m +n +ol =n +m +1v ik v jl ¯o 2kl .(54c )Contracting with the decay constants of Eq.(33)it is then easy to obtain the sum rule:m i =1m j =1F V i F V j δij M 2V i−m +n i =1m +nj =1F Πi F Πj δij M 2Πi=−2m k =1ml =1F k F l ¯m 2kl .(55)We also have the relationn+mi=1n+m j=1FΠi FΠj i Πa|ˆM2(3,3)|Πb j=−m i=1m j=1F V i F V j i V a|ˆM2(3,3)|V b j=δab m k=1m l=1F k F l¯m2kl.(56) In compact notation,combining Eq.(55)and Eq.(56),we haveV F2V M2V− A F2A M2A=−2ΠΠ′FΠFΠ′Π|ˆM2(3,3)|Π′ .(57)Comparing with the OPE in Eq.(13b)givesO d=4(3,3)=2 Π Π′FΠFΠ′ Π|ˆM2(3,3)|Π′ .(58)Of course this OPE coefficient vanishes in QCD.ThereforeˆM2(3,3)=0andˆM2=ˆM2(1,1)+ˆM2(2,2)(59)in large-N QCD.This chiral decomposition of the mass-squared matrix therefore implies SFSR2,V F2V M2V− A F2A M2A=0.(60) It is interesting that one can use the OPE to constrain the form of the meson mass-squared matrix.This result justifies the assumption of Coleman and Witten that the order parameter of chiral symmetry breaking transforms purely as(2,2)in the large-N limit[6].It is now straightforward tofind new relations that are independent of mixing matrices by contracting Eq.(35)with Eq.(54):m+n+oj=m+1G S jΠi G S jΠk(M2Πi+M2Πk−2M2S j)−m j=1G V jΠi G V jΠk(M2Πi+M2Πk−2M2V j)=0(61a)m+nk=1G V iΠk G V jΠk(M2V i+M2V j−2M2Πk)=0(61b)m+nk=1G S iΠk G V jΠk(M2V j−M2S i)=0.(61c)=SV=SV=SV=iikkkk-+--=SV 0-ii=SVkk-=SVkk-P Figure 8:Sum rules for four-point functions with no external currents(Eq.(64)).Dotted lines are pions.Crosses represent mass-squared inser-tions.The last sum rule is trivially satisfied via Eq.(38c).(In an appendix these sum rules are derived directly from the chiral algebra.)It is also straightforward to relate the axial couplings and masses to the symmetry breaking parameters.For instance,one findsi Πa |ˆM 2(2,2)|Πb jF 2π=δab [ml =1G V l Πi G V l Πj (M 2Πi +M 2Πj −2M 2V l)],(62)a result which we will return to below.In compact notation the complete set of mass-squared sum rules which follow from Eq.(57)and Eq.(61)are:Two-point functions (SFSR2):S G 2S πM 2S −VG 2V πM 2V =0(64a )SG S πG(0)SA i(M 2A i −2M 2S )−VG V πG (0)V A i (M 2A i −2M 2V )=0(64b )SG S πG (0)SP i (M 2P i −2M 2S )−VG V πG (0)V P i (M 2P i −2M 2V )=0(64c )ΣA, SΣA,V=Σ A, SΣA,V=A iA iV =iΣ VV iΣ+--Σ A, SΣA,V=-iiFigure 9:Sum rules for four-point functions with one external current(Eq.(65)).Dotted lines are pions.Crosses represent mass-squared inser-tions.S G(0)SA iG(0)SA k(M 2A i +M 2A k −2M 2S )−VG (0)V A i G (0)V A k (M 2A i +M 2A k −2M 2V )=0(64d )SG (0)SP i G (0)SP k (M 2P i +M 2P k −2M 2S )−V G (0)V P i G (0)V P k (M 2P i +M 2P k −2M 2V )=0(64e )SG (0)SP i G (0)SA k (M 2P i +M 2A k −2M 2S )−VG (0)V P i G (0)V A k (M 2P i +M 2A k −2M 2V )=0(64f )G V i πG V k π(M 2V i +M 2V k)+AG (0)V i A G (0)V k A (M 2V i +M 2V k −2M 2A )+ PG (0)V i P G (0)V k P (M 2V i +M 2V k−2M 2P )=0(64g )Four-point functions (one external current):。

Generalized Discrete Timed Automata:Decidable Approximations for Safety VerificationZhe Dang,Oscar H.Ibarra,and Richard A.KemmererSchool of Electrical Engineering and Computer ScienceWashington State University,Pullman,WA99164zdang@Department of Computer ScienceUniversity of California,Santa Barbara,CA93106ibarra,kemm@Abstract.We consider generalized discrete timed automata with general linearrelations over clocks and parameterized constants as clock constraints and withparameterized durations.We look at three approximation techniques(i.e.,the-reset-bounded approximation,the-bounded approximation,and the-crossing-bounded approximation),and derive automata-theoretic characteriza-tions of the binary reachability under these approximations.The characteriza-tions allow us to show that the safety analysis problem is decidable for general-ized discrete timed automata with unit durations and for deterministic generalizeddiscrete timed automata with parameterized durations.An example specificationwritten in ASTRAL is used to run a number of experiments using one of theapproximation techniques.1IntroductionAs a standard model for analyzing real-time systems,timed automata[3]have received enormous attention during the past decade.A timed automaton can be considered as a finite automaton augmented with afinite number of clocks.The clocks can be reset or progress at the same rate,and can be tested against clock constraints in the form of clock regions(i.e.,comparisons of a clock or the difference of two clocks against an integer constant,e.g.,,where and are clocks.).A fundamental result in the theory of timed automata shows that region reachability for timed automata is decidable[3]. This result has been very useful in defining various real-time logics,appropriate model checking algorithms,and tools[2,4,15,16,23,24,28]for verifying real-time systems (see[1]for a survey).However,not every real-time system can be modeled as a timed automaton.A complex(not necessarily large)real-time system might contain non-region clock con-straints,such as,where are clocks,and is a2parameterized(or unspecified)constant.Though there are practical needs to augment timed automata with more complex clock constraints,the“Turing computing”power of such augmented automata prevents automatic verification of properties such as reacha-bility[3,5].Therefore,decidable approximation techniques are of great interest,since these techniques would provide a user some degree of confidence in analyzing and de-bugging complex real-time specifications.In contrast to the most direct approximation techniques[5,10–12]that bound the number of transitions to afixed number,the ap-proximation techniques presented in this paper restrict the clock behaviors but do not necessarily bound the number of transition iterations to befinite.In this paper,we focus on timed automata with integer-valued clocks,i.e.,discrete timed automata,and extend them to generalized discrete timed automata by allowing general linear relations over clocks and parameterized constants as clock constraints. Furthermore,the duration of a transition can be a parameterized constant.These gen-eralizations have practical motivation.For example,many complex real-world speci-fications[6,8,10,21]written in the real-time specification language ASTRAL[6]use generalized clock constraints and parameterized durations in almost every specification. Therefore,the results presented in this paper may be useful in implementing an auto-matic specification debugger for complex real-time specification languages like AS-TRAL.We investigate three approximation techniques in this paper.The-reset-bounded approximation bounds the number of clock resets by a given positive integer for each clock.The-bounded approximation requires that,for each clock,the clock value is less than a given positive integer every time when the clock resets(but after the last reset,the clock can go unbounded.).Combining these two,the-crossing-bounded approximation requires that,for each clock,there are at most times that the clock re-sets after its value is greater or equal to.Given an approximation technique,we will focus on the binary reachability characterization of a generalized discrete timed au-tomaton.Binary reachability has recently been proposed and investigated for a number of timed systems[7,9,13,14],which makes it possible to reason about“non-region”properties for these systems.Wefirst show that,when a generalized discrete timed automaton has unit duration,the binary reachability under any one of the three approxi-mations is Presburger.Then by considering a generalized discrete timed automaton with parameterized durations,we show that the binary reachability under any one of the three approximations has a2DCM-padding when the machine is deterministic.Specifically, we show that the“padded language”for binary reachability can be accepted by a deter-ministic two-way counter machine with one reversal-bounded counter[19].The case for nondeterministic generalized discrete timed automata is open.These are good charac-terizations in the sense that the validity of Presburger formulas can be verified by these machines,and the emptiness problem for such machines is decidable.This allows us to establish,in principle,decidable results for the(Presburger)safety analysis problem for generalized discrete timed automata under the approximations.The2DCM-padding characterization is particularly interesting,since binary reachability is not necessarily semilinear.The paper is organized as follows.Section2gives the basic definitions and prelimi-nary results that are used in the paper.Section3presents the decidable characterization3 of binary reachability for generalized discrete timed automata with unit durations and with parameterized durations,under the approximations.Section4formulates the safety analysis problem and its decidability and proposes an open problem concerning nonde-terministic generalized discrete timed automata with parameterized durations.Section5 shows an example from an ASTRAL specification of the railroad crossing benchmark. Finally,in Section6,conclusions and future work are summarized.2Generalized Discrete Timed AutomataLet be afinite set of variables over the integers.An atomic linear relation on is defined as where and are integers.A linear relation on is con-structed from afinite number of atomic linear relations using and.denotes the set of all linear relations on.An atomic Presburger relation on is either an atomic linear relation on or a linear congruence mod,where and are integers.A Presburger formula can always be constructed from atomic Presburger relations using and.Let be the set of integers with denoting the set of nonnegative integers.A clock is a variable over.A generalized discrete timed automaton is a tuplewhere is afinite set of(control)states,is afinite set of clocks,is afinite set of parameterized constants,andis afinite set of edges.Each edge denotes a transition(or an edge)from to with enabling condition.is a parameterized constant indicating the duration of this transition.is the set clocks that are reset as a result of this transition.When is empty,the duration,which is a parameterized constant in or integer constant,must be positive.Clock resets take no time.Thus,when is not empty,the duration on this edge is simply0.The semantics is defined as follows.is called a configuration with being the state under this configuration.and denote the value of the clock and the value of the parameterized constant,respectively.Note that each clock and parameterized constant are nonnegative,i.e.,in.denotes a one-step transition along an edge in satisfying–Constant values do not change,i.e.,for each,,–The state is set to a new location,i.e.,–Each clock changes according to the edge given.When there is no clock reset on this edge,i.e.,,all the clocks synchronously progress by the amount of the duration,i.e.,for each,.In this case,the duration is positive,i.e.,.When,clocks in reset to0and all the other clocks areunchanged.Thus,clock resets take no time.That is,for each,and for each,.–The enabling condition is satisfied.That is,is true.We simply write if can reach by a one-step transition.is deterministic if for any configuration there is at most one such that.A pathsatisfies for each.Write if reaches through a path.In the definition of,there is no input.This is because the input is always one-way for4timed automata and,therefore,input symbols can be built into the control states.has parameterized durations.If we restrict each duration on an edge without clock resets to be1,then is called a generalized discrete timed automaton with unit durations.x-y+z>B{}, Cs1s22x-y<A{z}Fig.1.An example generalized discrete timed automatonExample1.Figure1shows a generalized discrete timed automaton with parameterized constants and clocks.The automaton has two transitions(that is from state to state)and(that is from state to state).As shown in the figure,has an enabling condition of,an empty clock reset set,and a parameterized duration.has an enabling condition of,a clock reset set,and a zero duration.According to the semantics,the following is an instance of an execution of the automaton when:5 being considered.The following are some known results about the state reachability problem.Theorem1.(1).The state reachability problem is decidable for standard timed au-tomata[3].(2).The state reachability problem is decidable for parameterized timed automata with only one parameterized clock(but can have many unparameterized clocks)[5].(3).The state reachability problem is undecidable for generalized discrete timed automata.Actually,Theorem1(3)follows from the following special cases.Theorem2.(1).The state reachability problem is undecidable for standard timed au-tomata when we allow“”operations in clock constraints,e.g.,[3].(2).The state reachability problem is undecidable for parameterized timed au-tomata(with more than2parameterized clocks)[5].On the other hand,binary reachability is the set of all configuration pairs such that(we use in the paper).Characterizations of the binary reachability of timed automata have recently been established.Theorem3.(1).The binary reachability of timed automata with real valued clocks is definable in the additive theory over reals and integers[7,9].(2).The binary reachability of timed automata with integer valued clocks is defin-able by a Presburger formula[7,14].Characterizations of binary reachability help us to reason about non-region proper-ties of timed systems.For instance,consider the following property:“for every config-uration there exists a configuration such that and the clock in is the sum of clocks and in.”Though constraint is not in the form of a clock region,this property can be automatically verified for timed automata[7,9, 14].However,for generalized discrete timed automata,the binary reachability is too strong to have an interesting characterization.In particular,even the membership problem for binary reachability,i.e.,deciding whether two given configurations and satisfy for a generalized discrete timed automaton,is undecidable.This follows from the fact that a two-counter machine can be simulated by a generalized discrete timed automaton,as shown in the the proofs of Theorem2(1)(2)[3,5]. Thus,the membership problem can be reduced to the halting problem for two-counter machines,which is undecidable.Theorem4.The membership problem for binary reachability is undecidable for gen-eralized discrete timed automata.This undecidability result for generalized discrete timed automata leads us to con-sider the following three approximations of.Let and be givenfixed positive integers.Thefirst approximation is-reset-bounded reachability.A pathis called-reset-bounded if each clock resets at most times.Write if6reaches through an -reset-bounded path.The second approximation is -bounded reachability.A path is called -bounded if for each ,each ,Write if reaches through a -bounded path.The third approximation is-crossing-bounded reachability.A path is called -crossing-bounded if there are at most many ’s such thatWrite if reaches through a -crossing-bounded path.Figure 2shows the behaviors of a clock under the three approximations with some bound and .Figure 2(a)is the case for -reset-bounded approximation,where the clock is reset for at most four times.Figure 2(b)is the case for -bounded approx-imation,where the clock value always stays below the bound before the last reset (but there could be many resets).Figure 2(c)is the case for -crossing-bounded approximation,where the clock crosses the bound but for at most four times before the lastreset.B(a)(b)(c)Fig.2.Behaviors of a clock under the three approximations The main results in this paper show that the three approximations of binary reach-abilityhave decidable characterizations,i.e.,they can be accepted by a class of machines with a decidable emptiness problem.Before we proceed to show the results,some further definitions are needed.A nondeterministic multicounter machine (NCM)is a nondeterministic machine with a finite set of (control)states ,and a finite number of coun-ters with integer counter values.Each counter can add 1,subtract 1,or stay unchanged.These counter assignments are called standard assignments .can also test whether a counter is equal to,greater than,or less than an integer constant.These7 tests are called standard tests.When an NCM is used as a language recognizer,we attach a separate one-way read-only input tape to the machine and assign a state in as thefinal state.accepts an input iff it can reach thefinal state.It is well-known that counter machines with two counters have undecidable halting problem.Thus,we have to restrict the behaviors of the counters.One such restriction is to limit the num-ber of reversals a counter can move.A counter is-reversal-bounded if it changes mode between nondecreasing and nonincreasing at most times.A counter is-strong-reversal-bounded if it changes mode between strictly increasing,strictly decreasing and unchanged at most times.For instance,the following sequence of counter values:exhibits only one counter reversal.However, it has11strong reversals.is(strong-)reversal-bounded if each counter in is-(strong-)reversal-bounded for some.We note that a(strong-)reversal-bounded does not necessarily limit the number of moves or the number of reachable configurations to befinite.Let denote the configuration of an NCM(without input tape) when it is in state,counter has value for.Similar to ,the binary reachability of is defined as all the pairs of configura-tions of such that can reach.It is known that the emptiness problem for reversal-bounded NCMs(when used as language recognizers)is decidable[18].When a reversal-bounded NCM uses linear relations as tests instead of using standard tests,we have the following result. Theorem5.Suppose is a nondeterministic strong-reversal-bounded multicounter machine without input tape that uses linear relations(on the counters and parameter-ized constants)as tests instead of using standard tests.Then its binary reachabilityis Presburger[20].The machines defined above,when used as language recognizers,have a one-way input tape.Suppose a two-way input is used instead.Let2DCM()denote the class of deterministic machines with a two-way input tape and-reversal-bounded counters. Then the emptiness problem for2DCM()when and is undecidable[18]. An interesting special case is when,i.e.,there is only one counter.A language is 2DCM-recognizable if it can be accepted by a2DCM().Theorem6.The emptiness problem for2DCM-recognizable languages is decidable [19].It is still open whether Theorem6holds for nondeterministic machines.That is, whether the emptiness problem for2NCM(1,r),which is a nondeterministic-reversal-bounded one counter machine with a two-way input tape,is decidable.Given a generalized discrete timed automaton,consider a subset of configuration pairs,.One can look at as some sort of reachability relation.For a given,if a2DCM()can be effectively constructed such that for every,is in iff there exists a such that accepts,then we say that has a2DCM-padding.From Theorem6,it is routine to show the following lemma.8Lemma1.(1).The emptiness problem for having2DCM-paddings is decid-able.(2).If and have2DCM-paddings,then so do the join and the composition.3Decidability ResultsDenote to be the binary reachability of a generalized discrete timed automaton through a path without clock resets.Note that itself can be considered as a nonde-terministic multicounter machine.However,tests in,which are linear relations on clocks and parameterized constants,are not standard.Assignments in in the form of with a parameterized constant are also not standard.But,if has only unit durations,the assignments are standard except when a clock reset occurs,i.e.,.Lemma2.Suppose is a generalized discrete timed automaton with unit durations. Then is Presburger.Proof.Let be a generalized discrete timed automaton with unit durations.Note that when a path of contains no clock resets,i.e.,all the clocks increase at the same rate 1,clocks(understood as counters)are0-strong-reversal-bounded.The theorem follows from Theorem5.9(3).-crossing-bounded reachability:Consider a-crossing-bounded path of.From the definition,each clock in the path can only cross the bound-ary at most times.Thus,the path can be divided into at most phases with thefirst phases being-bounded paths.The last one is either a-bounded path or a path with no clock resets.Phases are connected by a one-step transition with at least one clock reset.Thus,the result follows from(1)and(2)above.10resulting automaton is written as.Notice that is deterministic.That is,each node in cannot have more than one outgoing edge,from the definition of and the fact that is deterministic.Given two configurations and at control states and ,respectively,a key observation is:iff in and both and have labelThat is,if and only if reaches in with thefirst configuration and the last configuration having the same label.The convexity of each label ensures that all the intermediate configurations also have the same label.Sinceis deterministic,any path from to can be written as for some,with the length of,and bounded by the number of control states.That is,the path can be concatenated by a short path and a short cycle followed by another short path.In order to show has a2DCM-padding,we construct the required2DCM(1,).operates on an input tape with the following format:,where and are string encodings of the two configurations and.,and are short sequences of control states,which could be empty,with started by and ended by .From the above key observation,only needs to check:(1).the input tape is in the correct format.That is,actually has a path. This can be done by using afinite table look-up.Configurations and agree on each(positive)parameterized constant.This can be done with exactly counter reversals.(2).and have the same label.Note that is a linear relation.Checking satisfiability of a linear relation over configurations and can be done by accumu-lating the counter while reading the value of each clock and constant stored inon the two-way input tape.Since isfixed,this needs only counter reversals,where is the number of atomic linear relations in.(3).also needs to check that for each,the net clock changeis the same.This needs counter reversals.Also,makes sure that the net clock change agrees with the net clock change on a path for some.Denote,, to be the sum of the durations of the edges in sequences,respectively.Let .The counter infirst calculates.This can be done by at most1counter reversal.Then tests whether the counter,and doing the same test again by subtracting from the counter until the test succeeds.If the counter is equal to0,then accepts,else rejects.Note that has only one counter,thus cannot be stored.Since is a short path,instead subtracts each duration(that is a constant,stored in)one by one in.Thus,accepts some input as above iff.This completes the construc-tion of the required2DCM(1,)machine.4Verification of Safety PropertiesConsider and,two sets of configurations of a generalized discrete timed automaton definable by Presburger formulas.If,starting from a configuration in,can only reach configurations in,then is a safety property with respect to the initial condition .The following is an example:“starting from a configuration satisfying,cannot reach a configuration satisfying.”The safety analysis problem is to determine whether is a safety property with respect to the initial condition.Theorem9.The safety analysis problem is decidable for generalized discrete timed automata with unit durations for any one of the following approximations:-reset-boundedness,-boundedness and-crossing-boundedness.Proof.Let be a generalized discrete timed automata with unit durations.Without loss of generality,consider-reset-bounded approximation.The safety analysis problem is equivalent to the existence of and such that and From Theorem7,is ing the fact that and are Presburger and Presburger formulas are closed under quantification,the theorem follows.Remark:Theorem10can be strengthened.The class of languages accepted by deter-ministic two-way counter machines with one reversal-bounded counter is closed under Boolean operations[19].It follows that Theorem11remains valid even if the sets of configurations(property)and(initial condition)are sets accepted by these ma-chines.It is desirable to consider the decidability of the safety analysis problem for gener-alized discrete timed automata under some special form but without using the approxi-mations.One such form is parameterized timed automata with only one parameterized clock.As stated in Theorem1(2),the state reachability problem is decidable.However, surprisingly,the safety analysis problem,i.e.,“Deciding whether is a safety property with respect to the initial condition for a parameterized timed automaton with only one parameterized clock,where both and are definable by Presburger formulas”is still open.This problem is closely related to the open problem of the decidability of the emptiness problem of2NCM(1,).This is also an example showing that binary reachability is much stronger than state reachability.Another form that can be consid-ered is generalized discrete timed automata with clock constraints(is an integer). Thus,they are standard timed automata with parameterized durations.The state reach-ability problem for these automata is decidable.The idea is to translate a transition with a parameterized duration into a loop with unit duration and add a new clock,which is tested against,and use Theorem1(2)since is the only parameterized clock.But again,it is open whether the safety property analysis problem is decidable for these automata.This problem is simpler than the previous open problem,but the answer is currently not clear.It is also worthwhile to point out that the characterizations of the binary reachabil-ity of deterministic generalized discrete timed automata under the approximations,as shown in Theorem8,are not necessarily Presburger.In fact,there are nonsemilinear languages that are2DCM(1,)-recognizable[19].5A Verification ExampleIn practice,allowing generalized clock constraints and parameterized durations makes it possible to specify more complex real-time systems.In this section,we take an ex-ample specification[22]of the railroad crossing benchmark[17],which is written in ASTRAL[6].The specification specifies a system consisting of a set of railroad tracks that intersect a street where cars may cross the tracks.A gate is located at the crossing to prevent cars from crossing the tracks when a train is near.A sensor on each track detects the arrival of trains on that track.The critical requirement of the system is that whenever a train is in the crossing the gate must be down,and when no train has been in between the sensors and the crossing for a reasonable amount of time,the gate must be up.The complete ASTRAL specification of the railroad crossing system,which was written by Paul Kolano,can be found at /˜dang.The ASTRAL specification includes a global specification and two process specifi-cations:one is Gate and the other is Sensor.A transition system is specified inside each process specification along with local assumptions and local properties.ASTRAL adopts a modularized view of the specified system at the level of correctness proofs:the global properties in the global specification can be verified by using global assumptions and local properties(instead of the actual transition behaviors)of each process instance; local properties of a process instance can be verified by using the local assumptions,the local imported variable clause and the transition system of the instance,without look-ing at the transition behaviors of the other process instances(A reader need not worry about the possibility of circular proofs.The ASTRAL proof theory is sound,see[8]for details.).We take the Sensor process specification to see how parameterized durations and parameterized clock constraints are used in ASTRAL.Sensor reports whether a trainis beyond the crossing,which is indicated by a Boolean variable train R.The process specification has only two transitions enter R,which have pa-rameterized durations enter dur,respectively.Transition enter inR sets train R back to FALSE after the slowest train moves out of the crossing.That is,TRANSITION exit_IENTRY[TIME:exit_dur]train_in_R&now-Start(enter_R)>=RIImin-exit_dur EXITtrain_in_R=FALSE,where now indicates the current time,Start(enterR,and RIImin is a parameterized constant indicating the time for the slowest train to move out of the crossing.One of the two transitions is nonde-terministically chosen tofire as long as the ENTRY condition of the chosen transition is satisfied.However,if neither of them isfireable,the process idles:now progresses by one time unit and train R does not change.But,according to this semantics,transition enterR completes.This is not the intended behavior of Sensor.In fact,it is desirable that enter R specified as above only tells what happened(set train R to TRUE)when a train comes,instead of the time when a train comes.The pattern of a train’s arrival is con-trolled by the environment of the process.In this specification,transition enterR must be separated by at least RIImin many time units.Assuming the environment holds,the safety property of the process is specified as a schedule:(now>=response_time&Call(enter_R,now-response_time) ->train_in_R)&(now>=RIImin&Call(enter_R,now-RIImin)->˜train_in_R).Thefirst conjunct of the schedule says that a train will be sensed within enterR is placed(note that,according to the assumptions on the parameterized constants,responsedur,and satisfying RIImin>=response dur.).The second conjunct of the schedule says that the sensor will be reset when the slowest train is beyond the crossing.It is assumed that initially now is0and train R is FALSE.We manually translated Sensor into a generalized discrete timed automaton and computed the transitive closure of the one-step transition of the automaton using the Omega Library[25],which is a tool to manipulate Presburger formulas.Experiments were run on a Sun workstation with4CPUs and256M real memory and512M swap memory.Unfortunately,the closure,even when the durations(enterdur)were set to1,could not be computed.Then,we used the-bounded approxi-mation on the automaton with.This time,the binary reachability can be computed in about one minute of CPU time and using170M memory.But the other two approximation approaches were still too expensive to calculate.6Conclusions and Future WorkWe studied generalized discrete timed automata with general linear relations over clocks and parameterized constants as clock constraints and with parameterized durations.We focused on three approximation techniques and automata-theoretic characterizations of binary reachability under these approximations.The characterizations allow us to show that the safety analysis problem is decidable with respect to generalized discrete timed automata with unit durations,and deterministic generalized discrete timed automata with parameterized durations(modulo the approximations).We used an example spec-ification written in ASTRAL to run a number of experiments using one of the approx-imation techniques.The results of the experiments show that further improvements to the approximations have to be developed,since currently they are not practical for large specifications.For future work,we want to investigate how the approximation techniques proposed in this paper can be combined with existing image-approximation techniques[12]in debugging infinite state systems.Solutions to this problem would lead to an implemen-tation of an effective specification debugger for large real-time specifications.Another research issue is how to extend the results in this paper to the case of generalized timed automata with dense clocks.Recent ideas in[9]may provide some insights. References1.R.Alur,“Timed automata”,CAV’99,LNCS1633,pp.8-222.R.Alur,C.Courcoubetis and D.Dill,“Model-checking in dense real time,”Information andComputation,104(1993)2-343.R.Alur and D.Dill,“A theory of timed automata,”TCS,126(1994)183-2364.R.Alur,T.A.Henzinger,“A really temporal logic,”J.ACM,41(1994)181-2045.R.Alur,T.A.Henzinger and M.Y.Vardi,“Parametric real-time reasoning,”STOC’93,pp.592-601。

SIAM J.I MAGING S CIENCES c 2008Society for Industrial and Applied Mathematics Vol.1,No.1,pp.143–168Bregman Iterative Algorithms for 1-Minimization with Applications toCompressed Sensing∗Wotao Yin†,Stanley Osher‡,Donald Goldfarb§,and Jerome Darbon‡Abstract.We propose simple and extremely efficient methods for solving the basis pursuit problem min{ u 1: Au=f,u∈R n},which is used in compressed sensing.Our methods are based on Bregman iterative regularization,and they give a very accurate solution after solving only a very small number ofAu−f k 22for given matrix A and vector instances of the unconstrained problem min u∈R nμ u 1+12f k.We show analytically that this iterative approach yields exact solutions in afinite number ofsteps and present numerical results that demonstrate that as few as two to six iterations are sufficient in most cases.Our approach is especially useful for many compressed sensing applications where matrix-vector operations involving A and A can be computed by fast transforms.Utilizing a fastfixed-point continuation solver that is based solely on such operations for solving the above unconstrained subproblem,we were able to quickly solve huge instances of compressed sensing problems on a standard PC.Key words. 1-minimization,basis pursuit,compressed sensing,iterative regularization,Bregman distances AMS subject classifications.49,90,65DOI.10.1137/0707039831.Introduction.Let A∈R m×n,f∈R m,and u∈R n.The basis pursuit problem[23] solves the constrained minimization problem{ u 1:Au=f}(1.1)(Basis Pursuit)minuto determine an 1-minimal solution u opt of the linear system Au=f,typically underdeter-mined;i.e.,m<n(in many cases,m n),and Au=f has more than one solution.The basis pursuit problem(1.1)arises in the applications of compressed sensing(CS).A recent burst of research in CS was led by Cand`e s and Romberg[12],Cand`e s,Romberg,and Tao[14],Cand`e s and Tao[16],Donoho[35],Donoho and Tanner[36],Tsaig and Donoho [88],and others[80,86].The fundamental principle of CS is that,through optimization,the sparsity of a signal can be exploited for recovering that signal from incomplete measurements of it.Let the vector¯u∈R n denote a highly sparse signal(i.e.,k= ¯u 0:=|{i:¯u i=0}| n).∗Received by the editors September28,2007;accepted for publication(in revised form)January29,2008; published electronically March20,2008./journals/siims/1-1/70398.html†Department of Computational and Applied Mathematics,Rice University,Houston,TX77005(wotao.yin@ ).This author’s research was supported by an internal faculty research grant from the Dean of Engineering at Rice University.‡Department of Mathematics,UCLA,Los Angeles,CA90095(sjo@,jerome@).The research of these authors was supported by ONR grant N000140710810.§Department of Industrial Engineering and Operations Research,Columbia University,New York,NY10027 (goldfarb@).This author’s research was supported in part by NSF grant DMS06-06712,ONR grant N00014-03-0514,and DOE grant GE-FG01-92ER-25126.143144W.YIN,S.OSHER,D.GOLDFARB,AND J.DARBON This principle states that one can encode¯u by a linear transform f=A¯u∈R m for some m greater than k but much smaller than n,and then recover¯u from f by solving(1.1).It is proved that the recovery is perfect;i.e.,the solution u opt=¯u for any¯u whenever k,m,n,and A satisfy certain conditions(e.g.,see[13,32,39,44,80,100,101]).While these conditions are computationally intractable to check,it was found in[15,16]and other work that the types of matrices A allowing a high compression ratio(i.e.,m n)include random matrices with independent and identically distributed(i.i.d.)entries and random ensembles of orthonormal transforms(e.g.,matrices formed from random sets of rows of the matrices corresponding to Fourier and cosine transforms).Recent applications of 1-minimization can be found in[51,84,91,92]for compressive imaging,[61,68,70,69,96]for MRI and CT,[3,4,50,54,76,93]for multisensor networks and distributive sensing,[63,65,66,77,87]for analog-to-information conversion,[83]for biosensing,[97]for microarray processing,and[24,25,98,99]for image decomposition and computer vision tasks. 1-minimization also has applications in image inpainting and missing data recovery;see[42,82,101],for example.Also nonconvex quasi– p-norm approaches for 0≤p<1have been proposed by Chartrand[20,21]and Chartrand and Yin[22].Problem(1.1)can be transformed into a linear program and then solved by conventional linear programming solvers.However,such solvers are not tailored for the matrices A that are large-scale and completely dense,or are formed by rows taken from orthonormal matrices corresponding to fast transforms so that Ax and A x can be computed by fast transforms. This,together with the fact that f may contain noise in certain applications,makes solving the unconstrained problem(1.2)minu μ u 1+12Au−f 22more preferable than solving the constrained problem(1.1)(e.g.,see[26,27,29,37,41,49, 58,62,89]).Hereafter,we use · ≡ · 2to denote the2-norm.In section2.1,we give a review of recent numerical methods for solving(1.2).Because(1.2)also allows the constraint Au=f to be relaxed,it is used when the measurement f is contaminated by encoding errors such as noise.However,when there is no encoding error,one must assign a tiny value to μto heavily weigh thefidelity term Au−f 2in order for Au=f to be nearly satisfied. Furthermore,one can show that the solution of(1.2)never equals that of(1.1)unless they both have the trivial solution0.In this paper,we introduce a simple method based on Bregman iterative regularization[73],which we review in section2.2,forfinding a solution of problem (1.1)by solving only a small number of instances of the unconstrained problem(1.2).Our numerical algorithm,based on this iterative method,calls the fastfixed-point continuation solver FPC[55,56]of(1.2),which involves only matrix-vector multiplications(or fast linear transforms)and componentwise shrinkages(defined in(2.4)).Using a moderate value for the penalty parameterμ,we were able to obtain a very accurate solution to the original basis pursuit problem(1.1)for a very small multiple of the cost of solving a single instance of(1.2).Our results can also be generalized to the constrained problem(1.3)minu{J(u):Au=f}BREGMAN ITERATIONS FOR 1-MINIMIZATION145 for other types of convex functions J(refer to section5).Specifically,a solution of(1.3)can be obtained through afinite number of the Bregman iterations of(1.4)minu μJ(u)+12Au−f 2.In addition,in section5.3,we also introduce a two-line algorithm(given in(5.19)and(5.20))also involving only matrix-vector multiplication and shrinkage operators that generates asequence{u k}that converges rapidly to an approximate solution of the basis pursuit problem (1.1).In fact,the numerical experiments in[34]indicate that this algorithm converges to atrue solution if the parameterμis large enough.Finally,preliminary experiments indicatethat our algorithms are robust with respect to a certain amount of noise.This is also impliedby our theoretical results stated in Theorems2.1and5.5.The rest of the paper is organized as follows.In section2,we summarize the existingmethods for solving the unconstrained problem(1.2)and provide some background on ourBregman iterative regularization scheme.The main Bregman iterative algorithm is describedin section3.1;its relationship to some previous work[95]is presented in section3.2;and itsconvergence is analyzed in section3.3.Numerical results are presented in section4.Finally,we extend our results to more general classes of problems in section5,including a descriptionand analysis of our linearized Bregman iterative scheme,and conclude the paper in section6.2.Background.2.1.Solving the unconstrained problem(1.2).Several recent algorithms can efficientlysolve(1.2)with large-scale data.The authors of GPSR[48],Figueiredo,Nowak,and Wright[49],reformulate(1.2)as a box-constrained quadratic program,to which they apply the gra-dient projection method with Barzilai–Borwein steps.The algorithm 1 s[64]by Kim et al.[62]was developed for an 1-regularization problem equivalent to(1.2).The authors apply aninterior-point method to a log-barrier formulation of(1.2).The main step in each interior-point iteration,which involves solving a system of linear equations,is accelerated by using apreconditioned conjugate gradient method,for which the authors developed an efficient pre-conditioner.In the code SPGL1[90],Van den Berg and Friedlander apply an iterative methodfor solving the LASSO problem[85],which minimizes Au−f subject to u 1≤σ,by using an increasing sequence ofσ-values in their algorithm to accelerate the computation.In[71], Nesterov proposes an accelerated multistep gradient method with an error convergence rate O(1/k2).Under some conditions,the greedy approach StOMP[37]by Donoho,Tsaig,Drori, and Starck can also quickly solve(1.2).A method widely used by many researchers to solve(1.2)or general 1-minimizationproblems of the form(2.1)minuμ u 1+H(u)for convex and differentiable functions H(·)is an iterative procedure based on shrinkage(alsocalled soft thresholding;see(2.4)below).This type of method was independently proposedand analyzed by Figueiredo and Nowak in[46,72]under the expectation-minimization frame-work for wavelet-based deconvolution,De Mol and Defrise[31]for wavelet inversion,Bect et al.146W.YIN,S.OSHER,D.GOLDFARB,AND J.DARBON in[5]using an auxiliary variable and the idea from Chambolle’s projection method[17],Elad in [38]and Elad et al.[40]for sparse representation and other related problems,Daubechies,De-frise,and De Mol in[29]through an optimization transfer technique,Combettes and Pesquet [26]using operator-splitting,Hale,Yin,and Zhang[55]also using operator-splitting com-bined with a continuation technique in their code FPC[56],Darbon and Osher[27]through an implicit PDE approach,and others.In addition,related applications and algorithms can be found in Adeyemi and Davies[1]for image sparse representation,Bioucas-Dias[6]for wavelet-based image deconvolution using a Gaussian scale mixture model,Bioucas-Dias and Figueiredo for a recent“two-step”shrinkage-based algorithm[7],Blumensath and Davies[8] for solving a cardinality constrained least-squares problem,Chambolle et al.[19]for image denoising,Daubechies,Fornasier,and Loris[30]for a direct and accelerated projected gradi-ent method,Elad,Matalon,and Zibulevsky in[41]for image denoising,Fadili and Starck[43] for sparse representation-based image deconvolution,Figueiredo and Nowak[47]for image deconvolution,Figueiredo,Bioucas-Dias,and Nowak[45]for wavelet-based image denoising using majorization-minimization algorithms,and Reeves and Kingsbury[78]for image coding.While all of these authors used different approaches,they all developed or used algorithms based on the iterative scheme(2.2)u k+1←arg minu μ u 1+12δku−(u k−δk∇H(u k))2for k=0,1,...,starting from a point u0.The parameterδk is positive and serves as the step size at iteration k.Since the unknown variable u is componentwise separable in problem (2.2),each of its components u i can be independently obtained by the shrinkage operation, which is also referred to as soft thresholding:(2.3)u k+1i=shrink((u k−δk∇H(u k))i,μδk),i=1,...,n,where for y,α∈R,we define(2.4)shrink(y,α):=sgn(y)max{|y|−α,0}=⎧⎪⎨⎪⎩y−α,y∈(α,∞), 0,y∈[−α,α],y+α,y∈(−∞,−α).Among the several approaches that can be used to derive(2.2),one of the simplest is the following:first,H(u)is approximated by itsfirst-order Taylor expansion at u k,H(u k)+ ∇H(u k),u−u k .Then,since this approximation is only accurate for u close to u k,an 2-penalty term u−u k 2/(2δk)is added to the objective;the resulting step is(2.5)u k+1←arg minu μ u 1+H(u k)+ ∇H(u k),u−u k +12δku−u k 2,which is equivalent to(2.2)because their objectives differ by only a constant.It is easy to see that the larger theδk,the larger the allowable distance between u k+1and u k.It was proved in[55]that{u k}given by(2.2)converges to an optimum of(1.4)at a q-linear1rate 1q stands for“quotient”;{x k}converges to x∗q-linearly if lim k x k+1−x∗ / x k−x∗ exists and is less than1.BREGMAN ITERATIONS FOR 1-MINIMIZATION147 under certain conditions on H andδk.Under weaker conditions,they also established r-linear convergence of{u k}based on previous work by Pang[74]and Luo and Tseng[67]on gradient projection methods.Furthermore,it was also proved in[55]that under mild conditions,the support and signs of u k convergefinitely;that is,there exists afinite number K such that {i:u k=0}={i:u opt=0}and sgn(u k)=sgn(u opt)for all k>K,where u opt denotes the solution of(1.4).However,an estimate or bound for K is not known.To improve the efficiency of the iterations(2.2),various techniques have been applied to (2.2),including generalizing(2.3)by using more parameters[41],performing line searches[49], and using a decreasing sequence ofμ-values[55].The last technique is called path following or continuation.While our algorithm does not depend on using a specific code,we chose to use FPC[56],one of the fastest codes,to solve each subproblem in(2.2).In[27],[94],and other work,the iterative procedure(2.2)is adapted for solving the total variation regularization problem(2.6)minuμT V(u)+H(u),where T V(u)denotes the total variation of u(see[102]for a definition of T V(u)and its properties).Specifically,the regularization termμ u 1in(2.2)is replaced byμT V(u),yielding(2.7)u k+1←arg minu μT V(u)+12δku−(u k−δk∇H(u k))2.Each subproblem(2.7)can be efficiently solved,for example,by one of the recent graph/network-based algorithms[18,28,53].In[27]Darbon and Osher also studied an algorithm obtainedby replacingμT V(u)in(2.7)by its Bregman distance(see section2.2)and proved that ifH(u)=0.5 Au−f 2,then{u k}converges to the solution of min u{T V(u):Au=f}.Their algorithm and results are described in section5.3.In the next subsection,we give an intro-duction to Bregman iterative regularization.2.2.Bregman iterative regularization.Bregman iterative regularization was introducedby Osher et al.[73]in the context of image processing;it was then extended to wavelet-baseddenoising[95],nonlinear inverse scale space in[10,11],and compressed sensing in MR imaging[59].The authors of[73]extend the Rudin–Osher–Fatemi[81]model(2.8)minu μ|∇u|+12u−b 2,where u is an unknown image,b is typically an input noisy measurement of a clean image ¯u,andμis a tuning parameter,into an iterative regularization model by using the Bregman distance(2.10)below based on the total variation functional:(2.9)J(u)=μT V(u)=μ|∇u|.Specifically,the Bregman distance[9]based on a convex functional J(·)between points u and v is defined as(2.10)D pJ(u,v)=J(u)−J(v)− p,u−v ,148W.YIN,S.OSHER,D.GOLDFARB,AND J.DARBON where p∈∂J(v)is some subgradient in the subdifferential of J at the point v.BecauseD pJ (u,v)=D pJ(v,u)in general,D pJ(u,v)is not a distance in the usual sense.However,itmeasures the closeness between u and v in the sense that D pJ (u,v)≥0and D pJ(u,v)≥D pJ (w,v)for all points w on the line segment connecting u and v.Instead of solving(2.8)once,the Bregman iterative regularization procedure of Osher et.al.[73]solves a sequence of convex problems(2.11)u k+1←minu D p kJ(u,u k)+12u−b 2for k=0,1,...,starting with u0=0and p0=0(hence,for k=0,one solves the original problem(2.8)).SinceμT V(u)is not differentiable everywhere,the subdifferential ofμT V(u) may contain more than one element.However,from the optimality of u k+1in(2.11),it follows that0∈∂J(u k+1)−p k+u k+1−b;hence,they setp k+1:=p k+b−u k+1.The difference between(2.8)and(2.11)is in the use of regularization.While(2.8)regularizes u by directly minimizing its total variation,(2.11)regularizes u by minimizing the total variation-based Bregman distance of u to a previous solution u k.In[73]two key results for the sequence{u k}generated by(2.11)were proved.First, u k−b converges to0monotonically;second,u k also gets closer to¯u,the unknown noiselessimage,monotonically in terms of the Bregman distance D p kT V (¯u,u k),at least while u k−b ≥¯u−b .Numerical results in[73]demonstrate that forμsufficiently large,this simple iterative procedure remarkably improves denoising quality over the original model(2.8).Interestingly,not only for thefirst iteration k=0,but for all k,the new problem(2.11) can be reduced to the original problem(2.8)with the input b k+1:=b+(b k−u k)starting with b0=u0=0;i.e.,the iterations(2.11)are equivalent to(2.12)u k+1←minu J(u)+12u−b k+1 2,where b k+1=b+(b k−u k),and can be carried out using any existing algorithms for(2.8).The iterative procedure(2.12)has an intriguing interpretation:Letωrepresent the noise in b,i.e.,b=¯u+ω,and letμbe large.At k=0,b k−u k=0,so(2.12)decomposes the input noisy image b into u1+v1.Sinceμis large,the resulting image u1is oversmoothed(by total variation minimization)so it does not contain any noise.Consequently,u1can be considered to be a portion of the original clean image¯u.The residual v1=b−u1=(¯u−u1)+ω,hence,is the sum of the unrecovered“good”signal(¯u−u1)and the“bad”noiseω.We wish to recover (¯u−u1)from v1.Intuitively,one would next consider letting v1be the new input for(2.8) and solving(2.8).However,Bregman iterative regularization turns out to be both better and “nonintuitive”:it adds v1back to the original input b.The new input of(2.12)in the second iteration isb+v1=(u1+v1)+v1=u1+2(¯u−u1)+2ω,which,compared to the original input b=u1+(¯u−u1)+ω,contains twice as much of both the unrecovered“good”signal¯u−u1and the“bad”noiseω.What is remarkable is thatBREGMAN ITERATIONS FOR 1-MINIMIZATION149 the new decomposition u2is a better approximation of¯u than u1(forμlarge enough);one explanation is that u2not only inherits u1but also captures a part of(¯u−u1),the previously uncaptured“good”signal.Of course,as the convergence results indicate,u k will eventually pick up the noiseωsince{u k}converges to b=¯u+ω.However,a high quality image can befound among the sequence{u k}:the image u k that has u k−b closest to ¯u−b is cleaner and has a higher contrast than the best image that could be obtained by solving(2.8)one single time,with the bestμ.Formally Bregman iterative regularization applied to the problem(2.13)minuJ(u)+H(u)is given as Algorithm1in which the Bregman distance D p kJ(·,·)is defined by(2.10).Algorithm1(Bregman iterative regularization).Require:J(·),H(·)1:Initialize:k=0,u0=0,p0=0.2:while“not converge”do3:u k+1←arg min u D p k J(u,u k)+H(u)4:p k+1←p k−∇H(u k+1)∈∂J(u k+1)5:k←k+16:end whileWe conclude this section by citing some useful convergence results from[73]that are used in section3.3.Assumption1.J(·)is convex,H(·)is convex and differentiable,and the solutions u k+1in step3of Algorithm1exist.Theorem2.1.Under Assumption1,the iterate sequence{u k}satisfies the following:1.Monotonic decrease in H:H(u k+1)≤H(u k+1)+D p kJ (u k+1,u k)≤H(u k).2.Convergence to the original in H with exact data:If˜u minimizes H(·)and J(˜u)<∞,then H(u k)≤H(˜u)+J(˜u)/k.3.Convergence to the original in D with noisy data:Let H(·)=H(·;f)and supposeH(˜u;f)≤δ2and H(˜u;g)=0(f,g,˜u,andδrepresent noisy data,noiseless data,perfect recovery,and noise level,respectively).Then D p k+1J (˜u,u k+1)<D p kJ(˜u,u k)aslong as H(u k+1;f)>δ2.3.Bregman iterations for basis pursuit.3.1.Formulations.The main purpose of this paper is to show that the Bregman iterative procedure is a simple but very efficient method for solving the basis pursuit problem(1.1),as well as a broader class of problems of the form(1.3),in both theory and practice.Below we first give the details of the algorithm,describe our motivation,and then prove that in afinite number of iterations,u k becomes a minimizer of u 1among{u:Au=f}.We solve the constrained problem(1.1)by applying Algorithm1to(1.2)for J(u)=μ u 1150W.YIN,S.OSHER,D.GOLDFARB,AND J.DARBONand H(u)=12Au−f 2:Version1:u0←0,p0←0,(3.1)For k=0,1,...dou k+1←arg minu D p kJ(u,u k)+12Au−f 2,(3.2)p k+1←p k−A (Au k+1−f);(3.3)Version2:f0←0,u0←0,(3.4)For k=0,1,...dof k+1←f+(f k−Au k),(3.5)u k+1←arg minu J(u)+12Au−f k+12.(3.6)Given u k and p k in Version1,u k+1satisfies thefirst-order optimality condition: 0∈∂J(u k+1)−p k+∇H(u k+1)=∂J(u k+1)−p k+A (Au k+1−f). Therefore,(3.7)p k+1=p k−A (Au k+1−f)∈∂J(u k+1);hence,D p k+1J (u,u k+1)is well defined.Clearly,if u k+1i=0,then one can pick any p k+1i∈[−1,1]and still have a well defined D p k+1J (u,u k+1).However,the choice of p k+1in(3.7)is notonly simple but also crucial for the sequence{u k}to converge to the minimizer u opt of the constrained problem(1.1).Theorem3.1.The Bregman iterative procedure Version1(3.1)–(3.3)and Version2(3.4)–(3.6)are equivalent in the sense that(3.2)and(3.6)have the same objective functions(up toa constant)for all k.Proof.Let u k and¯u k denote the solutions to Versions1and2,respectively.The initial-ization(3.1)gives D p0J (u,u0)=J(u),while(3.4)gives f1=f.Therefore,at iteration k=0,(3.2)and(3.6)solve the same optimization problem,min u J(u)+12Au−f 2.We note that this problem,as well as those for all other iterations k,may have more than one solution.We do not assume that in this case,u1(Version1)is equal to¯u1(Version2). Instead,we use the fact from[55]that A (f−Au)is constant for all optimal solutions u;i.e., A (f−Au1)=A (f−A¯u1).According to(3.3),p0=0,and f=f1,we havep1=p0−A (Au1−f)=A (f−Au1)=A (f−A¯u1)=A (f1−A¯u1).Next,we use induction on p k=A (f k−A¯u k).Given p k=A (f k−A¯u k),we will show the following:(i)the optimization problems in(3.2)and(3.6)at iteration k are equivalent,BREGMAN ITERATIONS FOR 1-MINIMIZATION 151(ii)A (Au k +1−f )=A (A ¯u k +1−f ),and (iii)p k +1=A (f k +1−A ¯u k +1).Clearly,part (i)proves the theorem.Part (i):From the induction assumption it follows thatD p k J (u,u k )+12 Au −f 2=J (u )− p k ,u +12 Au −f 2+C 1=J (u )− f k −A ¯u k ,Au +12 Au −f 2+C 2=J (u )+12 Au −(f +(f k −A ¯u k )) 2+C 3=J (u )+12Au −f k +1 2+C 3,where C 1,C 2,and C 3stand for terms constant in u ;hence,(3.2)and (3.6)have the same objective function (up to a constant).Part (ii):A (Au k +1−f )=A (A ¯u k +1−f )follows from part (i)and the result in [55].Part (iii):It follows from the induction assumption,as well as (3.3),(3.5),and part (ii),thatp k +1=p k −A (Au k +1−f )=p k −A (A ¯uk +1−f )(3.8)=A (f k −A ¯u k )−A (A ¯u k +1−f )=A f +(f k −A ¯uk )−A ¯u k +1 =A (f k +1−A ¯uk +1).(3.9)Remark .When J is not strictly convex,the subproblems in Versions 1and 2may both have more than one solution.The above proof shows,however,that even if Versions 1and 2generate different intermediate solutions at a certain iteration,they remain equivalent thereafter.Each iteration of (3.6)is an instance of (1.2),which can be solved by the code FPC[56].Although our convergence result below holds for any strictly positive μ,we choose μso that (1.2)is solved efficiently by FPC and the total time of the Bregman iterations is nearly optimal.3.2.Motivation.In [95],Xu and Osher applied Bregman iterative regularization to wavelet-based denoising.Briefly,they considered(3.10)min u μ u 1,1+12f −u 2L 2,where u 1,1is the Besov norm defined in [33];if u = j ˜u j ψj and f = j ˜f j ψj ,for a waveletbasis {ψj },they solved min {˜u j }μ j |˜u j |+12 j|˜f j −˜u j |2.It was observed in [95]and elsewhere that this minimization procedure is equivalent to shrink-age;i.e.,˜u j =shrink(˜f j ,μ),for all j ,where shrink(·,·)is defined in (2.4).152W.YIN,S.OSHER,D.GOLDFARB,AND J.DARBON What is interesting is that Bregman iterations gives(3.11)˜u k j=⎧⎪⎨⎪⎩˜fj,|˜fj|>μk−1,k˜f j−μsign(˜f j),μk≤|˜f j|≤μk−1, 0,|˜f j|≤μk.So soft shrinkage becomesfirm shrinkage[52]with thresholdsτ(k)=μk andτ(k−1)=μ(k−1).In[10,11]the concept of nonlinear inverse scale space was introduced and analyzed,whichis basically the limit of Bregman iteration as k andμincrease with kμ→t.This iterativeBregman procedure then approaches hard thresholding:(3.12)˜u j(t)=˜fj,|˜fj|>1t, 0,|˜f j|≤1t.For Bregman iterations it takes(3.13)k j=smallest integer≥μ|˜f j|iterations to recover˜u j(k)=˜f j for all k≥k j.This means that spikes return in decreasing order of magnitude and sparse data comes back very quickly.Next,we consider the trivial example of minimizing u 1subject to a u=f,where 0=a∈R n and f∈R.Obviously,the solution is u opt=(f/a j)e j,where e j is the j th unit vector and a j is the component of a with the largest magnitude.Without loss of generality, we suppose a≥0,f>0,and the largest component of a is a1>0,which is strictly larger than the rest(to avoid solution nonuniqueness);hence,u opt=(f/a1)e1.Let f k>0;then the solution of the Bregman iterative subproblemmin u μ u 1+12(a u−f k)2is given by(3.14)u k=0,μ≥f k a1,f k a1−μa21e1,0<μ<f k a1.The Bregman iterations(3.6)start with f1=f.Ifμ≥f1a1,then u1=0so f2=f+(f1−a 0)=2f;hence,as long as u i remains0,f i+1=(i+1)f.Therefore,we have u j=0and f j+1=(j+1)f for j=1,...,J,forJ=max{k:μ≥f k a1}=μfa1.Ifμ<f1a1,J=0.In both cases,u J+11=((J+1)fa1−μ)/a21sof J+2=f+(f J+1−a u J+1)=(J+2)f−a1(J+1)fa1−μa21=f−μa1andu J +2=f J +2a 1−μa 21e 1=f a 1e 1;i.e.,u J +2=u opt .Therefore,the Bregman iterations give an exact solution inμf max i {|a i |}+2steps for any problem with a one-dimensional signal f .We believe that these simple examples help explain why our procedure works so well in compressed sensing applications.3.3.Convergence results.In this section,we show that the Bregman iterative regular-ization (3.1)–(3.3)(or,equivalently,(3.4)–(3.6))described in section 3.1generates a sequence of solutions {u k }that converges to an optimum u opt of the basis pursuit problem (1.1)in a finite number of steps;that is,there exists a K such that every u k for k >K is a solution of (1.1).The analytical results of this section are generalized to many other types of 1and related minimization problems in section 5.We divide our analysis into two theorems.The first theorem shows that if u k satisfies the linear constraints Au k =f ,then it minimizes J (·)=μ · 1;the second theorem proves that such a u k is obtained for a finite k .Theorem 3.2.Suppose an iterate u k from (3.2)satisfies Au k =f ;then u k is a solution of the basis pursuit problem (1.1).Proof .For any u ,by the nonnegativity of the Bregman distance,we haveJ (u k )≤J (u )− u −u k ,p k (3.15)=J (u )− u −u k ,A (f k −Au k ) (3.16)=J (u )− Au −Au k ,f k −Au k (3.17)=J (u )− Au −f,f k −f ,(3.18)where the first equality follows from (3.9).Therefore,u k satisfies J (u k )≤J (u )for any u satisfying Au =f ;hence,u k is an optimal solution of the basis pursuit problem (1.1).Theorem 3.3.There exists a number K <∞such that any u k ,k ≥K ,is a solution of the basis pursuit problem (1.1).Proof .Let (I j +,I j −,E j )be a partition of the index set {1,2,...,n },and defineU j :=U (I j +,I j −,E j )={u :u i ≥0,i ∈I j +;u i ≤0,i ∈I j −;u i =0,i ∈E j },(3.19)H j :=min u 12Au −f 2:u ∈U j .(3.20)There are a finite number of distinct partitions (I j +,I j −,E j ),and the union of all possible U j ’s is H ,the entire space of u .At iteration k ,let (I k +,I k −,E k )be defined in terms of p k as follows:(3.21)I k +={i :p k =μ},I k −={i :p k =−μ},E k ={i :p k ∈(−μ,μ)}.。

Passage 3-01Labour should be ordered by 0900 hrs for 2nd shift (1500 to 2300 hrs) on same day and by 1300 hrs for 1st shift (0700 to 1500 hrs) for next day. By1100 hrs, for 3rd shift (2300 to 0500 hrs) in same day.Under normal circumstances, no work is performed during meal hours 1100 to 1200 hrs and 1830 to 1930 hrs. unless the ship is classified as a key vessel or the agent orders work during the meal hours.C 001. If you want to order labour in this port for the third shift in same day, the order should bemade by ________ .A. any timeB. 0900 hoursC. 1100 hoursD. 1300 hoursB 002.________is performed during meal hours under usual circumstances.A. Cargo workB. No cargo workC. Ordering work by AgentD. loading and dischargingA 003. Loading and discharging can be performed during meal hours _____A. if the ship is classified as a key vesselB. if agent orders work beforehandC. under normal circumstancesD. either A or BC 004. What does the word “key” mean ________?A. ladenB. smallC. pivotalD. bigPassage 3-02Before arrival in the United Kingdom, the master will have informed his owners or agents of the approximate time of the vessel's arrival at the pilot station for the port of destination. The vessel should be flying her ensign and also her signal letters and the requisite pilot signal when approaching the pilot station. The international signals, as well as any local port signals, can be found in the Sailing Directions, which is also known as the "Pilot Book"When a pilot is required most ports now require due notice of the vessel's ETA to be sent in by radio. However, this does not relieve the ship's obligation to display the pilot signal ("G" by any of the methods of signaling ) until the pilot is aboard when "H" flag will be flown. If the master or first mate of the vessel has a pilotage certificate for the district then the above is unnecessary, in such case the pilot flag (white and red horizontal halves, as on the pilot vessel) will be flown.B 005. While the pilot is on ship, ________should be displayed on the top of the ship mast.(1) "G" flag (2) "H" flagA. only (1)B. only (2)C. both (1) and (2)D. (1) plus (2)A 006. What's the meaning of the "ETA"________ .A. Estimated time of ArrivalB. Expected time of ArrivalC. Both A and BD. Neither A nor BA 007. Which statement under below is correct________ .A. The Sailing Directions Book is as same as Pilot BookB. The Sailing Directions Book and the Pilot Book are different booksC. The Pilot Book is a part of the Sailing Directions BookD. The Sailing Directions Book contains the Pilot BookA 008. Before arrival in the United Kingdom, the master should inform ________ tohis owner or agentsA. ETAB. ETDC. GMTD. LTPassage 3-03The basic concept of the GMDSS is that on-shore Search and Rescue authorities, in addition to shipping in the immediate vicinity of a ship in distress, will be alerted rapidly to an inc ident so that they can assist in co-ordinating a search and rescue operation with the minimum of delay. The system will also provide for urgency and safety communications, and the dissemination of Maritime Safety Information (MSI) including Navigational Warnings and Weather Messages. GMDSS applies to all ships of 300 G.R.T. or larger, to all passenger ships and all ships on international voyages, which are subject to the SOLAS Convention 1974, as amended 1988. GMDSS has been adopted by the IMO.A 009. What's the meaning of the SOLAS ________ .A. International Regulations for Preventing Collision At Sea.B. International Convention for the Safety of Life At SeaC. Salvage Operations of Life At SeaD. All above are wrongB 010. The "G.R.T." refers to ________ .A. Gross Register TonnageB. Gross Registered TonnageC. Gross TonnageD. Gross River TonnageA 011. GMDSS applies to ________ .A. all ships of 300 G.R.T. and over 300 G.R.T.B. all passenger shipsC. all ships on international routeD. all above are correctC 012. The full name of GMDSS is ________ .A. Global Maritime Distress Satellite SystemB. Global Maritime Distress and Safety SystemC. Global Marine Distress Satellite SystemD. Global Marine Distress and Safety SystemPassage 3-04The axial thrust of the propeller is the force working in a fore and aft direction. This force causes the ship to move ahead through the water or to go astern. Because of her shape，a ship will move ahead through the water more easily than going astern.The transverse thrust is the sideways force of the propeller as it rotates. The transverse effect of the propeller blades at the top near the surface of the water is not strong enough to counteract the opposite effect of the lower blades. For right-handed propellers this cants the ship's stern to starboard and her bow to port，when the ship is going ahead. The effect is small and can be corrected by the rudder. When the engines are put astern，the effect is the opposite and the stern cants to port. This effect is stronger and cannot easily be corrected. V essels with left-handedpropellers behave in the opposite way.A 013.The force that causes the ship to move ahead through the water or to go astern is known as________.A. axial thrustB. transverse thrustC. the transverse effect of the propeller blades at the top near the surface of the waterD. the transverse effect of the lower blades of the propeller near the bottom of the waterB 014.A left-handed propellers，when the ship is going ahead，will cant ship's stern to________.A. starboardB. portC. to move aheadD. move asternA 015.The transverse thrust of the propeller is stronger when the ship is________.A. going a sternB. going aheadC. stoppedD. making no way through the waterA 016.The transverse thrust of the propeller can mainly be overcome by ________.A. the rudderB. the propeller itselfC. the nautical instrumentD. wind and tidePassage 3-05The task of a helmsman is to steer the ship precisely according to the instructions of the officer of the watch. On passage, these instructions take the form of a course to be held, regardless of wind, sea or other sources of deviation. At other times, the instructions will call for the ship to be turned to port or starboard, or to be prevented from turning, or to be lined up on a particular heading.A good helmsman uses the least rudder deflection to maintain a steady course. He must learn quickly at the beginning of each trick whether the ship is carrying port or starboard helm.Most ships now set and steer courses in degrees on the full scale of 0-360degrrs. This system replaces the traditional use of the compass scale divided into cardinals and points. Nevertheless, every seaman and navigator must be familiar with the older system, since the necessary for obtaining certain qualifications (e.g. the Certificate of Efficiency as lifeboatmen); points are also used for indication wind direction.C 017. Nowadays, the “helmsman” is also called ________ onboard.A. Ordinary seamanB. BosunC. Able seamanD. CarpenterB 018.What was used for indication the ship's course in traditional navigation________ .A. Adegrees on the full scale of 0-360B. cardinals and points of compassC. sounding in metersD. distance in milesA 019. When dose the helmsman has to steer the ship precisely ________ .A. according to the instructions of the officer of the watchB. when the ship to be turned or to be prevented from turningC. when the ship to be lined up on a particular headingD. all above are correctB 020. According to "points are also used for indication wind direction", please predict how toindicate the set of the current________ .A. degreesB. pointsC. centigradeD. Beaufort scalePassage 3-06While every effort is made to ensure that the data provided through the Notices to Mariners service is accurate，the user needs to be aware of the risks to corruption of data. It is important that the user should only use the data on suitable equipment and that，other applications should not be running on the user's machine at the same time. Users should exercise their professional judgement in the use of data，and also consult the Mariners Handbook (NP100) for further details. The user needs to be aware that there is a possibility that data could be corrupted during transmission，or in the process of display or printing on the user's equipment，or if converted to other software formats，and is accordingly advised that the UKHO cannot accept responsibility for any such change，or any modifications or unauthorised changes，made by licensees，or other parties.D 021.The data may become corrupted in any of the following process except _______.A. during transmissionB. in the display or printing on the user's equipmentC. in converting to other software formatsD. in air mail delivery to the readersA 022.The use of the data is advised to consult _______ for further details.A. Mariners HandbookB. Sailing DirectionsC. Guide to Port EntryD. Notices to MarinersD 023.Of the following items _______ is not mentioned for which UKHO will accept noresponsibility.A. change in the process of display or printingB. unauthorised changes made by licensees or other partiesC. modifications made by licensees or other partiesD. professional amendmentsC 024.It is implied that _______.A. the data are incorrectB. the data are to be corrected intensivelyC. although the data are accurate enough，you are still advised to use it with cautionD. not to use it if you have not enough time or proper equipment to effect necessarycorrectionPassage 3-07The container ship is different from the conventional type and is an innovation noted for easier handling and quicker turnover of cargoes. Cargoes to be carried by this type of ship are pre-packed into containers before being loaded aboard the ship.Containers are sealed after being packed with cargoes. Made of metal or other durable materials，they are watertight after sealing and can therefore be stowed on deck whilst being carried. One of the features of container ships is that some of the containers are usually stowed on deck.The container ship is becoming increasingly popular in trading circles，and the trend is that the tonnage thereof will grow at a faster pace in future.C 025.What does "innovation" in the first paragraph mean? ________.A. making changesB. the introduction of an antigenic substance into the body against a specific diseaseC. The act of introducing something new.D. revolutionD 026.Containers are sealed after being packed with cargoes.A. filledB. loadedC. stuffedD. closed officially or under the supervision of notary publicC 027.Of the following，________ is not the feature of the container ship?A. Some of the containers are usually stowed on deck.B. It is easy for handling and quick turnover of the cargoC. The container ship is becoming increasingly saferD. Cargoes are pre-packed into the containerB 028.The tonnage of container ship is ________.A. decreasingB. increasingC. remaining the sameD. changingPassage 3-08Nautical charts are indispensable to mariners. They，however，are subject to frequent changes，such as those of navigational aids，of waterways due to the dredging and construction，of depths of water，and of removal or appearance of wrecks. In order to keep up-to-date and reliable，nautical charts have to undergo correction. Changes of importance are generally promulgated by weekly edition of Notices to Mariners，which enable mariners to correct the charts by hand. If major changes make it impracticable to do so，the Notices will provide a reproduction of a small area，which is also called block，to be pasted onto the chart in its correct position.C 029.Nautical charts need correction because ________.A. navigational aids are sometimes indispensable.B. there are always some mistakesC. wrecks may appear or be removedD. they could never be reprintedA 030.Correction to charts are made by crew members in accordance with ________.A. Notices to MarinersB. Sailing DirectionsC. Guide to Port EntryD. SupplementC 031.In the passage，Blocks are ________.A. large scale chartsB. representations of chartsC. reproductions of portions of chartsD. small scale chartsA 032.The purpose of correction to charts is to ________.A. keep them up-to-dateB. make the charts brand-newC. keep the charts available to all mariners in the worldD. keep the charts free from mistakesPassage 3-09Corrections to Sailing Directions are given in Section Ⅳ. Those in force at the end of the year are reprinted in the Annual Summary of Notices to Mariners. A list of corrections in force is published in Section Ⅳ of the Weekly Edition for the last week of each month.It is recommended that corrections be kept in a file with the latest list of corrections in force on top. The list should be consulted when using the parent book to see if any corrections affecting the area under consideration are in force.It is not recommended that corrections be stuck in the parent book or current supplement，but，if this is done，when a new supplement is received care must be taken to retain those corrections issued after the date of the new supplement，which may be several months before its receipt on board.B 033.________ are reprinted in the Annual Summary of Notices to Mariners.A. The Sailing DirectionsB. The corrections to Sailing DirectionsC. The effective corrections to Notices to MarinersD. The Weekly EditionA 034.The parent book is ________.A. The Sailing DirectionB. The corrections to Sailing Directions in forceC. the Annual Summary of Notices to MarinersD. the Weekly EditionD 035.It is recommended that corrections to the Sailing Directions be ________.A. made by handB. consulted at the last week of each monthC. stuck in the parent book or current supplementD. kept in a file with the latest list of corrections in force on topA 036.If the corrections be stuck in the parent book or current supplement，________.A. when a new supplement is received，those corrections issued after the date of the newsupplement must be retainedB. the parent book must be consultedC. the current supplement must be consultedD. the Annual Summary of Notices to Mariners must be usedPassage 3 -10The amount of detail shown on a chart varies with the scale of the chart. On a large scale chart，for example，full details of all lights and fog signals are shown，but on smaller scales the order of reduction of information in elevation，period，range，until on an ocean chart of the area only lights with a range of 15 miles or more will normally be inserted，and then only their light-star and magenta flare. On the other hand，radio beacons are omitted from large scale charts where their use would be inappropriate，and，unless they are long range beacons，from ocean charts.B 037.Ocean charts are ________ ones.A. large scaleB. small scaleC. inappropriateD. omittedA 038.What cannot be found in the large scale charts? ________.A. Radio beacons of small rangeB. Full details of all lights.C. ElevationsD. Full details of fog signalsC 039.The light-star and magenta flare are shown on ________.A. large scale charts onlyB. small scale charts onlyC. both small and large scale chartsD. neither small nor large scale chartsA 040.The title of this passage should be ________.A. Lights and Beacons on ChartsB. Characteristic of Lights and BeaconsC. Corrections to Small and Large Scale chartsD. Navigational Charts PublicationPassage 3-11DALIAN OBSY GALE WARNING 190600ZCOLD FRONT WILL PASS BOHAI SEA BOHAI STRAITS NORTH AND CENTRAL HUANGHAI SEA CAUSING GALE WINDS TOMORROW AFTERNOON AND EVENING STOP.SYNOPTIC SITUA TION 190600ZLOW 994 HPA A T 48N 118E MOVING SE 8 KTS WITH COLD FRONT FROM CENTER PASSISNG 44N 128E HIGH 1013HPA A T 38N 124E STA TIONARY STOP24HOURS WEA THER FORECAST FROM 191000ZBOHAI SEA BOHAI STRAITS NORTH AND CENTRAL HUANGHAI SEA PARLY CLOUDY BECOMING OVERCAST TOMORROW WITH RAIN SW WINDS FORCE 7 TO 8 TOMORROW A TERNOON AND EVENING SEA ROUGH BECOMING VERY ROUGH STOP.A 041.The COLD FRONT will pass Bohai Sea，Bohai Straits，North and central Huanghai Sea onA. The 20thB. The 19thC. The 18thD. The 6thB 042.________ is stationary at 38N 124E.A. Low 994 HpaB. High 1013 HpaC. Cold frontD. Warm frontD 043.The winds are expected tomorrow to be_______?A. roughB. very roughC. SE 8 knotsD. SW 7-8 in forceA 044.What is the weather like tomorrow in this area? _______.A. It will be partly cloudy becoming overcast with rain and SW force 7-8 windsB. LOW 994 HPA at 48N 118E is moving SE 8 KTS with COLD FRONT from centerpassing 44N 128EC. HIGH 1013HPA at 38N 124E will be stationaryD. It will rain the whole dayPassage 3-12In some parts of a chart where the spaces are rather blank and there are no symbols of any kind，there may be Cautions，Warnings，Notes，etc.，which should be taken into account while using a chart. All of those Cautionary Notes give the mariner facilities to ensure safe navigation，such as to avoid running aground in shallow waters and making damages to nearby fishing gears，and to keep off any hazards in areas where submarine frequently exercises. Furthermore，they are of goodhelp to mariners，as to the reliability of the navigational aids especially in congested waters or narrow channels，to prevent any possible accidents.D 045.What is the main topic of this passage? ________.A. Regulations of the harborB. Details in the Sailing DirectionsC. Rules of the terminalD. Description on Admiralty ChartsA 046.According to the passage，you must pay attention to ________ while using a chart.A. Cautions，Warnings and NotesB. Reports，Symbols and ChartsC. Explanations，accounts and answersD. Damages，hazards and injuriesC 047.Cautionary Notes are helpful for mariners ________.A. to run aground in shallow watersB. to make damages to nearby fishing gearsC. to keep off hazards in areas where submarine exercisesD. to keep the reliability of the aids to navigation in congested waters or narrow channelsD 048.Cautions，Warnings，Notes，etc. are likely inserted in some parts of a chart where ________.A. submarine frequently exercisesB. there are fishing gearsC. the waters is congested and the channels are narrowD. the spaces are rather blank and there are no symbols of any kindPassage 3-13Logbooks required by law，to be filled out by masters or officers on duty of every ship，the forms of which must be proved by the shipping companies or marine authorities.Logbooks are used to record the events occurring during the ship's stay in a harbor，at anchorage，or underway，and they are also requested to produce evidences in case officials inquire about accidents.On completion of the voyage the logbook must be submitted to the superintendent of the owner or the marine authorities for justification，checking or approval. Therefore，everything recorded in the logbook must be true and accurate.When a misentry has been made in the log，a red line would be drawn on those parts. The correct entry with signature should be made near or above them. No erasures or cuts are to be allowed.B 049.The best title for the passage is " ________ "A. The forms of logbooksB. The use of logbooksC. Characteristics of logbooksD. How to check logbooksD 050.When a misentry has been made in the log，________.A. erasures or cuts are to be allowed.B. it is to be corrected out by masters or officers on duty of every shipC. it is to be produced in case officials inquire about accidents.D. a red line would be drawn on those parts，with correct entry with signature being madenear or above them.B 051.The forms of logbooks must be proved by ________.A. officials who inquire about accidents.B. the shipping companies or marine authorities.C. masters or officers on duty.D. the superintendent of the owner.A 052.The logbook must be submitted ________ to the superintendent of the owner or the marineauthorities for justification，checking or approval.A. on completion of the voyageB. in a harborC. at anchorageD. underwayPassage 3-14For navigation，radar is of incredible value. It provides the navigator with his position，his distance from ships or obstructions nearby and other accurate information to prevent collision and ensure the safety of the ship. Radar can display all objects within its working range clearly，either in clear weather or in thick fog. In addition，if the radar information is correctly interpreted，the navigator can easily work out the speed and direction of an approaching object and take proper measures to keep his ship from any danger.Shore-based radar also plays an important role in shipping. If ship's radar is in trouble，the radar observer at the stations will use VHF radio to alert them to other traffic in the vicinity as well as to advise their position. Up to now，many radar surveillance systems have been installed in most large seaports. They are intended to smooth and control the flow of traffic to and from the harbor.B 053.For navigation，the radar is ________.A. of no valueB. very importantC. so expensive that people don't know how much it isD. valuelessC 054.Which of the following statements about radar's function for marine purposes is incorrect?________.A. It provides the navigator the ship's positionB. It provides information to protect ships from collisionC. It displays all the objects at sea clearlyD. It displays the observer's distance from ships and obstructions nearbyA 055.If the ship's radar is in trouble，the shore-based radar ________.A. may provide the ship of her positionB. should be installed with surveillance systemsC. shall advise the ship to use VHFD. will be put into use immediatelyD 056.Radar surveillance systems ________.A. may provide all ships of their technical conditionsB. should be installed with VHFC. shall be correctly interpreted，D. are intended to smooth and control the flow of traffic to and from the harbor.Passage 3-15Communications over relatively short distances can be made by visual or sound signals. V isual signals can be sent by using flags or an Aldis lamp. An Aldis lamp is an electric lamp used for flashing messages in Morse code. The traditional method of signaling from one ship to another is by using flags. There are different colored flags for each letter of the alphabet. There are alsopennant-shaped flags for numbers，and a long pennant，known as an answering or code pennant. Three other flags，which are burgee-shaped，are known as substitutes. These show that the flat or pennant is being repeated. Besides standing for a letter of the alphabet，each flag，when hoisted along，has another meaning. For example，the "W" flag also means: "I require medical assistance". Flags can also be hoisted in combinations of two，three or four. Siren，whistle，bell or other sound signals can be used in fog and similar circumstances when visual signals can not be seen.D munications over relatively short distances may be made by ________.A. visual signalsB. sound signalsC. Morse CodeD. Either visual or sound signalsA 058.An Aldis lamp is used for ________.A. transmitting Morse codeB. flashing flagsC. sending flag signalsD. sending sound signalsA 059.Burgee-shaped flags are used as substitutes to show ________.A. "repeating"B. "answering"C. "code" pennantD. "I requiring medical assistance"D 060.________ are used in fog and similar circumstances when visual signals can not be seen.A. V isual signalsB. SubstitutesC. Pennant-shaped flagsD. The ship's siren，whistle or bellPassage 3-16When the senders of goods have large shipment to make，and especially when bulk cargo is concerned，it is advisable that they have some ships at their disposal. Some of the big companies set up a fleet of their own，but the rest may find it more profitable to hire instead of building or buying ships. This is called "chartering". The chartering of the ship is usually done through the intermediary of brokers，who，when hired，will go through all the necessary formalities on behalf of the charterer. In London there is a special center "the Baltic Exchange"，where the brokers operate in much the same way as stock and share brokers on a stock exchange. But it is easy for home shippers to hire Chinese or foreign ships through China National Chartering Corporation，which takes care of chartering business on orders from various import and export corporations.D 061.When large shipment is concerned，________ is not the way for the sender to have ships attheir disposal.A. to charter shipsB. to build shipsC. to buy shipsD. to scrape shipsB 062.In chartering all the necessary formalities are performed through ________.A. the intermediary of agentsB. the intermediary of brokersC. the charterersD. the "Baltic Exchange"D 063.The function of "the Baltic Exchange" is ________.A. to deal with stocksB. to exchange cargoesC. to operate on sharesD. to charter shipsD 064.China National Chartering Corporation takes care of chartering business for home shippers."To take care of " means ________.A. to pay attention toB. to be concerned withC. to be liable forD. to take charge ofPassage 3-17A tropical storm is not so extensive as the depression of higher latitudes but，within 75 miles or so of the center，the wind is often far more violent，and the high and confused seas near the center may cause considerable damage to large and well-found ships，while small vessels (for example，destroyers) have foundered. The danger is still greater when ships are caught in restricted waters without adequate room to maneuver. Within 5 to 10 miles of the center the wind is light or moderate and variable，the sky is clear or partially so，and there is a heavy，sometimes mountainous，confused swell. This area is known as the "eye" of the storm. After passing through the relatively windless center of the storm the wind will suddenly，and with great violence，commence to blow from a direction opposite to that experienced on the other side of the windless center. Due to torrential rain visibility near the storm center is almost nil.D 065.Within ________ of a tropical storm center，the wind is violent.A. no more than 75 milesB. not more than 75 milesC. 75 miles or a greater distanceD. about 75 milesD 066.Among the following，________ one may not be found in the "eye" of the storm?A. The visibility is moderate or goodB. The wind is light or moderateC. The sky is clear or partly cloudyD. The swell is low or moderate.C 067.In the passage，"a well-found ship " means ________.A. a ship has been found in any placeB. a ship has been found in good visibilityC. a ship with all the necessary equipment properly maintainedD. a ship in huge sizeA 068.The visibility near tropical storm center is ________.A. V ery poorB. PoorC. ModerateD. GoodPassage 3-18By turning the GAIN control clockwise，the gain of the receiver increases and the observing range of the target expands. Adjust this control so that the best pictures may be displayed on the screen，according to the range scale in use. In the short range，it is advisable to operate the equipment with this control set at a setting where the receiver gain is rather lowered a little. In the long range，it is advisable to operate the equipment with this control set at a setting where the receiver gain is rather increased a little. With too little gain，the small targets are missed and there is a decrease in the detected range. With excessive gain，since the screen becomes brighter because the noise increases，the contrast between echoes and background noise reduces，making target observation more difficult. In the crowded regions，the gain may be reduced to clear the picture.A 069.Switching from short range to long range，you will have to _______.。

老托福阅读100篇passage33试题及答案为了帮助大家备考托福阅读，提高成绩，下面小编给大家带来老托福阅读100篇passage 33试题及答案，希望大家喜欢!老托福阅读100篇passage 33试题及答案PASSAGE 33Researchers in the field of psychology have found that one of the best ways to make an important decision, such as choosing a university to attend or a business to invest in, involves the utilization of a decision worksheet. Psychologists who study optimization compare the actual decisions made by people to theoretical ideal decisions to see how similar they are. Proponents of the worksheet procedure believe that it will yield optimal, that is, the best decisions. Although there are several variations on the exact format that worksheets can take, they are all similar in their essential aspects. Worksheets require defining the problem in a clear and concise way and then listing all possible solutions to the problem. Next, the pertinent considerations that will be affected by each decision are listed, and the relative importance of each consideration or consequence is determined. Each consideration is assigned a numerical value to reflect its relative importance. A decision is mathematically calculated by adding these values together. The alternative with the highest number of points emerges as the best decision.Since most important problems are multifaceted, there are several alternatives to choose from, each with unique advantages and disadvantages. One of the benefits of a pencil and paper decision-making procedure is that it permits people to deal with more variables than their minds can generally comprehend andremember. On the average, people can keep about seven ideas in their minds at once. A worksheet can be especially useful when the decision involves a large number of variables with complex relationships. A realistic example for many college students is the question What will I do after graduation? A graduate might seek a position that offers specialized training, pursue an advanced degree, or travel abroad for a year.A decision-making worksheet begins with a succinct statement of the problem that will also help to narrow it. It is important to be clear about the distinction between long-range and immediate goals because long-range goals often involve a different decision than short-range ones. Focusing on long-range goals, a graduating student might revise the question above to What will I do after graduation that will lead to successful career?1. What does the passage mainly discuss?(A) A tool to assist in making complex decisions.(B) A comparison of actual decisions and ideal decisions(C) Research on how people make decisions(D) Differences between long-range and short-range decision making2. The word essential in line 7 is closest in meaning to(A) introductory(B) changeable(C) beneficial(D) fundamental3. The word pertinent in line 9 is closest in meaning to(A) relevant(B) preceding(C) insightful(D) responsive4. Of the following steps, which occurs before the others in making a decision worksheet?(A) Listing the consequences of each solution(B) Calculating a numerical summary of each solution(C) Deciding which consequences are most important(D) Writing down all possible solutions5. According to decision-worksheet theory, an optimal decision is defined as one that(A) has the fewest variables to consider(B) uses the most decision worksheets(C) has the most points assigned to it(D) is agreed to by the greatest number of people6. The author develops the discussion in paragraph 1 by means of(A) describing a process(B) classifying types of worksheets(C) providing historical background(D) explaining a theory7. The author states that On the average, people can keep about seven ideas in their minds atonce (lines 17-18) to explain that(A) most decisions involve seven steps(B) human mental capacity has limitations(C) some people have difficulty making minor as well as major decisions(D) people can learn to keep more than seven ideas in their minds with practice8. The word succinct in line 24 is closest in meaning to(A) creative(B) satisfactory(C) personal(D) concise9. Which of the following terms is defined in the passage ?(A) Proponents (line 5)(B) Optimal (line 5)(C) Variables (line 17)(D) Long-range goals (line 25)10. The word it in line 24 refers to(A) worksheet(B) problem(C) distinction(D) decision11. The word revise in line 26 is closest in meaning to(A) ask(B) explain(C) change(D) predictPASSAGE 33 ADADC ABDBB C托福阅读怎么抓住定位词首先介绍一下，什么是定位词?其实很简单，打个比方，你和朋友约好了去酒吧，朋友和你说酒吧在沈阳新东方正对面，这个酒吧你是不知道地点的，也就是你的目的地;而新东方却很熟知，那么你只需找到新东方便可以找到酒吧了。

高二英语数学建模方法单选题20题答案解析版1.In a mathematical modeling project, we need to find the optimal solution. The word "optimal" means _____.A.goodB.bestC.acceptablemon答案：B。

“optimal”的意思是“最佳的”，选项A“good”是“好的”，不如“最佳的”准确；选项C“acceptable”是“可接受的”，也不符合；选项D“common”是“常见的”，与“optimal”意思相差甚远。

2.When we build a mathematical model, we often use some assumptions. The word "assumptions" means _____.A.factsB.guessesC.proofsD.results答案：B。

“assumptions”的意思是“假设”，比较接近“guesses”（ 猜测）；选项A“facts”是“事实”；选项C“proofs”是“证明”；选项D“results”是“结果”，都与“假设”意思不同。

3.In a mathematical modeling process, we need to validate the model. The word "validate" means _____.A.testB.createC.changeD.ignore答案：A。

“validate”的意思是“验证”，与“test”（ 测试）意思相近；选项B“create”是“创建”；选项C“change”是“改变”；选项D“ignore”是“忽略”，都不符合“验证”的意思。

4.A mathematical model should be accurate and reliable. The word "reliable" means _____.efulB.trustworthyC.interestingD.difficult答案：B。

the z-score and the standard normal table第一节讲了probability的概念和隐藏的一个前提，random sampling。

random sampling又分为random sampling with replacement (independent random sampling)和random sampling without replacement。

最后讲了probability和frequency distributions的关联，就在于，probabilities and proportions are equivalent。

第二节开始，在之前（第二章）作为一个宋兵甲出现过的小角色normal distribution开始隆重登场了。

因为normal distribution是可以用方程来表示的，还可以求面积。

而面积意味着可以计算分布其中的数值的比例，比例意味着probability的计算。

这里又有一个现成的表格，the unit normal table，可以辅助你计算面积或是概率。

The unit normal table lists relationships between z-score locations and proportions in a normal distribution…Similarly, if you know th proportions, you can use the table to find the specific z-score location.具体细节不再多说。

第三节是深入第二节的内容，各种计算问题。

前三节实际上都是理论，第四节讲的是应用。

probability, proportion, normal distribution, the unit normal table的具体应用，体现在binomial distribution上。

Shelling Procedure and Optimization by Simulated Annealing For Sphere Packing SummaryWe provide two models (model A and model B) for Gamma Knife treatment planning, one is a Mixed Integer Programming (MIP) model, and the other is an optimal sphere-packing model.Based on dose distribution and the requirements of Gamma Knife unit, a Mixed Integer Programming (MIP) model is constructed and discussed.Another model is a sphere packing approach to Gamma Knife Treatment Planning problem. Based on the theory of digital image processing and simulated annealing, an algorithm to obtain an initial sphere-packing plan is designed, and a method from the thought of simulated annealing is used to optimize the solution.Due to the complexity of the MIP model and the large variables, the calculating would cost too much time. Many researchers have developed or have been developing other programming models. There are the same problems of time complexity. So we have to establish model B to speed up the computation and find the best or approximate best plan for the Gamma Knife Treatment.In model B, we simulate some images from CT/MRI, and then translate them into digital-image as matrices by image processing. Then we get an initial result through Shelling Procedure and adjust further.Our algorithm is recursion procedure for sphere packing. It determines the candidate location of the sphere center heuristic from largest sphere radius to the smallest one. As we know, when a ball moves inside a large close region, the track of the sphere center forms a close surface with the shape similar to the boundary of the region. We adopt heuristic method to choose one voxel as the center of the sphere. The idea of finding the largest distance between two voxels is that we should use the space efficiently. The spheres located in these two points can reduce the waste voxels near the boundary, and can spare more space for further utilization.The solution obtained by Shelling Procedure must be adjusted further. From the thought of Simulated Annealing, we realize a program to maximum the sum of volume of spheres and to minimum the number of them.Simulated computation shows model B satisfies the constraints of the Gamma Knife unit treatment planning and can be further studied to use in the real systems.BackgroundGamma Knife Unit[6,7,8,9,10]The Gamma Knife is a highly specialized treatment unit that provides an advanced stereotactic approach to the treatment of tumor and vascular malformations within the head. The Gamma Knife delivers a single, high dose of gamma ray emanating from 201 Cobalt-60 unit sources. Inside a shielded treatment unit, beams from 201 cobalt-60 radioactive sources are focused so that they intersect at a certain point in space, producing an ellipsoidal region of high radiation dose referred to as a shot. (Figure 1)Figure 1: A shot of radiation is formed at the intersection of 201 beamsA brief historyIn 1968, Professor Lars Leskell of the Karolinska Institute in Stockholm, Sweden and Professor Borge Larsson of the Gustaf Werner Institute at the University of Uppsala, Sweden developed the Gamma Knife. As far back as the 1940's, Leskell recognized the need for an instrument to target deep-seated intracranial structures without the risks of invasive open skull surgery. Currently, there are about 200 Gamma Knife machines worldwide.Treatment Procedure1.Fix patient’s head.e “magnetic resonance imaging” (MRI) or “computed tomography” (CT) scan thepatient’s head to find the location and the volume of the tumor.3.Develop the patient's treatment plan. Find a optimal multiple shots plan due to theirregularity and size of tumor shapes and the fact that the focusing helmets are only available in four sizes (4, 8, 14 and 18mm).4.Deliver an efficient high dose of radiation to the target volume.Treatment GoalThe plan aims to deliver a high dose of radiation to the intracranial target volume with minimum damage to the surrounding normal tissue. The treatment goals can vary from one neurosurgeon to the next, so a planning tool must be able to accommodate several differentrequirements. Among these requirements, the following are typical, although the level of treatment and importance of each may vary.1. A complete 50% isodose line coverage of the target volume. This means that thecomplete target must be covered by a dose that has intensity at least 50% of the maximum delivered dosage. This can be thought of as a “homogeneity” requirement. 2.To minimize the nontarget volume that is covered by a shot or the series of deliveredshots. This requirement is clear and can be thought of as a “conformity” requirement. 3.To limit the amount of dosage that is delivered to certain sensitive structures close to thetarget. Such requirements can be thought of as “avoidance” requirements.In addition to these requirements, it is also preferable to use as small number of shots as possible to reduce the treatment time.Problem AnalysisThe goal of stereotactic radiosurgery for a brain tumor is to deliver the desired dosage to the target, and only the target. This is not possible in reality. So they do the next best thing, which is to deliver enough dosage to the target, to avoid as much normal tissue as possible, and to deliver as little radiation as possible to whatever normal tissue must be affected. There are two additional important criteria–dose homogeneity and dose conformality. That is, we do not want ‘hot spots,’ which have been experimentally determined to cause complications; and we do want rapid falloff of dose levels outside the actual tumor. As for simplification, we can firstly consider the focused beams an ideal sphere solid. That means, the dose inside the sphere is homogeneous and the dose level of sphere of different size is the same. And the dose outside the sphere is 0, take D=4mm sphere for example (Figure 2).Figure 2: diameter of 4mm sphere dose distribution curveThe figure shows that a shot merely acts upon in the spherical region; the outside of the sphere is not affected by the shot. Therefore, the target volume can be filled with several shots, and the problem can be reduced to sphere packing problem.However, the dose distribution of a shot in practice is not a simple step function with 0 dose outside the sphere and the dose inside the sphere is not homogeneous, that is, the dose gradient exists. In this situation, we can reduce the target volume to a smaller size and we also use sphere packing plan to solve this problem.Basic Assumptionsz 201 beams of radiation simultaneously intersect at the isocenter forming an idealspherical dose distribution.z A shot is a non-elastic, solid 3D sphere.z The shot sphere diameters are the same size as the diameters of collimator (diameter of 4,8, 14 and 18mm).z The tumor volume is not large. The length, width and height of a tumor region are from20mm to 40mm.z The dose distribution is that the dose level inside the shot sphere is high enough to killcancer cells, and the dose outside is low enough.z A treatment plan is considered acceptable if 50% isodose curve encompasses the target. z We consider isodose curve encompasses the target region tightly and closely.Constructing the ModelsModel A: Model Based on Mixed Integer ProgrammingDose distribution modelLet x be the distance from the isocenter of a dose sphere, r be a measure of the radius of the sphere, A sum of error functions has been noted in the literature to approximate this dose distribution. Then the radius dose distribution for x can be expressed as:∑=−⋅−n i ii i r x erf m 1(1σ (11=∑=n i i m , n =2 for simply) where ∫∞−−=xt dt e x erf 2/221)(π. According the above expression, we could draw the images that dose distribution (Figure 3)and effective area of dose (Figure 4) approximately.Figure 3: Dose distribution Figure 4: Effective dose distribution areaDescription of ModelThe complete dose distribution can be calculated as a sum of contributions from each shot delivered, once the location of the center of that shot (s s s z y x ,,) is known, and the length of time of delivery s t ; w is known as four kinds of diameter. In practice, this means that for all (i, j, k ), we restrict the shot locations to be within the target area, set shot location on grid, (we generate grid on target that values 1mm (Figure 5), and use binary variables to indicate if a pair of (shot location, shot size) is used or not. We choose a grid of possible shot location, and pre-calculate for each grid location, every pixel and each r and use the optimization algorithm to decide whether or not to use a shot as a particular location. Since the only shots that are considered are shots that lie within the target, so we could easily determine whether or not a shot lies within the target by the terms of track bounds of target. When the description of the dose is determined an optimization model can be formulated. The basic variables of the optimization include the number of shots of radiation that will be delivered, along with the coordinates of the center location of the shot (s s s z y x ,,), the dose of any point in the target ),,(k j i D s , the time s t that each shot is exposed, the T of bound on the length of time that a particular shot can be exposed, and s ψ(heuristically. If use the shot s, then, s ψ=1, else s ψ=0).Figure 5: Grid on figureNow, the problem is reduced to solving a MIP (Mixed Integer Programming ):),,(:),,,,,(}1,0{2),,(1:.),,(),,(:k j i D k j i z y x D N t k j i Dose t S k j i D t k j i Dose Min s s s s s s S s s s s s s T T =∈≤≤≤≤≤∗=∑∑∈−−ψψψψDue to the complexity of the MIP model and the large variables, the calculating would cost too much time. Other programming models were constructed by many researchers [6,7,8,9,10], there were the same problems of time complexity. So we have to establish another model to simplify the problem.Model B: Model Based on Shelling Procedure andOptimization by Simulated AnnealingDescription of ModelUnder the assumption that a shot is a non-elastic, solid 3D sphere, the radiosurgical treatment planning can be deduced an optimization of packing unequal spheres into a three-dimensional bounded region. Given an input (R; V; S; L), where R is a 3D bounded region, V a positive integer, S a multiset of spheres, and L a location constraint on spheres. We want to find a packing of R using the minimum number of spheres in S such that the covered volume is at least V, the location constraint L is satisfied; and the number of points on the boundary of R that are touched by spheres is maximized. Wang [14] shows that not only finding an optimal solution to the problem is computationally intractable, but also optimization of the related problems is NP-hard. Therefore, some sort of approximation is needed. The paper [2,3] proposes a model under the assumption that spheres are available with unlimited supply, the 3D bounded region is a polytope, and there are no location constraints. The model is a nonconvex optimization problem with quadratic constraints and a linear objective function. The computation complexity of the model of [3] is very high.Shelling Procedure to produce an initial solutionOur treatment planning for a gamma knife unit is based on sphere-packing problem. We establish a shelling procedure to solve sphere-packing problems.Before we depict our procedure, it is necessary to know the following principles:The sequence of packing is from the sphere of the largest diameter to the smallest one. Definition: Considering a voxel as a three-dimensional box. The length, width and the height of each voxel is 1unit.The major steps of this algorithm are as follows:Step1:Transfer the data of the boundary of the target to three-dimensional 0-1 arrays. The voxel value in the target region is 1, and the voxel value outside the region is 0.Step2: Find a reasonable sphere center.Pack from the largest sphere to the smallest one. First, we pack the sphere of diameter D=18mm into the region. As we know, when a small ball moves inside a large close region, the track of the sphere center forms a close surface with the shape similar to the boundary of the region. However, the shape similarity is small as the diameter increases.The method to obtain the close surface of a sphere center is similar to the methods in [1,7]. Find two voxels on the close surface, the distance between which is the largest. We choose one voxel as the center of the sphere in a heuristic way.The idea of finding the largest distance between two voxels is that we should utilize the space efficiently. The spheres located in these two points can reduce the waste of boundary voxels, and can spare more space for further utilization.Step3: Remove the spherical region. Find whether the sphere of D=18mm can be packed into the remaining region, go on packing the sphere until the remaining region cannot be contain the sphere of this diameter. Repeat step2 using a smaller sized sphere until the remaining region cannot contain any sphere of all sizes.Step4: Find whether the total volume of all sized spheres meets the given requirements, and the sphere number is reasonable.Step5: Store the location of centers of all the spheres and their diameters, thus form an initial solution, which is a feasible Gamma Knife treatment plan.The solution may not be the best plan. A further optimal algorithm will be depicted in the next section.A two-dimensional example of above procedure:As a two-dimensional situation, the close region can be seen as a close curves. Spheres can be reduced to given sized circles. According to above algorithm, we select a close plane region arbitrarily and make a trial of above procedure; the procedure is shown in Figure 7.Figure 7: the procedure of circle packing using above algorithmProcedure details of Figure 7:1. Generate the 0-1 array of given region.2. Find the tracks of the largest circle center (Figure 8), if the track does not exist, that meansthe profile cannot totally contain the circle, we may choose a smaller sized circle, and find the track again.Figure 8: track of circle center3. Then, we may find two points on the track, the distance between which is the largest.Record their locations. Select one location, place the circle, and remove the circle area from initial region; repeat the above procedure, find the largest circle or circles (given four sizes) that can be placed in the remaining region until none of the four sized circles can be placed in the remaining region; record the result.4. Compute the area of all circles in each situation; if the computed area is more enough,such as 90% of origin area and the shot number is less, the better the result is.Optimization by Simulated AnnealingNotice that simulated annealing method can be applied to optimize problems. The initial solution obtained by Shelling Procedure must be further adjusted. The objective function is using as few spheres as possible to occupy maximum space (more than 90% of original space). Denote the objective function as:∑=−=N i i n r V F 1334π where n F is the volume of the residual fractions. V is the target volume. i r is the radius of the sphere i, N is the number of spheres. The solution space is a 4-dimensional space consisting of three continuous and one discrete variable parameters. Denote one solution as:(){}N i r z y x S i i i i ,,2,1,,,,"==where ()i i i z y x ,, is the coordinate of the center of the sphere i. The solution space is so large that they cannot be explored exhaustively. The adjustment algorithm is as follows:Step 1: Calculate the objective function value of the initial solution. Then we obtain 0F andS .Step 2: Obtain the shift matrix W . It must be have the same degree of S .W =(){}N i k q p h i i i i ,,2,1,,,,"=Where i i p h ,and i q can be obtained by random way, and all of these data obey the standard normal distribution from –ε to ε. When the temperature t m is high, we should set ε larger. Where i k is the special data which will be described as follows:Step 3: Obtain new potential solution by W S +.S ′=(){}N i q z p y h x i i i i i i ,,2,1,,,"=+++Where i k +i r must be one element of the set ( 0, 4,8,14,18). And i k +i r can only be near i r . For example, if i r = 18, then i k +i r can select 14 or 18. If i r =4, then i k +i r can be 0, 4 or 8. Step 4: Test if S ′ is the proper solution, in other words, if S ′ obeys the restriction: (1) the distance between each two center of the spheres is lager than or equal to the sum of the radius of these two spheres. (2) the distance from the center of any sphere to the boundary of the target volume is larger than or equal to the sphere’s radius. If S ′ meets the restriction, continue steps 5, or go to step 2.Step 5: Calculate the value of the new solution, then we obtain 0F ′. If 0F ′ is smaller than 0F , in other words, S ′ is better than S , so let S = S ′, 1F =0F ′. Then we have the new solution, go to step 6.Use Shelling Procedure to estimate if we can add one or more sphere to fill lacuna. Note that this way may lead to adding the dimension of the solution matrix. If 0F ′ is larger than 0F , then if()1,0exp 00random t F F m >⎟⎟⎠⎞⎜⎜⎝⎛−′ S ′ is better than S , go to step 6. Otherwise go to step 2.Step 6: When new better solution is found, if ζ>−+1n n F F (where ζ is a small number we preassigned.) go to step 1. We should reduce the temperature if ζ<−+1n n F F (which means we have found the best solution at this temperature), but how to realize it? We have many ways. In our program, we select Lundy and Mess’ method:mm m t t t ×+=+β11Where f ft t M t t ××−=00β0t and f t are two critical value. 0t is the lowest temperature and f t is the highest. Both of them have been designated at advance. M is the largest number of adjustable times. In paper [13], the iterate method is one time per each temperature, so the total of iterative times is M .Step 7: Let 1+=m m t t , set ε a lower value, go to step 1. One example on 2-D is as follows (Figure 9):(a) initial result (b) final resultFigure 9: comparison of initial result and the final resultFigure 9 (b) is adjusted from (a). The area of all circles in initial result to total area of origin shape is less than 84% and the shots number is 10, while the area of all circles in final result to total area of origin shape is more than 90% and the shot number is 9. So we get a satisfied result.Evaluation of Our ModelsModel A was established on the basis of dose distribution. We find that many optimization of Gamma Knife Radiosurgery are also based on this distribution. So Model A is available in theory. But in fact the model is a MIP problem. It generates a large enormous data, and is difficult to calculate within time variable. But in some case MIP is flexible and practical on account of that MIP could find global optimization.Shelling Procedure and Simulated Annealing cooperating with each other made our model and method preferable. Shelling Procedure has an advantage on the initial solution: we search for the possible area of spheres’ center. The area is very important to our work. All of the solutions are to be contained in this area. Based on this fact, we find better solutions quickly, but we cannot use this method to optimize the result which is produced by itself. Simulated Annealing needn’t deal with the difference of every target volume. We only need a few data. Computer can treat all other business, and the efficiency can meet the requirements. And itTeam # 167 Page 11 of 11can improve the solution undoubted. Furthermore, we can watch the total movement course, for example, our program can realize it like a cartoon movie. The process is similar to the genuine steel annealing. Like all Simulated Annealing, how to define solution space and objective function is very hard. It is inadvisable to use Simulated Annealing to find the initial solution, the cooperation of shelling procedure and the Simulated Annealing is a efficient way to solve this problem.References1.Taeil YiI, A Tree for a Brain Tumor. Florida Section Newsletter. February 2001,V olume22, Issue 2.2. A. Sutou and Y. Dai, A Study of the Global Optimaization Approach to Spherical PackingProblems, Dept. of Math. and Comp. Sciences Research Reports B-361, Tokyo Institute of Technology, May 2000, submitted, revised July 2001.3. A. Sutou and Y. Dai, Global Optimization Approach to Unequal Sphere PackingProblems in 3D, Journal of Optimization and Applications.4. A.Rosenfeld and A.C.Kak. Digital Picture Processing, V.2, 2nd Edition, Academic Press,1982.5.Kenneth R. Castleman,Digital Image Processing, Prentice Hall 2000.6.M. C. Ferris, MATLAB and GAMS interfacing optimization and visualization soft ware,Technical Report Mathematical Programming Technical Report98-19, Computer Sciences Department, University of Wisconsin, Madison, Wisconsin, 1998.7.M. C. Ferris, J. Lim, and D. M. Shepard, An optimization approach for Radiosurgerytreatment planning, SIAM Journal On Optimization, Forthcoming,(2002).8.M. C. Ferris, J. Lim, and D. M. Shepard, Radiosurgery optimization via nonlinearprogramming, Annals of Operations Research, Forthcoming, (2002).9.Michael C. Ferris, Jinho Lim, and David Shepard, Radiosurgery Optimization viaNonlinear Programming, Annals of Operations Research, vol 119, pp 247-260, 2003. 10.Jinho Lim, Optimization in radiation treatment planning, PhD Thesis, December 2002,University of Wisconsin - Madison.D. P. Bertsekas, Network Optimization: Continuous and Discrete Models, Athena Scientic, 1998.11.L. Luo, H. Shu, W. Yu, Y. Yan, X. Bao, and Y. Fu, Optimizing computerized treatmentplanning for the Gamma Knife by source culling, International Journal of Radiation Oncology, Biology and Physics, 45 (1999), 1339-1346.12.H. Z. Shu, Y. L. Yan, X. D. Bao, Y. Fu, and L. M. Luo, Treatment planning optimizationby quasi-newton and simulated annealing methods for gamma unit treatment system, Physics, Medicine and Biology, 43 (1998), 2795-2805.i K K, Chan J W M. Developing a simulated annealing algorithm for the cutting stockproblem, Computers Ind. Engng, 1996.14.J. Wang, Packing of unequal spheres and automated radiosurgical treatment planning,Journal of Combinatorial Optimization, 3 (1999), 453-463.。