Etienne Mabille

Marc Boyer

Loic Féjoz

Stephan Merz

Loïc Fejoz

Certifying Network Calculus in a Proof Assistant

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Result Certification for Network Calculus

In critical embedded systems such as those used in avionics, ensuring the quality of the software (using tests, qualification, and certification) represents a significant fraction of the overall cost of development. Several complementary methods exist to ensure the quality of systems, including rigorous development processes, testing of components and systems, and the use of formal methods. The latter have long been considered unable to scale up to industrial applications, and the cost of full-fledged proof remains prohibitively high in most cases. Nevertheless, advances in formal methods research have resulted in techniques that have proved practical and effective for specific problems. The Network Calculus [START_REF] Boudec | Network Calculus[END_REF] and supporting tool sets have found widespread use for determining bounds on message delays and for dimensioning buffers in embedded networks, such as in the design and certification of the Airbus A380 AFDX backbone [START_REF] Aeec | Arinc[END_REF][START_REF] Frances | Using network calculus to optimize AFDX network[END_REF][START_REF] Grieu | Analyse et évaluation de techniques de commutation Ethernet pour l'interconnexion des systèmes avioniques[END_REF]. However, results produced by existing tools based on the Network Calculus have to be trusted: although the underlying theory is generally well understood, implementation errors may result in faulty network designs, with unpredictable consequences.

For application domains subject to strict regulatory requirements, it appears mandatory to ensure the correctness of the results computed by the tools supporting network designs. In this contribution, we suggest a technique for certifying these results based on a "proof by instance" approach. In a nutshell (cf. Figure 1), the tool outputs a trace of the calculations it performs, as well as their results. The validity of the trace (both w.r.t. the applicability of the computation steps and the numerical correctness of the result) can be established offline by a trusted checker. For tools used at design time, we argue that this approach has several advantages over proving the calculator correct:

• Instrumenting the calculator for generating a trace is much easier (and hence less expensive) than attempting a full-fledged correctness proof, and checking the correctness of a computation is essentially trivial.

• The calculator is treated as a black box: it can be implemented using any software development process, programming language, and hardware by a tool provider separate from the checker. It can be updated without having to be requalified, as long as it still produces certifiable computation traces.

• Proof by instance is a good match for industrial processes based on testing. Nevertheless, a design that has been checked is guaranteed to be correct for any inputs matching the hypotheses of the model.

Checking the trace of a Network Calculus tool is similar to checking a mathematical proof: applying a given rule requires establishing its hypotheses, beyond pure calculation. We suggest to implement the checker by taking advantage of the trusted kernel of a proof assistant, specifically Isabelle/HOL [START_REF] Nipkow | HOL. A Proof Assistant for Higher-Order Logic[END_REF]. A companion paper describing this work from the perspective of interactive theorem proving appears in [START_REF] Mabille | Towards certifying network calculus[END_REF].

Network Calculus Theory

Books Article

Network Topology (flows and servers)

Network calculus

Mathematical background

Network calculus [START_REF] Boudec | Network Calculus[END_REF] is a theory for computing upper bounds in networks. Its mathematical background is a theory of the set of functions

F = { f : R ≥0 → R ≥0 ∪ {∞} | x ≤ y =⇒ f (x ) ≤ f (y) } (1)
that form a dioid under the operations ⊓ and + defined as pointwise minimum and addition. For practical applications, four families of functions are particularly interesting, which are defined as

δ d (t) = 0 if t ≤ d , ∞ otherwise, β R,T (t) = 0 if t ≤ T and R(t -T ) otherwise, and γ r ,b (t) = 0 if t ≤ 0 and rt + b otherwise.
Operations of interest on F include convolution ⊕, deconvolution ⊘, and the sub-additive closure f ⊕ .

(f ⊕ g)(t) = inf 0≤u≤t (f (t -u) + g(u)) (2)
(f ⊘ g)(t) = sup 0≤u (f (t + u) -g(u)) (3)
f ⊕ = δ 0 ⊓ f ⊓ (f ⊕ f ) ⊓ (f ⊕ f ⊕ f ) ⊓ • • • (4)
An important operator related to performance evaluation is the horizontal deviation h(•, •), defined by:

h(f , g) = sup s≥0 (inf { τ ≥ 0 | f (s) ≤ g(s + τ ) }) (5)
The horizontal deviation and some common curves are illustrated on Figure 2. These operators have important mathematical properties: for example, convolution is associative and commutative. It also "preserves" inequalities:

f ≤ g =⇒ f ⊕ h ≤ g ⊕ h.

Flow and server modelling

A flow is represented by its cumulative function R ∈ F, where R(t) is the total number of bits sent by this flow up to time t. R(t + s) -R(t) ≤ α(s), meaning that, from any instant t, the flow R will produce at most α(s) new bits of data in s time units. An equivalent condition, expressed in the

A flow R has function α ∈ F as arrival curve (denoted R α) if ∀t, s ≥ 0 : t γ r ,b β R,T h(γ r ,b , β R,T ) t γ r ,b δ D δ D ′ h(γ r ,b , δ D ) h(δ D , δ D ′ )
(⊓, +) dioid, is R ≤ R ⊕ α.
A server S is a relation between an input cumulative function R, and an output cumulative function

R ′ (denoted R S -→ R ′
) such that R ′ ≤ R (representing the intuition that the flow crosses the server, and that the output is produced after the input). Such a server has a service curve

β if R ′ ≥ R ⊕ β holds.
The maximal delay of a flow R crossing a server S is defined as the maximal horizontal deviation between the input and any possible output.

d (R, S ) = max R S -→R ′ h(R, R ′ ) (6)

Some network calculus results

The delay incurred by the flow R with arrival curve α at a server S with service curve β can be bounded by the horizontal deviation between curves α and β (cf. eq. 5, and also Fig. 2).

d (R, S ) ≤ h(α, β) (7)
In this case, network calculus provides several results:

• α ⊘ β is an arrival curve for the output flow R ′ ,

• δ d(α,β) is another service curve for the server S , and

• α ⊘ δ d(α,β) also is an arrival curve for the output flow R ′ (by simple application of the two previous items)

Having several results allows to choose one or another, depending on the context. For example, it can be shown that α ⊘ β ≤ α ⊘ δ d(α,β) , meaning that one result will produce smaller (i.e., better) numerical bounds. Conversely, the computation of α ⊘ δ x is usually much simpler than that of α ⊘ β.

The two following results express that if α is an arrival curve for some flow R, so is α ⊕ and also any

α ′ ≥ α. R α =⇒ R α ⊕ (8) R α, α ≤ α ′ =⇒ R α ′ (9)
If a flow crosses a sequence of two servers, S and S ′ (i.e. R S -→ R ′ S ′′ --→ R ′′ ), of respective service curves β and β ′ , then, the system made of the server sequence offers a service curve β ⊕ β ′ . This result is known as the pay burst only once phenomenon (PBOO).

The benefit of using a formal framework is that many results have a quite simple proof, based on operator properties. For example, equation ( 8) can be proved using associativity of convolution:

R ≤ R ⊕ α ≤ (R ⊕ α) ⊕ α = R ⊕ (α ⊕ α) ≤ . . . ≤ R ⊕ (α ⊕ • • • ⊕ α).
The pay burst only once result uses the same kind of proof:

R ′′ ≥ R ′ ⊕ β ′ ≥ (R ⊕ β) ⊕ β ′ = R ⊕ (β ⊕ β ′ ).

Using network calculus to compute bounds

This presentation gives a flavor of the essence of Network Calculus: a tool box of individual theorems that can be assembled to compute results. Consider a configuration with a flow R crossing two servers S 1 , S 2 in sequence: R

S1 -→ R ′ S2 -→ R ′′ .
Assume that R has arrival curve α and that each server S i offers a service of curve β i . Then, the delay of R in S 1 can be bounded by d = h(α, β 1 ). Among others, α ′ = α ⊘ δ d can be chosen as an arrival curve for R ′ . Its sub-additive closure (α ′ )

⊕ is also an arrival curve for R ′ (cf. eq. 8), but may be too expensive to compute. A simpler approximation is given by α ′ ⊓ δ 0 ≥ (α ′ ) ⊕ (using eq. 9), and the delay of R ′ in S 2 can be bounded by h(α ′ ⊓ δ 0 , β 2 ). The end-to-end delay can also be bounded by the sum of bounds on local delays, i.e. h(α,

β 1 ) + h(α ′ ⊓ δ 0 , β 2 ).
Other results can be used to compute bounds: for example, the PBOO computes h(α, β 1 ⊕β 2 ) as another bound on the global delay. The implementor of an analysis algorithm will choose among the available bounds based on a compromise between the accuracy of the result, the complexity of implementation, and the required computing time.

The Proof Assistant Isabelle/HOL

Our objective in this work is to represent the parts of Network Calculus that underly the computation of bounds in a trusted proof assistant. Isabelle [13] is a software framework in which the syntax and proof rules of logics can be formally defined, and that provides generic tools for machine-supported reasoning. In the tradition of the LCF family of proof assistants [START_REF] Gordon | Edinburgh LCF: A Mechanised Logic of Computation[END_REF], it is based on a small kernel written in the ML programming language that in particular defines a type thm representing theorems. Starting from user-defined axioms and proof rules specific to the particular logic, new theorems can only be created by applying rules to already existing theorems. In this way, the soundness of derived theorems only depends on that of the user-provided axioms and proof rules (which can easily be inspected) and the correct implementation of the kernel that implements rule application. The Isabelle framework provides many semi-automatic tools for proof search, such as a tableau-based theorem prover, a rewriting engine or decision procedures for fragments of arithmetic, and these tools can be instantiated for different logics. Moreover, users can program new proof methods (or tactics) from scratch or reusing the generic tools. All these methods provide justifications that are checked by the Isabelle kernel, and only then the result is accepted as a theorem. In this sense, result certification is extensively used within the Isabelle framework itself.

Isabelle/HOL [START_REF] Nipkow | HOL. A Proof Assistant for Higher-Order Logic[END_REF] is the encoding in Isabelle of higher-order logic, an expressive language for formalizing mathematical theories. HOL is a typed logic whose type system is close to that of standard functional programming languages. It provides base types such as bool , nat or real , function and product types such as real ⇒ bool or nat × nat, type constructors including lists or sets such as (real × bool ) set, and type variables written ′ a, ′ b etc. The function type constructor associates to the right, so real ⇒ real ⇒ bool is the type of functions that expect two real-valued arguments and return a boolean value. Function application is denoted by juxtaposition, and functions are first-class objects in HOL: if f is a function of the above type and r is an expression of type real , then f r denotes a unary function of type real ⇒ bool . For any expression, Isabelle will infer its most general type, and type annotations are therefore usually unnecessary. They may nevertheless be used for purposes of documentation or for restricting an expression to a less general type, as in 0 :: nat.

New types can be introduced as algebraic data types based on constructors, as in functional programming. For example, a data type of trees could be defined with constructors Leaf and Node that represent a leaf and an inner node of a tree. Isabelle/HOL also provides the typedef construct for introducing a type based on a characteristic predicate. For example, typedef even = {n :: nat . ∃m. n = m + m} introduces a type isomorphic to the set of even natural numbers. Since types must be non-empty in HOL for logical consistency, such definitions come with a proof obligation that requires the user to prove the existence of some value satisfying the predicate.

Developments in Isabelle are organized in theory files that contain definitions of types and operators as well as statements and proofs of lemmas and theorems. For example, definition is-upper-bound :: real set ⇒ real ⇒ bool where is-upper-bound S r ≡ ∀x ∈ S . x ≤ r introduces a predicate that indicates if a real number is an upper bound of a set of real numbers. Isabelle/HOL also supports operators defined inductively, as in the following definition of binomial coefficients

fun binom n 0 = 1 | binom 0 (Suc k ) = 0 | binom (Suc n) (Suc k ) = binom n k + binom n (Suc k )
where Suc denotes the successor function over natural numbers. Since HOL is a logic of total functions, inductive definitions must be shown to be terminating. In many cases, such as for the above function binom, this termination proof can be found automatically by built-in heuristics.

Lemmas and theorems assert facts, typically involving previously defined operators. For example, we could state the following lemma about upper bounds: lemma upper-bound-union:

assumes is-upper-bound S r and is-upper-bound T r shows is-upper-bound (S ∪ T ) r

All variables that appear in the statement of a theorem are implicitly universally quantified. The user must provide a proof that can be checked using Isabelle's proof methods. For a simple lemma such as the one above, one would simply expand the definition of the operator is-upper-bound, then invoke Isabelle's default automatic prover. In more complex cases, proofs could proceed by induction, possibly based on a generalization of the statement of the theorem. Finally, Isabelle/HOL also provides facilities for generating executable code in functional programming languages (including ML, Haskell, and Scala) from operator definitions. Programs generated in this way enjoy a high degree of trustworthiness since they are based on the same definitions that underly formal proofs about these operators. Although we do not currently rely on code generation for our result certifier, we may consider doing so in the future, since a stand-alone program can execute much more efficiently than when all computations are checked by the Isabelle kernel, as is the case for our current prototype.

Encoding Network Calculus in Isabelle

Mathematical Bases

In order to develop a result certifier within Isabelle, we need to formalize the theory underlying Network Calculus, to the extent that it underlies the algorithms whose results we intend to certify. As a side benefit of this formalization, we obtain a rigorous development of Network Calculus, including all possible corner cases that may be overlooked in pencil-and-paper proofs. We now give an outline of our formal development.

The fundamental notion in Network Calculus is that of a flow, modeled as a non-decreasing function f : R ≥0 → R ≥0 ∪ {∞}. We define a suitable type in Isabelle as

typedef ndf = {f :: ereal ⇒ ereal . (∀r ≤ 0. f r = 0) ∧ mono f }
where ereal is a pre-defined type corresponding to R ∪ {∞}, and mono f holds if f is a weakly monotonic function. Note that we extended the domain of f to arbitrary real numbers, fixing the result to be zero over negative reals, as this turned out to simplify the subsequent definitions. Over this function type, we define operations such as +, * , and ≤ by pointwise extension over the arguments, and establish basic algebraic properties: for example, the resulting structure forms an ordered commutative monoid with 0 and 1.

We also introduce specific classes of functions such as δ d (t), β R,T (t) and γ r ,b (t) mentioned in Section 2.1 and characteristic operations on flows. For example, the operation of convolution introduced in equation ( 2) is represented in Isabelle as

definition convol :: ndf ⇒ ndf ⇒ ndf (infix ⊕) where f ⊕ g ≡ Abs-ndf (λt. if t < 0 then 0 else Inf { f • (t -u) + g • u | u. 0 ≤ u ∧ u ≤ t ∧ u = ∞ })
where Abs-ndf denotes the injection from functions over type ereal to flows of type ndf , and where • denotes the evaluation of a flow at an argument of type ereal .

In our Isabelle theories we prove characteristic properties of flows and operations. For example, a flow is called sub-additive if its value at argument x + y is at most the sum of its values at x and y, formally and we prove that the convolution of two sub-additive flows is again sub-additive:

definition is-sub-additive :: ndf ⇒ bool where is-sub-additive f ≡ (∀x y. 0 ≤ x ∧ 0 ≤ y -→ f • (x + y) ≤ f • x + f • y) Producer ✲ 1 MBit/s Switch1 θ 1 = 1µs ✲ 10 MBit/s Switch2 θ 2 = 20µs ✲ 5 MBit/s Consumer in γ r ,b S 1 β R 1 ,θ 1 mid S 2 β R 2 ,θ 2 out
lemma convol-sub-add-stable: assumes is-sub-additive f and is-sub-additive g shows is-sub-additive (f ⊕ g) Similar definitions and corresponding properties are formalized for the other fundamental constructs of Network Calculus, mirroring the informal presentation in Section 2.1.

Servers and service curves

Network Calculus represents a simple server as a left-total relation between (input and output) flows such that the output flow is not larger than the input flow:

typedef server = {s :: (ndf × ndf ) set. (∀in. ∃out. (in, out) ∈ s) ∧ (∀(in, out) ∈ s. out ≤ in)}
and we define what it means for a flow to be constrained by an arrival curve α and for a server to provide minimum service β:

R α ≡ R ≤ R ⊕ α S β ≡ ∀(in, out) ∈ S : in ⊕ β ≤ out.
The delay between two flows f and g at argument t corresponds to the least τ such that the flow g evaluated at time t + τ exceeds flow f at time t:

definition delay :: ndf ⇒ ndf ⇒ eral ⇒ ereal where delay f g t ≡ Inf { τ. 0 ≤ τ ∧ f • t ≤ g • (t + τ ) }
and the maximal delay of a flow R crossing server S , noted d (R, S ) in equation ( 6), is the maximal delay between the input R and any possible output of the server. definition worst-delay-server :: ndf ⇒ server ⇒ ereal where worst-delay-server R S

≡ Sup { delay R R ′ t | t R ′ . t = ∞ ∧ (R, R ′ ) ∈ Rep-server S }
Again, we prove results relating these constraints to bounds on delays and backlogs. For example, the following theorem, corresponding to equation ( 7) of Section 2.3, provides a bound on the delay of a simple server

theorem d-h-bound :
assumes in α and S β shows worst-delay-server in S ≤ h-dev α β

where the horizontal deviation h-dev is the name of the function defined in equation ( 5). Moreover, we define constructs such as composing servers in sequence or packetization. Finally, these concepts are extended to multiple-input multiple-output servers that take vectors of flows as input and output.

Certifying a Simple Network Computation

In order to illustrate the use of our theories on a simple example, let us consider the producer-consumer setup shown in Fig. 3. The producer sends at most one frame every T = 20 ms. We assume that the payload is of maximal size 980 bytes. Assuming an overhead of 20 bytes per frame, the maximum frame size is MFS = 8000 bits. This flow is sent to a consumer via two switches with switching delays θ 1 = 1 µs and θ 2 = 20 µs. The physical links between the producer, the switches, and the consumer have bandwidths of 1, 10 and 5 MBit/s, respectively.

The NC model appears in the lower part of Fig. 3. Flow in is constrained by the arrival curve α in = γ r ,b where b equals MFS and r = MFS T = 8000 20×10 3 = 2 5 . The service curves are given by the function β Ri ,θi where the bandwidths are R 1 = 10 bit/µs and R 2 = 5 bit/µs, and the delays are θ 1 and θ 2 .

We are interested in the maximal delays that frames may incur. Using theorem d-h-bound , the delay at server 1 is bounded by h(α in , β 10,1 ), which evaluates to 801 µs. As explained in Section 2, the arrival curve of flow mid can be computed as

α mid = (α in ⊘ δ 801 ) ⊓ δ 0 = (γ 2 5 ,8000 ⊘ δ 801 ) ⊓ δ 0 = γ 2 5 , 41602 5
.

Continuing for the second server, its delay is at most h(α mid , β 5,20 ) = 42102 25 µs. Consequently, the overall delay incurred by frames equals 801 µs + 42102 25 µs = 62127 25 µs.

These computations are performed by a Java prototype and certified using Isabelle.

The trace shown in the appendix A consists of output of the prototype interleaved with Isabelle lemmas whose assertions and proofs were automatically generated.

Industrial effort

Formal methods are often presented as a way to increase software quality. We believe that they can also contribute to decreasing the costs of software development in case the software has to meet high standards of qualification.

The avionics industry has a long experience of ensuring high quality software: the DO-178B (and the new DO-178C) standards guide the development process, aiming at increasing quality by imposing rigorous requirements on development, documentation, verification, and configuration. Enforcing these requirements has a significant cost: following a stringent process increases the cost of software production by a factor that is estimated between 4 and 10. We expect that proof by instance can significantly decrease the cost of producing a tool set for computing bounds on network delays, implementing network calculus.

Such a tool is mainly composed of two elements:

1. a (min, +) library, implementing the basic operations on functions (including min, +, ⊕, ⊘, • ⊕ , h(•, •), etc.) and, 2. a network calculus engine, applying network calculus theorems on subsystems and producing numerical results.

Each component can have different implementations, handling more complex and accurate functions and algorithms.

From the COINC project [START_REF]ARC COINC home[END_REF] we know that it takes about 1 man-year to develop a state of the art (min, +) library, implementing the very generic algorithms of [START_REF] Bouillard | An algorithmic toolbox for network calculus[END_REF]. Developing a full AFDX analyser based on network calculus represents an effort of about 3 man-years (including the (min, +) library), based on our experience on the RTaW-PEGASE tool [START_REF] Boyer | PEGASE, a robust and efficient tool for worst case network traversal time[END_REF].

To evaluate the feasibility of proof by instance for a network calculus tool, a prototype was developed: this prototype, presented in the previous section, implements a simple version of the (min, +) library and some elementary network calculus algorithms.

The development was carried out during a 9-month internship. This effort includes the development of the network calculus engine, but also the modelling of a network calculus library in Isabelle (cf. "Network calculus in Isabelle" box in Figure 1), and trace generation (cf. "Computation traces" box in Figure 1).

Since the (min, +) library and the network calculus algorithms used in the current version are quite basic, the accuracy of the tool is not very good: it compute bounds that are, globally, two times bigger than a state of the art tool, such as the RTaW-PEGASE tool [START_REF] Boyer | PEGASE, a robust and efficient tool for worst case network traversal time[END_REF]. Moreover, the Isabelle library is not fully proved: some results are presently considered as axioms, and we plan to fully prove them, further increasing our confidence in the correctness of the calculations.

Nevertheless, the low effort required for the current prototype confirms our intuition that proof by instance is a good approach for producing high-quality software at a reasonable cost. The quality of the result is ensured by the Isabelle kernel that implements conceptually simple operations and underlies many significant formal developments [START_REF] Klein | seL4: formal verification of an operating-system kernel[END_REF].

Conclusion

We argue that the computations of tool sets used in the applications of formal methods such as Network Calculus should be certified in order to increase the confidence in the correctness of the designs. The work reported here attempts to evaluate the feasibility of such certification, based on a formalization of Network Calculus in the interactive proof assistant Isabelle/HOL.

Developing a Network Calculus engine that is able to handle an AFDX configuration requires about three years of implementation. The effort for developing a qualified version of such an engine, using state-ofthe-art techniques (documentations, testing, peer-review, etc.) is higher by a factor of 4 or 10.

The proof-by-instance approach promises to reduce this effort while increasing the confidence in the results produced by the software. We have so far invested less than 1 development year for encoding the fundamental concepts of Network Calculus in Isabelle, and for instrumenting an existing tool to produce a full formal proof of the correctness of bounds for one single switch. We estimate that the overall effort for producing the proof for a realistic network should be between 2 and 3 years. In particular, it will be important to speed up the computations on real numbers inside the proof assistant.

In other words, we believe that result certification could lower the overhead for developing a trustworthy version of a Network Calculus tool to a factor of 2 or 3, while significantly improving its quality.

[13] Makarius Wenzel, Lawrence C.

A. Example Trace

1 ################################# 2 # Used algorithm: SFA FIFO 3 ################################# 4 # Time unit: microsecond 5 # Frame size unit: bit 6 ################################# 7 delta0 := delay(0) 8 ################################# 9 # Input flows 10 ################################# 11
# unique_flow entering at System_1>port-port1

12 lmax_unique_flow := 8000 lemma AC_sub_add_closure: "sub_add_closure (gamma (2/5) 8000) = gamma (2/5) 8000" using sub_add_closure_gamma by auto assumes "(R, R') ∈ Rep_server S" and "S beta 10 1" and AC_in: "R ≺ gamma (2/5) 8000" shows "R' ≺ gamma (2/5) (41602/5)" proofhave R'_cont: "continuous_on (range ereal) (Rep_ndf R')" sorry have "worst_delay_server R S ≤ ereal 801" using assms delay_server1 by auto with Out_AC_Delay have "R' ≺ gamma (2 / 5) 8000 °⊘°delta (ereal 801)" using assms gamma_real PInfty_neq_ereal(1) R'_cont by auto thus ?thesis using deconv_gamma_delta[of "ereal 801" "2/5" "8000"] by auto qed

13 14 ################################# 15 # EndSystem : System_1>port-port1 16 #################################
1 ################################# 2 # Server : theRouter>port-P3 3 #################################

Figure 1 :

 1 Figure 1: Proof by instance process

Figure 2 :

 2 Figure 2: Common curves and delay.

Figure 3 :

 3 Figure 3: A simple system and its Network Calculus representation.

#

  Packet arrival curve for unique_flow on System_1->theRouter 20 unique_flow_System1portport1_P := star(uaf([(0,0)]](0,8000)2/5(+Infinity,+Infinity)[)) 21 22 assert(unique_flow_System1portport1_P = uaf([(0,0)]](0,8000)2/5(+Infinity,+Infinity)[))

4 # Cumulated arrivals 5 cumA_theRouterportP3 := zero 6 # unique_flow on System_1->theRouter 7 cumA_theRouterportP3 := cumA_theRouterportP3 + unique_flow_System1portport1_P 8 assert 1 # Common service 2 S_theRouterportP3

 4567812 (cumA_theRouterportP3 = uaf([(0,0)]](0,8000)2/5(+Infinity,+Infinity)[))% lemma AC_add_zero: "(gamma (2/5) 8000) + (0::ndf) = gamma (2/5) 8000" by auto lemma delay_server1: assumes "(R, R') ∈ Rep_server S" and "S beta 10 1" and "R ≺ gamma (2/5) 8000" shows "worst_delay_server R S ≤ ereal(801)" prooffrom assms d_h_bound have "worst_delay_server R S ≤ h_dev (gamma (2/5) 8000) (beta 10 1)" by auto with hor_dev_beta_gamma show "worst_delay_server R S ≤ ereal (801)" by auto qed

6 # 1 # Common service 2 S_theSecondRouterportP2 2 # WCTT 3 # 4 # 5 #

 6122345 ############################## Server : theSecondRouter>port-P2 7 ################################# 8 = uaf([(0,0)]](0,41602/5)2/5(+Infinity,+Infinity)[))% lemma AC'_add_zero: "gamma (2/5) (41602/5) + (0::ndf) = gamma (2/5) (41602/5)" by auto assumes "(R', R'') ∈ Rep_server S2" and "S2 beta 5 20" and "R' ≺ gamma (2/5) (41602/5)" shows "worst_delay_server R' S2 ≤ ereal(42102/25)" prooffrom assms d_h_bound have "worst_delay_server R' S2 ≤ h_dev (gamma (2/5) (41602/5)) (beta 5 20)" by auto with hor_dev_beta_gamma show "worst_delay_server R' S2 ≤ ereal (42102/25)" by auto qed1 ################################# ################################ Boundfor unique_flow received by System_2 internal delay in sending server System_1>port-port1

  Paulson, and Tobias Nipkow. The Isabelle framework. In Otmane Aït Mohamed, César A. Muñoz, and Sofiène Tahar, editors, 21st Intl. Conf. Theorem Proving in Higher-Order Logics (TPHOLs 2008), volume 5170 of LNCS, pages 33-38, Montréal, Canada, 2008. Springer.

© 2013 by E. Mabille, M. Boyer, L. Fejoz and S. Merz. Published by the EUCASS association with permission.

# unique_flow on theRouter->theSecondRouter


Van Luong Le

Fabien Lauer
email: fabien.lauer@loria.fr

Gérard Bloch
email: gerard.bloch@univ-lorraine.fr

Selective ℓ 1 minimization for sparse recovery

Keywords: Convex relaxation, sparsity, system identification, hybrid systems, switched systems

Motivated by recent approaches to switched linear system identification based on sparse optimization, the paper deals with the recovery of sparse solutions of underdetermined systems of linear equations. More precisely, we focus on the associated convex relaxation where the ℓ1-norm of the vector of variables is minimized and propose a new iteratively reweighted scheme in order to improve the conditions under which this relaxation provides the sparsest solution. We prove the convergence of the new scheme and derive sufficient conditions for the convergence towards the sparsest solution. Experiments show that the new scheme significantly improves upon the previous approaches for compressive sensing. Then, these results are applied to switched system identification.

I. INTRODUCTION

This paper considers a convex relaxation approach to recover sparse solutions of underdetermined systems of linear equations, with an application to hybrid dynamical system identification. Linear hybrid systems are dynamical systems switching between multiple linear subsystems. More precisely, we consider single-input single-output systems in switched autoregressive with external input form as

y i = θ ⊤ λi ϕ i + v i , (1)
where λ i ∈ {1, . . . , s} is the discrete state or mode with s the number of submodels, {θ j } s j=1 are the parameter vectors of the submodels, v i ∈ R is a noise term and ϕ i = [y i-1 , . . . , y i-na , u i-1 , . . . , u i-n b ] ⊤ ∈ R p is the regression vector with p = n a + n b , where n a and n b are the model orders. Then, the identification problem is, given a collection S = {(ϕ i , y i )} N i=1 , to estimate: (i) the number of submodels s, (ii) the parameter vectors {θ j } s j=1 , and (iii) the switching sequence {λ i } N i=1 . Related work. Many hybrid system identification approaches have been proposed over the last decade [START_REF] Garulli | A survey on switched and piecewise affine system identification[END_REF]. Here, we focus on methods based on convex optimization and which offer several guarantees, the first of which being that the estimates correspond to a global solution of the formulated optimization problem. In addition, when a convex relaxation of a nonconvex optimization problem is considered, theoretical guarantees of equivalence between the two formulations can be obtained. More particularly, we follow the approach of [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF], which iteratively estimates each parameter vector individually by maximizing the sparsity of the error vector. This sparse optimization problem is solved via its ℓ 1 -norm convex relaxation and the sparsity of the solution is improved by the iteratively reweighted ℓ 1 minimization scheme developed in the compressive sensing literature [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. Note that many results on ℓ 1 -norm relaxations were also developed in this field [START_REF] Davenport | Introduction to compressed sensing[END_REF], [START_REF] Fornasier | Compressive sensing[END_REF], [START_REF] Donoho | Compressed sensing[END_REF].

Contribution. This paper proposes in Sect. II a new iterative method based on ℓ 1 -norm minimization for the recovery of sparse solutions. As in the method of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] (recalled in Sect. II-A), we consider a weighted form of the ℓ 1 convex relaxation. But, instead of updating the weights in a soft manner, in Sect. II-B we explicitly set a weight to zero at each iteration. The proposed scheme offers three major advantages when compared with the one of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]: (i) it converges in a finite number of steps, (ii) theoretical guarantees of convergence towards the sparsest solution can be obtained, and (iii) experiments in Sect. IV-A show that it allows the sparsest solution to be recovered in a larger range of sparsity level. The advantages of this new sparsity enhancing scheme are used in Sect. III to improve the approach of [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF] for hybrid system identification, as witnessed by experiments in Sect. IV-B.

II. SPARSE RECOVERY

Consider an underdetermined system of linear equations, Ax = b, with a full row rank matrix A ∈ R m×n and a non-zero vector b ∈ R m , where m ≪ n. We are interested in sparse solutions of this system, i.e., solutions with few nonzero components, which are specifically obtained by solving

min x∈R n x 0 , s.t. Ax = b, (2)
where the ℓ 0 -pseudo norm is defined as

x 0 = |{i : x i = 0}|.
More precisely we concentrate on instances of (2) with a unique minimizer and assume that the following assumption holds in the rest of the paper.

Assumption 1. Problem (2) has a unique minimizer.

The following theorem shows that Assumption 1 holds in many cases.

Theorem 1 (Uniqueness via the spark [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF]). If a system of linear equations Ax = b has a solution x obeying

x 0 < spark(A) 2 ,
where spark(A) is the smallest number of linearly dependent columns of A, this solution is necessarily the sparsest possible.

Even when having a unique solution, (2) remains a nonconvex optimization problem which is intractable for large n due to its combinatorial search nature. Nonetheless, it has been the focus of many works over the last decade in various fields and particularly in the context of compressive sensing [START_REF] Davenport | Introduction to compressed sensing[END_REF], [START_REF] Fornasier | Compressive sensing[END_REF], [START_REF] Donoho | Compressed sensing[END_REF].

As discussed in [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF] and references therein, a common alternative for (2) is to consider the convex relaxation based on the ℓ 1 -norm. This leads to

min x∈R n x 1 , s.t. Ax = b, (3)
where

x 1 = n i=1 |x i | is the ℓ 1 -norm of x.
This problem is convex and can typically be solved efficiently [START_REF] Boyd | Convex Optimization[END_REF]. As shown in [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], for a column normalized matrix A, problems (2) and (3) are equivalent if, for a solution of Ax = b, the condition

x 0 < 1 2 1 + 1 µ(A) (4)
holds, where µ(A) is the mutual coherence of A [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF] which is defined, for any matrix

A = [A 1 , . . . , A n ] ∈ R m×n , by µ = µ (A) = max 1≤i,j≤n,i =j A ⊤ i A j A i 2 A j 2 . ( 5
)
In the case of an unnormalized matrix A, a similar equivalence is obtained by considering a weighted version of (3), i.e., min

x∈R n W A x 1 , s.t. Ax = b, where W A = diag ( A 1 2 , . . . , A n 2 ).
However, in many applications, the matrix A cannot be freely chosen and the sufficient condition (4) might be violated. Yet, the problem defined in (2) may have a unique solution x 0 , as stated by Theorem 1. Thus, since 1 + 1 µ(A) ≤ spark(A) [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], there is a range of problems with

1 2 1 + 1 µ(A) ≤ x 0 0 ≤ spark(A) 2 .
For these problems, recovering the sparsest solution is fore a well-defined problem, but not directly solvable through the ℓ 1 convex relaxation: a solution x 1 to (3) may have more nonzero elements than x 0 .

A. The classical iteratively reweighted approach

In order to improve the sparsity of the solutions of (3), a reweighted ℓ 1 minimization scheme is proposed in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. At each iteration l, the following problem is solved:

x (l) 1 ∈ arg min x∈R n W l x 1 , s.t. Ax = b, (6)
where

W l = diag(w (l) 1 , . . . , w (l)
n ) is a weighting diagonal matrix which penalizes differently the entries of x. At the first iteration, the weights are equal, i.e., W 1 = I n . Then, W l is updated with

w (l+1) i = 1 x (l) 1i + ǫ ,
where the w (l+1) i

are the weights at the (l + 1)th iteration, x (l) 1i

is the ith element of the solution of (6) at the lth iteration and ǫ > 0 is a parameter preventing a division by zero. Note that the choice of ǫ has an influence on the convergence of x (l)

1 . Important open issues, highlighted in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF], regarding this scheme are: (i) What are smart and robust rules for selecting the parameter ǫ? (ii) Under what conditions does the algorithm converge? Though preliminary results on the convergence are given in [START_REF] Xu | Breaking through the thresholds: an analysis for iterative reweighted ℓ 1 minimization via the Grassmann angle framework[END_REF], these rely on a condition involving both the matrix A and the solution to [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]. On the contrary, the following presents another reweighting mechanism with a convergence analysis relying only on conditions on A.

Note that other methods have been proposed to recover sparse solutions, including the greedy algorithm of [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF] and the one in [START_REF] Novara | Sparse identification of nonlinear functions and parametric set membership optimality analysis[END_REF] that mixes ℓ 1 minimization with a greedy approach.

B. Selective ℓ 1 minimization

In this section, we propose a new method for updating the weighting matrix W l in (6) in order to improve the sparsity of the solution. The new method, named Selective

ℓ 1 Minimization (Sℓ 1 M), is given in Algorithm 1, where x ∞ = max i |x i | is the ℓ ∞ -norm of x. Algorithm 1 Sℓ 1 M Require: A ∈ R m×n and b ∈ R m . -Initialize l = 0, W 1 = I n . repeat -Set l = l + 1. -Get any x (l)
1 in the solution set of ( 6). -Select the smallest index q (l) ∈ arg max i=1,...,n

w (l) i |x (l) 1i |. -Calculate W l+1 with w (l+1) i = w (l) i , if i = q (l) , 0, if i = q (l) . until W l x (l) 1 ∞ = 0 or W l+1 x (l) 1 ∞ = 0. return x * 1 = x (l) 1 .
Algorithm 1 relaxes the optimization of the nonzero variables by setting their weights w i to 0 in the cost function of [START_REF] Donoho | Compressed sensing[END_REF], thus putting more weight on the other variables that are pulled towards 0. When the stopping criterion is met, we have W l x (l) 1 0 ≤ 1. Hence, if it returns at iteration l < n, the algorithm yields a sparse solution.

The following provides an analysis of the proposed iterative scheme. The convergence in a finite number of steps is proved and a condition on the matrix A and the sparsity level guaranteeing the convergence towards the desired solution is derived.

1) Convergence in a finite number of steps:

The following theorem, based on Lemmas given in Appendix A, guarantees that the algorithm Sℓ 1 M converges in a finite number of steps.

Theorem 2. The solution x * 1 , returned by Algorithm 1, is found in at most m + 1 iterations and x * 1 0 ≤ m, where m is the number of rows in A.

Proof: If Algorithm 1 does not converge after m iterations, Lemma 3 implies that x (l) 1q (l) = 0, ∀l ≤ m. Then, according to Lemma 4, the columns A q (l) , l ∈ {1, . . . , m}, of A are linearly independent and b ∈ R m can be expressed as a linear combination of these columns. Therefore, the minimum value of the sum to (6) at iteration m + 1 satisfies the stopping criterion of the algorithm as

W (m+1) x (m+1) 1 ∞ = 0.
In addition, if Algorithm 1 converges at the lth iteration,

x (l+1) 1 0 ≤ l -1. Thus, x * 1 0 ≤ m. It is worth noting that if x * 1 0 < spark(A)

2

, according to Theorem 1, x * 1 is the sparsest solution of (2).

2) Convergence towards the sparsest solution:

In order for the iterative algorithm to converge to the solution x 0 of (2) under Assumption 1, we need to ensure that the variables removed from the weighted sum in the cost function of (6) correspond to nonzeros in x 0 . The following proposition shows under which conditions the choice of the index q (l) in Algorithm 1 corresponds to a nonzero element in the solution of (2), x 0 , for which we define the two sets

I 0 = {i : x 0i = 0}, I 1 = {i : x 0i = 0}. ( 7
) Proposition 1. If x 0 -x (l) 1 1 < µ+1 2µ W l x 0 ∞
, where µ is defined as in [START_REF] Fornasier | Compressive sensing[END_REF] with a column normalized matrix A, then q (l) ∈ I 1 for all q (l) ∈ arg max i w

(l) i x (l) 1i . Proof: Assume q (l) / ∈ I 1 , x 0q (l) = 0, then x 0 -x (l) 1 q (l) = |x (l) 1q (l) | = max i w (l) i x (l) 1i . Let j ∈ arg max i w (l) i |x 0i | and |x 0j | = W l x 0 ∞ > 2µ µ+1 x 0 -x (l) 1 1
. Then, since x

(l) 1j ≥ |x 0j |-x 0 -x (l) 1 j , Lemma 5 in Appendix applied to δ = x 0 -x (l) 1 leads to x (l) 1j > µ µ+1 x 0 -x (l) 1 1
. Thus

x (l) 1q (l) = x 0 -x (l) 1
q (l) ≥ x (l) 1j > µ µ + 1 x 0 -x (l) 1 1
.

(8) On the other hand, Lemma 5 implies that (x 0x

(l) 1 ) q (l) ≤ µ µ+1 x 0 x (l) 1 1
, which contradicts [START_REF] Boyd | Convex Optimization[END_REF]. Thus, the assumption q (l) / ∈ I 1 is wrong and we conclude that q (l) ∈ I 1 .

3) Influence of the weighting matrix on sparse recovery:

The following lemma shows that, with a good choice of W such that w i = 0, i ∈ I 1 , the solution to (6) is exactly x 0 . Lemma 1. Given a diagonal matrix W , with entries w i ≥ 0, if for all nonzero δ ∈ Ker(A) the following condition holds

i∈I1 w i |δ i | < i∈I0 w i |δ i |, ( 9
)
where I 0 and I 1 are defined by [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], then the solution x 0 to (2) under Assumption 1 uniquely solves problem (6), i.e.,

x 0 = arg min x∈R n W x 1 s.t. Ax = b. Proof: x 0 uniquely solves problem (6) if, for all nonzero δ ∈ Ker(A), W x 0 1 < W (x 0 + δ) 1 , or equivalently if i∈I1 w i (|x 0i | -|x 0i + δ i |) < i∈I0 w i |δ i | ,
which is implied by the condition (9) since

i∈I1 w i (|x 0i | -|x 0i + δ i |) ≤ i∈I1 w i |δ i |.
If the condition in Proposition 1 is satisfied for h iterations, the condition ( 9) is relaxed to

i∈I1 |δ i | - i∈{q (l) } h l=1 |δ i | < i∈I0 |δ i |. ( 10
)
When this condition is satisfied, the solution of ( 6) with W l = W h+1 is the sparsest vector x 0 thanks to Lemma 1. We see that the left-hand side (LHS) of ( 10) can decrease to zero after |I 1 | iterations. This also shows that the algorithm converges after at most |I 1 | + 1 iterations if the indexes q (l) are all wellchosen in I 1 .

4) Sparse recovery condition:

We now show, in Theorem 3, a condition on the matrix A and the sparsity level which can guarantee the convergence towards the desired solution. But this requires the following definition using the notation δ T for the subvector of δ containing its components of indexes in T .

Definition 1 (Definition 1 in [START_REF] Fornasier | Compressive sensing[END_REF]). A matrix A ∈ C m×n is said to satisfy the null space property (NSP) of order k with constant γ ∈ (0, 1) if

δ T 1 ≤ γ δ T c 1 ,
for all sets T ⊂ {1, . . . , n} with |T | ≤ k, T c = {1, . . . , n} \T and for all δ ∈ Ker(A).

The following theorem extends Theorem 1 in [START_REF] Fornasier | Compressive sensing[END_REF] to cases where the sparsity level x 0 0 equals k + h with h > 0 and provides a sparse recovery condition for the iterative Algorithm 1 (proof given in Appendix B).

Theorem 3. Given a matrix A that satisfies the NSP of order k with constant γ ∈ (0, 1 2 ), if x 0 is such that Ax 0 = b and

x 0 0 ≤ k + h with the integer h ∈ [1, k] satisfying γ < 1 -(4h -1)µ 1 + (4h + 1)µ , ( 11
)
where µ is defined as in [START_REF] Fornasier | Compressive sensing[END_REF], Algorithm 1 converges to x 0 in at most h + 1 iterations.

Note that the NSP of A and the value of γ can be difficult to determine directly, but can be related to the easier to handle restricted isometry property (see [START_REF] Fornasier | Compressive sensing[END_REF] for details).

III. HYBRID SYSTEM IDENTIFICATION

We now turn to the problem of hybrid system identification, i.e., of estimating the parameter vectors

{θ j } s j=1 in (1) from a data set S = {(ϕ i , y i )} N i=1 .
We follow the approach of [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF] in which we introduce the proposed reweighted scheme (Sℓ 1 M). This yields Algorithm 2 for the estimation of a single parameter vector.

Algorithm 2 requires a number of iterations, N s , in order to deal with noisy data (for which true sparsity cannot be obtained). However, with knowledge of the number of modes, s, we can set N s = s-1 s N , since in this case the largest fraction of points of a mode in the data set is at least N/s.

As in [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF], the identification procedure obtains the submodels one by one. After applying Algorithm 2 to estimate a parameter vector θj , the data points verifying the error condition, |y iϕ ⊤ i θj | ≤ δ where δ is a fixed threshold, are associated to this submodel and removed from the data set. Then, the next parameter vectors are iteratively estimated from reduced data sets until all data points are removed, at which point the estimated number of modes ŝ is given by the number of submodels obtained. Note that, if the noise is unbounded, -Select an index in the maximal absolute error set (break the tie arbitrarily if necessary):

q ∈ arg max i w i |e i (θ * )| .
-Set the qth entry on the diagonal of W to w q = 0. end while return θ * and W . e.g., Gaussian, the procedure should stop before that, i.e., when a small and predefined fraction of the data remains, to avoid creating irrelevant submodels for small groups of points corrupted by large noise terms.

IV. NUMERICAL EXPERIMENTS

A. Compressive sensing example

We consider a classical example of sparse signal recovery used in many works to show the efficiency of the proposed method. The goal is to recover a sparse signal x of length n with x 0 = k. The k nonzero positions are chosen randomly, and the nonzero values are randomly drawn according to a zero-mean unit-variance Gaussian distribution. The sensing matrix A ∈ R m×n is a Gaussian matrix, i.e., with entries following a zero-mean Gaussian distribution with variance 1/m. The sparsity level k is increased from 10 to 60 to see the capacity in signal recovery.

To compare with the classical reweighting method, we set the experiment as in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] with n = 256, m = 100. For each value of k, we run 500 trials to estimate the probability of perfect signal recovery (be successful if x 0 -x0 ∞ ≤ 10 -3 ). Figure 1 reports the successful recovery probability, Pr(recovery), for the unweighted ℓ 1 -norm minimization (Unweighted ℓ 1 ), the Orthogonal Matching Pursuit (OMP) [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], the reweighted scheme of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] (Reweighted), and the proposed one (Sℓ 1 M). We see that the requisite oversampling factor for perfect recovery [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF], min k m/k s.t. Pr(recovery)= 1, decreases from approximately 4 for Unweighted ℓ 1 or 3 for the method of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] with ǫ = 1 to 100/40 = 2.5 for our method. Moreover, in our method there is no hyperparameter to tune, whereas ǫ can influence the results for the classical scheme [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF].

From Theorem 1, we must have k < spark(A)/2 ≤ m/2 = 50 to guarantee the uniqueness of the solution to (2). This explains why all methods have a small successful recovery probability with the sparsity level k close to 50. Nonetheless, our method shows a successful recovery probability greater than 0.9 at k = 45.

These improvements are paid for by the computational cost of the proposed method, which requires a larger number of Fig. 1. Empirical probability of successful recovery versus the sparsity level k for the OMP method [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], the classical reweighting of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] with various ǫ and the proposed one (Sℓ 1 M). iterations compared with [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. This trade-off is consistent with the results of the OMP [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], which obtains a low probability of success but at a much lower computational cost.

B. Hybrid system identification example

We now consider the switched linear system with 3 modes and n a = n b = 2 given as an example in Sect. 4 of [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF] to test the identification procedure. Training sets of N = 600 points are generated with a uniformly distributed random sequence of λ i ∈ {1, 2, 3} and an additive Gaussian noise with σ v = 0.1.

We compare the Normalized Parametric Error,

N P E = s j=1 θj -θ j 2 2 / s j=1 θ j 2 2
, of the original method of [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF] (SO) using the reweighting of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] with the one of the same method based on Sℓ 1 M and Algorithm 2 (labeled Sℓ 1 M for short). More precisely, over 100 trials with different input, switching and noise sequences, we report the mean and standard deviation of the NPE. Since the methods estimate the parameter vectors one by one until the data set is empty, the number of modes cannot be fixed. If ŝ > s, the NPE is computed with the s best parameter vectors that yield the smallest NPE. The threshold δ used to assign data points to a submodel is varied in the range [2σ v , 7σ v ]. Due to the unbounded Gaussian noise, both methods are stopped with 5% of the data left unassigned to a mode.

Figure 2 shows that the proposed method (with N s set as suggested in Sect. III) yields a model with a smaller error and a better estimate of the number of modes than the SO method. Moreover, this is true over a large range of values for δ.

V. CONCLUSIONS

The paper proposed a new iterative algorithm to improve the sparsity of the solution of an ℓ 1 -norm relaxation. Compared with the state-of-the-art scheme of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF], the proposed algorithm benefits from the absence of hyperparameters and a finite convergence in a number of steps at most equal to the number of linear equations. In addition, a sparse recovery condition, guaranteeing the convergence towards the sparsest solution, was proved and experiments showed that the new scheme can recover the sparsest solution in more difficult cases. Finally, we presented an application to hybrid system identification, where the increased sparse recovery capacity of the method translates into more accurate parameter estimates.

APPENDIX

Throughout the appendix, for a vector x ∈ R n and an index set T ⊂ {1, . . . , n}, we denote x T the subvector of x containing its components of indexes in T and x T 1 = i∈T |x i |. Proof: We know that x 1T 1 ≤ x1T 1 or equivalently

A. Useful lemmas and definitions

Lemma 2. Given a solution x 1 of min x∈R n x T 1 , s.t. Ax = b, if x 1i = 0 for some i ∈ T ,
x 1 T 1 + |x 1i | ≤ x1 T 1 + |x 1i | while x 1 T 1 ≥ x1 T 1 . Therefore |x 1i | ≥ |x 1i | > 0. Lemma 3. For b = 0, after l iterations, if W l x (l) 1 ∞ = 0, Algorithm 1 yields x (l) 1 such that x (l)
1q (i) = 0, ∀i = {1, . . . , l}. Proof: Lemma 3 is a consequence of Lemma 2 and the fact that

W 1 x (1) 1 ∞ = x (1) 1q (1) = 0 if b = 0. Lemma 4. In Algorithm 1, if, ∀l ≤ m, x (l)
1q (l) = 0, then the columns A q (i) , i = 1, . . . , l, of the full row rank matrix A are linearly independent.

Proof: To prove that, we show that at the jth iteration, ∀j ≥ 2, if x (l) 1q (l) = 0, ∀l ≤ j, then A q (j) is linearly independent of A q (l) j-1 l=1 . Assume that this is not true, i.e., A q (j) = (j-1)

l=1 β l A q (l) . Then, b = Ax (j) 1 = j-1 l=1 x (j) 1q (l) + x (j) 1q (j) β l A q (l) + i / ∈{q (l) } j l=1 x (j)
1i Ai, and we can get another solution x * (j) of b = Ax whose elements are given as the coefficients of the columns of A in the above. In particular, x * (j) q (j) = 0 and W j x *(j)

1 = i / ∈{q (l) } j l=1 x (j) 1i < |x (j) 1q (j) | + i / ∈{q (l) } j l=1 x (j) 1i = W j x (j) 1 1
. But this contradicts the definition of x (j)

1 by ( 6) and we conclude that A q (j) is linearly independent of A q (l) j-1 l=1 . Lemma 5. For any δ ∈ R n and A ∈ R m×n such that Aδ = 0 and all columns of A are unit vectors, the bound

|δ i | ≤ µ µ + 1 δ 1 ,
where µ is the mutual coherence of A, holds.

Proof: If Aδ = 0, then A ⊤ A -I n δ = -δ. Thus,

|δ i | ≤ n j=1 |δ j | A ⊤ A -I n i,j ≤ n j=1 µ |δ j | -µ |δ i | .
Rearranging the terms yields the sought statement.

Definition 2 (taken from [START_REF] Fornasier | Compressive sensing[END_REF]). Let Σ k = {x ∈ R n : x 0 ≤ k}. The best k-term approximation error in terms of ℓ p -norm of a vector x is defined as

σ k (x) p = inf z∈Σ k x -z p .

B. Proof of Theorem 3

We now give a condensed proof of Theorem 3, while full details can be found in chapter 4 of [START_REF] Le | Hybrid dynamical system identification: geometry, sparsity and nonlinearities[END_REF]. The proof is decomposed in two main steps: (i) showing that q (l) ∈ I 1 , ∀l ≤ h, (ii) showing that, after Step (i), the algorithm converges to the unique solution x 0 at iteration h + 1.

Step (i). First, we will prove by induction that ∀l ≤ h,

x 0 -x (l) 1 1 < µ + 1 2µ W l x 0 ∞ , (12)
from which we will use Proposition 1 to conclude that q (l) ∈ I 1 , ∀l ≤ h. To prove [START_REF] Novara | Sparse identification of nonlinear functions and parametric set membership optimality analysis[END_REF], we follow a path similar to that of the proof of Theorem 1 in [START_REF] Fornasier | Compressive sensing[END_REF]. We define the sets T (j) q = q (i) j-1 i=1 , T (j) s = {1, . . . , n} \T (j) q , T (j) as the set of index of k entries of x 0 with largest magnitude such that T (j) ∩ T (j) q = ∅ and T

(j) r = {1, . . . , n} \ T (j) ∪ T (j) q . Let δ (j) = x 0 -x (j) 1 . We have δ (l) 1 = δ (l) T (l) 1 + δ (l) T (l) q 1 + δ (l) T (l) r 1
, where the first term on the right-hand side (RHS) can be bounded as follows. By the fact that A satisfies the NSP (Definition 1) and that δ (l) ∈ Ker(A), the definitions of T (l) and T

(l) r imply δ (l) T (l) 1 ≤ γ δ (l) T (l) q 1 + δ (l) T (l) r 1 . (13)
Therefore, δ (l) 1 is bounded by

δ (l) 1 ≤ (1 + γ) δ (l) T (l) q 1 + δ (l) T (l) r 1 . (14)
Now, we bound the second term of the RHS of (14) as follows.

For l ≥ 2, we have

x 0T (l) 1 + x 0T (l) q 1 + x 0T (l) r 1 = x0 1 ≥ x (1) 1 1 ≥ x (2) 1T (2) s 1 + x (1) 1q (1) ≥ x (3)
1T

(3) s 1 + x (1)
1q

1) + x (2) 1q (2) ≥ . . . ≥ x (l) 1T (l) s 1 + l-1 i=1 x (i) 1q (i) = x (l) 1T (l) 1 + x (l) 1T (l) r 1 + l-1 i=1 x (i) 1q (i) . (
Then, by using the triangle inequalities, x

l) 1T (l) 1 ≥ x 0T (l) 1 -δ (l) T (l) 1 and x (l) 1T (l) r 1 ≥ δ (l) T (l) r 1 -x 0T (l) (
T (l) r 1 ≤ δ (l) T (l) 1 + 2 x 0T (l) r 1 + x 0T (l) q 1 - l-1 i=1 x (i) 1q (i) . (15)
The second term of the RHS of (15) can be computed as a best k-term approximation error (see Definition 2), i.e.,

x 0T (l)

r 1 = σ k (W l x 0 ) 1 . (16)
For the two last terms on the RHS of (15), the triangle inequality yields x 0T (l)

q 1 - l-1 i=1 x (i) 1q (i) ≤ l-1 i=1 δ (i)
q (i) . Then, introducing this inequality with ( 13) and ( 16) in (15) gives, by keeping all the terms with δ (l)

T (l) r on the LHS, δ (l) T (l) r 1 ≤ 1 1 -γ γ δ (l) T (l) q 1 + 2σ k (W l x0) 1 + l-1 i=1 δ (i) q (i) . (17)
Thus, introducing (17) in ( 14) leads to

δ (l) 1 ≤ 1 + γ 1 -γ 2σ k (W l x0) 1 + l-1 i=1 δ (i) q (i) + δ (l) T (l) q 1 . (18)
By applying Lemma 5 to δ (j) , ∀j ∈ {1, . . . , l}, we have ∀i ∈ {1, . . . , n} , δ

Then, using these bounds in (18) leads to

δ (l) 1 ≤ 2(µ + 1) µ (4h -l + 1) σ k (W l x0) 1 + 1 4h -l + 1 l-1 i=1 δ (i) 1 . (20)
Now, we prove by induction that the inequality (12) holds ∀l ∈ {1, . . . , h}. Let denote x0j (i) the ith largest absolute value of W j x 0 . For l = 1, Theorem 1 in [START_REF] Fornasier | Compressive sensing[END_REF]

yields x0 -x (1) 1 1 ≤ 2(1 + γ) 1 -γ σ k (x0)1. (21)
Since σ k (x 0 ) 1 = σ k (W 1 x 0 ) 1 ≤ hx 01 (k + 1) ≤ hx 01 (1) = h W 1 x 0 ∞ , and (19), the result (21) leads to [START_REF] Novara | Sparse identification of nonlinear functions and parametric set membership optimality analysis[END_REF] with l = 1. Now, assume that ( 12) is true until l-1 with l ≥ 2. To prove that it is true for l, we need to bound the sum l-1 i=1 δ (i) 1 involved in (20). We can show, as detailed in [START_REF] Le | Hybrid dynamical system identification: geometry, sparsity and nonlinearities[END_REF], that each term in the sum is bounded by

δ (j) 1 < µ + 1 2µ W l x 0 ∞ , ∀j ∈ {1, . . . , l -1} (22)
via (20) and the fact that, for j < l ≤ h ≤ k,

σ k (W j x 0 ) 1 ≤ (h -j + 1) W l x 0 ∞ . (23)
In addition, similar steps can be used to complete the induction on l and prove that (12) holds ∀l ∈ {1, . . . , h}. Then, by using Proposition 1 we conclude that q (l) ∈ I 1 , ∀l ≤ h with h ≤ k.

Step 2. Now, we prove that in the (h+1)th iteration, solving (6) yields x 0 . By using Lemma 1, x 0 uniquely solves problem [START_REF] Donoho | Compressed sensing[END_REF] if for all nonzero δ ∈ Ker(A), the following condition holds: (24)

Indeed, the LHS of (24) can be rewritten as , and the NSP of order k (Definition 1) is applied with T = I 1 \ T (h+1) q and |T | = k. Then, we apply once more the NSP of order k with T = T (h+1) q and |T | = h ≤ k for the RHS of (25) to get i∈I1\T (h+1)

q |δ i | < γ i∈I0 |δ i | + γ 2 i∈I0 |δ i | + γ 2 i∈I1\T (h+1) q
|δ i |, which we rewrite as

i∈I1\T (h+1) q |δ i | < γ 1 -γ i∈I0 |δ i |. (26)
Thus, the assumption γ < guarantees that i∈I1\T (h+1) q |δ i | < i∈I0 |δ i | and that (24) holds. Hence, we obtain the unique solution x 0 in h + 1 iterations.

  i / ∈{q (l) } m l=1 |x 1i | is 0 and the solution

  x

Algorithm 2

 2 Estimation of a single parameter vector Require: A data set S = {(ϕ i , y i )} N i=1 and a number of iterations N s . -Initialize l = 0, W = I N . while l < N s do -Set l = l + 1. -Obtain θ * = arg min θ W e (θ) 1 , where e (θ) = y 1ϕ ⊤ 1 θ, . . . , y Nϕ ⊤ N θ ⊤ is the error vector.

Fig. 2 .

 2 Fig. 2. Parametric error (NPE) and estimated number of modes for different values of the threshold δ used in the switched system identification example.

  then, for a solution, x1 , of min x∈R n x T 1 , s.t. Ax = b, where T = T \i, we have x1i = 0.

≤ µ µ+1 δ (j) 1 .

 1 With j = i, i = 1 . . . l-1,


J Kristo Ersen

Francois Laroussinie

Kim G Larsen

Paul Pettersson

Wang Yi

A Compositional Proof of a Real{Time Mutual Exclusion Protocol

In this paper, we apply a compositional proof technique to an automatic veri cation of the correctness of Fischer's mutual exclusion protocol. It is demonstrated that the technique may avoid the state{ explosion problem. Our compositional technique has recently been implemented in a tool CMC 5 , which veri es the protocol for 50 processes within 172.3 seconds and using only 32MB main memory. In contrast all existing veri cation tools for timed systems will su er from the state{ explosion problem, and no tool has to our knowledge succeeded in verifying the protocol for more than 11 processes.

Introduction

It is well{known that the major problem in applying automatic veri cation techniques to analyze nite{state concurrent systems is the potential combinatorial explosion of the state space arising from parallel composition. In the last few years, there has been a number of automatic veri cation tools for real{time systems [START_REF] Bengtsson | Uppaal in 1995[END_REF][START_REF] Ho | Automated Analysis of an Audio Control Protocol[END_REF][START_REF] Daws | Verifying ET-LOTOS programs with KRO-NOS[END_REF]. Experiences with these tools show that the state{explosion problem is even more serious in verifying timed systems. As such a system must satisfy certain timing constraints on its behaviour, a model{checker must keep track of not only the part of state{space explored, but also timing information associated with each state (i.e. possible clock values), which is both time and space{consuming.

During the last decade, various techniques have been developed to avoid the state{explosion problem in verifying nite{state systems, either by symbolic representation of the states space using BDDs 5], by application of partial order methods [START_REF] Godefroid | A Partial Approach to Model Checking[END_REF][START_REF] Valmari | A Stubborn Attack on State Explosion[END_REF] which suppresses unnecessary interleavings of transitions, or by application of abstractions and symmetries 6, 7, 9]. These techniques have been further extended to deal with timed systems, e.g. [START_REF] Bengtsson | Uppaal in 1995[END_REF][START_REF] Ho | Automated Analysis of an Audio Control Protocol[END_REF], 17], 8]. However, when applying these techniques to parallel systems such as Fischer's protocol, a potential explosion in the global state{space remains. In 2], a compositional y Basic Research in Computer Science, Centre of the Danish National Research Foundation. 5 CMC: Compositional Model Checking veri cation technique is developed by Andersen 2] for nite{state systems. [START_REF] Laroussinie | Compositional Model Checking of Real Time Systems[END_REF][START_REF] Kim | Compositional and Symbolic Model-Checking of Real-Time Systems[END_REF], the technique has been further extended to deal with real{time systems modelled as networks of timed automata, which allows components of a real{time system to be gradually moved from the system description into the speci cation, thus avoiding any global state{space construction and even examination. Essential to the technique is that intermediate speci cations are kept small using e cient minimization heuristics.

In this paper, we apply this technique to give a compositional proof for Fischer's mutual exclusion protocol. In particular, it is shown that state{explosion is avoided in the veri cation of the protocol: the size of the correctness proof we o er only grows polynomially in the size of the number of processes in the protocol. A similar compositional technique has recently been implemented using C++ in a tool called CMC, Compositional Model Checking. This tool gives further experimental evidence of the potential of the technique: using only 172.3 seconds and 32MB main memory CMC automatically veri es the mutual exclusion property for the acyclic version of Fischer's protocol with 50 processes.

The paper is organized as follows: In the next section we brie y introduce our modelling and speci cation languages for real{time systems, and the formal description of Fischer's mutual exclusion protocol. Section 3 describes the compositional quotienting method and simpli cation techniques for logical formulas. In section 4, we present the proof for the mutual exclusion property of Fischer's protocol. In section 5 we report on the experimental results obtained using the CMC tool and compare the performance with that of our existing tool{suite 3]. Finally, in section 6 we give some concluding remarks and illustrate future work.

Real{Time Systems

In this section, we brie y introduce our modelling and speci cation languages for real{time systems, that have been studied previously in the literature, e.g. [START_REF] Yi | Automatic Veri cation of Real-Time Communicating Systems By Constraint-Solving[END_REF][START_REF] Laroussinie | Compositional Model Checking of Real Time Systems[END_REF][START_REF] Kim | Compositional and Symbolic Model-Checking of Real-Time Systems[END_REF][START_REF] Kim | Diagnostic Model-Checking for Real-Time Systems[END_REF]. For details, we refer to 15].

Models

We use timed transition systems as a basic semantic model for real-time systems. A timed transition system is a labelled transition system with two types of labels: atomic actions and delay actions (i.e. positive reals), representing discrete and continuous changes of real-time systems. Assume a nite set of actions Act ranged over by a; b etc, and a nite set of atomic propositions P ranged over by p; q etc. We use R to stand for the set of non-negative real numbers, for the set of delay actions f (d) j d 2 Rg, and L for the union Act . De nition 1. A timed transition system over Act and P is a tuple, S = hS; s 0 ; ! ; V i, where S is a set of states, s 0 is the initial state, ! S L S is a transition relation, and V : S ! 2 P is a proposition assignment function that for each state s 2 S assigns a set of atomic propositions V (s) that hold in s. u t

We use synchronization functions to describe concurrency and synchronizations between timed transition systems. A synchronization function f is a partial function (Act f0g) (Act f0g) , ! Act, where 0 denotes a distinguished noaction symbol 6 . Now, let S i = hS i ; s i;0 ; ! i ; V i i, i = 1; 2, be two timed transition systems and let f be a synchronization function. Then the parallel composition S 1 j f S 2 is the timed transition system hS; s 0 ; !; V i, where s 1 j f s 2 2 S whenever s 1 2 S 1 and s 2 2 S 2 , s 0 = s 1;0 j f s 2;0 , ! is inductively de ned as follows:

{ s 1 j f s 2 c ! s 0 1 j f s 0 2 if s 1 a ! 1 s 0 1 , s 2 b ! 2 s 0 2 and f(a; b) = c { s 1 j f s 2 (d) ! s 0 1 j f s 0 2 if s 1 (d) ! 1 s 0 1 and s 2 (d) ! 2 s 0 2
and nally, the proposition assignment function V is de ned by V (s

1 j f s 2 ) = V 1 (s 1 ) V 2 (s 2 ).
The type of systems we are studying is a particular class of timed transition systems that are syntactically described by networks of timed automata [START_REF] Yi | Automatic Veri cation of Real-Time Communicating Systems By Constraint-Solving[END_REF][START_REF] Laroussinie | Compositional Model Checking of Real Time Systems[END_REF][START_REF] Kim | Compositional and Symbolic Model-Checking of Real-Time Systems[END_REF][START_REF] Kim | Diagnostic Model-Checking for Real-Time Systems[END_REF]. A timed automaton 1] is a standard nite-state automaton extended with a nite collection of real-valued clocks. Let C be a nite set of real-valued clocks ranged over by x; y etc. We use B(C) ranged over by g (and latter D), to stand for the set of formulas that can be an atomic constraint of the form:

x n or x y n for x; y 2 l 0 is the initial node, E N B(C) Act 2 C N corresponds to the set of edges, and nally, V : N ! 2 P is a proposition assignment function. In the case, hl; g; a; r; l 0 i 2 E it is written, l g;a;r ! l 0 . u t

The semantics of a timed automaton is given in terms of clock assignments. A clock assignment u for C is a function from C to R. Let R C denote the set of clock assignments for C. For u 2 R C , x 2 C and d 2 R, u + d denotes the time assignment which maps each clock x in C to the value u(x) + d. For C 0 C, C 0 7 ! 0]u denotes the assignment for C which maps each clock in C 0 to the value 0 and agrees with u over CnC 0 . A semantical state of an automaton A is a pair (l; u) where l is a node of A and u a clock assignment for C. The initial state of A is (l 0 ; u 0 ) where u 0 is the initial clock assignment mapping all clocks in C to 0. The semantics of A is given by the timed transition system S A = hS; 0 ; !; V i, where S is the set of states of A, 0 is the initial state (l 0 ; u 0 ), ! is the transition relation de ned as follows:

{ (l; u) a !(l 0 ; u 0 ) if there exist r; g such that l g;a;r ! l 0 , g(u) and u 0 = r ! 0]u { (l; u) (d) !(l 0 ; u 0 ) if (l = l 0 ), u 0 = u + d and V is extended to S simply by V (l; u) = V (l).

hs; ui j = c ) c(u) hs; ui j = p ) p 2 V (s) hs; ui j = cp _ ' ) hs; ui j = cp or hs; ui j = ' hs; ui j = ' ^ ) hs; ui j = ' and hs; ui j = hs; ui j = 8 8' ) 8d; s 0 : s (d) ! s 0 ) hs 0 ; u + di j = ' hs; ui j = a] ' ) 8s 0 : s a ! s 0 ) hs 0 ; ui j = ' hs; ui j = x in ' ) hs; u 0 i j = ' where u 0 = fxg ! 0]u hs; ui j = Z ) hs; ui j = D(Z) Finally, for two timed automata A and B and a synchronization function f, the parallel composition A j f B denotes the timed transition system S A j f S B .

Speci cations

To specify safety and bounded liveness properties of timed systems, we use the timed modal logic L s , studied in [START_REF] Laroussinie | From Timed Automata to Logic | and Back[END_REF][START_REF] Kim | Compositional and Symbolic Model-Checking of Real-Time Systems[END_REF][START_REF] Kim | Diagnostic Model-Checking for Real-Time Systems[END_REF]. Let K be a nite set of clocks, called formula clocks, and Id a set of identi ers. The set of formulas of L s over K, Id, Act, and P is generated by the following syntax with ' and ranging over L s :

' ::= cp j cp _ ' j ' ^ j 8 8' j a] ' j z in ' j Z where cp may be an atomic clock constraint c in the form of x n or x y n for x; y 2 K and natural number n, or an atomic proposition p 2 P, a 2 Act (an action), z 2 K and Z 2 Id (an identi er). The meaning of the identi ers is speci ed by a declaration D assigning a formula of L s to each identi er. When D is understood we write Z def = ' for D(Z) = '.

Given a timed transition system S = hS; s 0 ; !; V i described by a network of timed automata, the L s formulas are interpreted in terms of an extended state hs; ui where s 2 S is a state of a timed transition system, and u is a clock assignment for K.

Let D be a declaration. Formally, the satisfaction relation j = D between extended states and formulas is de ned as the largest relation satisfying the implications of Table 1. For simplicity, we shall omit the index D and write j = instead of j = D whenever it is understood from the context. Finally, a network of timed automata A satis es a formula ' written A j = ' when h(l 0 ; u 0 ); v 0 i j = ' where l 0 is the initial node of A, and u 0 and v 0 are the assignments with u 0 (x) = 0 for all automaton clocks x and v 0 (z) = 0 for all formula clocks z. Note that (l 0 ; u 0 ) is the initial state of A.

Fischer's Protocol Revisited

As an example of networks of timed automata, we study Fischer's mutual exclusion protocol. The reason for choosing this example is that it is well{known and well{studied by researchers in the context of real{time veri cation. More importantly, the size of the example can be easily scaled up by simply increasing The protocol is to guarantee mutual exclusion in a concurrent system consisting of a number of processes, using clock constraints and a shared variable. We shall model each of the processes as a timed automaton, and the protocol as a network of timed automata. Each of the processes is assumed to have a local clock. The idea behind the protocol is that the timing constraints on the local clocks are set so that all processes can change the global variable to its own process number, then read the global variable later and if the shared variable is still equal to its own number, enter the critical section. Each process P i with i being its identi er, has a clock x i . Let A k = f:= i j i = k + 1:::ng, T k = f= i j i = k + 1:::ng, F k = f6 = i j i = k + 1:::ng, and S k = A k T k F k . We model the shared variable as a timed automaton V over the set of atomic actions S 0 f:= 0; = 0g, where V = hN; h 0 ; E; V i with N = fV 0 :::V n g, h 0 = V 0 , E = fhV i ; tt; := j; ;; V j i j i; j = 0:::ng fhV i ; tt; = i; ;; V i i j i = 0:::ng fhV i ; tt; 6 = j; ;; V i i j i 6 = jg, and we simply assume V is de ned by V (V i ) = ; for all i n. The automaton for a typical process P i is shown in Fig 1.

A i tt fx i g = 0 6 = i x i < 1 x i > 2 := i = i fx i g fg := 0 tt tt fg CS i C i B i fg
We assume that the proposition assignment function is de ned in such a way that at(l 0 ) 2 V (l) if l 0 = l and :at(l 0 ) 2 V (l) if l 0 6 = l for all nodes l and l 0 . Now, the whole protocol is described as the following network:

FISCHER n (P 1 j f1 (P 2 j f2 (P 3 j f3 :::j fn 1 P n ):::)j g V where j fi and j g are the interleaving and full synchronization operators, induced by synchronization functions f i and g respectively, de ned by f i (a; 0) = a when a 2 f:= i; = i; 6 = ig and f i (0; a) = a when a 2 S i , and g(a; a) = a. Note that in P i j fi (P fi+1 :::), P i is allowed to perform f:= i; = i; 6 = ig and the righthand side is allowed to perform all actions with indices higher than i that is, S i .

Intuitively, the protocol behaves as follows: The constraints on the shared variable V ensure that a process must reach B{node before any process reaches C{node; otherwise, it will never move from A{node to B{node. The timing constraints on the clocks ensure that all processes in C{nodes must wait until all processes in B{nodes reach C{nodes. The last process that reaches C{node and sets V to its own identi er gets the right to enter its critical section.

We need to verify that there will never be more than one process in its critical section. An instance of this general requirement can be formalized as an invariant property: M 12 (:at(CS 1 )_ :at(CS 2 )) ^Va2S0 a] M 12 ^8 8M 12 So we need to prove the theorem FISCHER n j = M 12 3 Compositional Model{Checking Model{checking of real{time systems may be carried out in a symbolic fashion e.g. 11, 19]. However, when applying these techniques to parallel networks such as FISCHER n a potential explosions in the global symbolic state{space may seriously hamper the technique.

In 13, 15] we presented a compositional veri cation technique, which allows components of a real{time system to be gradually moved from the system description into the speci cation, thus avoiding any global state{space construction and even examination. Essential to the technique is that intermediate speci cations are kept small using e cient minimization heuristics. Our technique may be seen as a real{time extension of the compositional technique presented and experimentally applied by Andersen 2] for ordinary nite{state systems. In this section we give a brief review of the technique in 13, 15].

Quotient Construction

The main ingredient in our compositional veri cation technique is the so{called quotient construction, which allows components of a network to be moved into the speci cation. More precisely, given a formula ', and two timed automata A and B we may construct a formula (called the quotient) ' f B such that

A j f B j = ' if and only if A j = ' f B (1)
The bi{implication indicates that we are moving parts of the parallel system into the formula. Clearly, if the quotient is not much larger than the original formula, we have simpli ed the task of model{checking, as the (symbolic) semantics of A is signi cantly smaller than that of A j f B. More precisely, whenever ' is a formula over K, B is a timed automaton over C and l is a node of B, we de ne the quotient formula ' f l over C K in Table 2 on the structure of ' 7 8 . Note that the quotient construction for identi ers introduces new identi ers of the form X l . The new identi ers and their de nitions are collected in the (quotient) declaration D B . We recall from 15] the following important theorem, which justi es the construction: 7 For g = c1 ^: : : cn a clock constraint we write g ) ' as an abbreviation for the formula :c1 _ : : : _ :cn _ '. This is an Ls{formula as atomic constraint are closed under negation. 8 In the rule for a]', we assume that all nodes l of a timed automaton are extended with a 0{edge l tt;0;; ! l. Theorem 3. Let A and B be two timed automata and let l 0 be the initial node of B. Then A j f B j = D ' if and only if A j = DB ' f l 0

c f l = c p f l = t t ; p 2 V (l) p ; p 6 2 V (l) ('1 ^'2) f l = ('1 f l) ^('2 f l) (8 8') f l = 8 8 ' f l (x in ') f l = x in (' f l) (c _ ') f l = (c f l) _ (' f l) (p _ ') f l = (p f l) _ (' f l) X f l = X l where X l def = D(X) f l ( a]') f l = l g;c;r ! l 0 ^f(b; c) = a g ) b](r in ' f l 0 )

Minimizations

It is obvious that repeated quotienting leads to an explosion in the formula (in particular in the number of identi ers). The crucial observation made by Andersen in the (untimed) nite{state case is that simple and e ective transformations of the formulas in practice may lead to signi cant reductions.

In presence of real{time we need, in addition to the minimization strategies of Andersen, heuristics for propagating and eliminating constraints on clocks in formulas and declarations. Below we describe the transformations considered:

Reachability: When considering an initial quotient formula ' f l 0 not all identi ers in D B may be reachable. Application of an \on-the-y" technique will insure that only the reachable part of D B is generated.

Boolean Simpli cation Formulas may be simpli ed using the following simple boolean equations and their duals: ff ^' ff, tt ^' ', x in ff ff.

Constraint Propagation: Constraints on formula clocks may be propagated using various distribution laws (see Table 3). In some cases, propagation will lead to trivial clock constraints, which may be simpli ed to either tt or ff and hence made applicable to Boolean Simpli cation. As can be seen in Table 3 certain operations are to be performed on constraints during propagation. These operations include the following: D " = fu + d j u 2 D and d 2 Rg frgD = f r 7 ! 0]u j u 2 Dg D # = fu j9d 2 R : u + d 2 Dg It may be shown that the set of constraints B(K) is closed under the above operations, and that they together with inclusion{ and emptyness{checking may be computed e ciently (in cubic time in the number of clocks) (see e.g. 15]).

;

) ' t t D ) c t t ; if D c D ) ( a]') a](D ) ') D ) ('1 ^'2) (D ) '1) ^(D ) '2) D ) (x in ') x in (fxgD ) ') D ) (p _ ') p _ (D ) ') D ) (c _ ') (D ^:c) ) ' D ) (8 8') 8 8(D " ) ') ; if D # D D ) X D ) D(X) Table 3. Constraint Propagation
Constant Propagation: Identi ers with identi er-free de nitions (i.e. constants such as tt or ff) may be removed while substituting their de nitions in the declaration of all other identi ers. Trivial Equation Elimination: Equations of the form X def = a]X are easily seen to have X = tt as solution and may thus be removed. More generally, let S be the largest set of identi ers such that whenever X 2 S and X def = ' then ' tt=S] 9 can be simpli ed to tt. Then all identi ers of S can be removed provided the value tt is propagated to all uses of identi ers from S (as under Constant Propagation). The maximal set S may be e ciently computed using standard xed point computation algorithms. Equivalence Reduction: If two identi ers X and Y are semantically equivalent (i.e. are satis ed by the same timed transition systems) we may collapse them into a single identi er and thus obtain reduction. However, semantical equivalence is computationally very hard 10 . To obtain a cost e ective strategy we approximate semantical equivalence of identi ers as follows: Let R be an equivalence relation on identi ers. R may be extended homomorphically to formulas in the obvious manner: i.e. ('

1 ^'2 )R(# 1 ^#2 ) if ' 1 R# 1 and ' 2 R# 2 , (xin')R(xin#)
and a]'R a]# if 'R# and so on. Now let = be the maximal equivalence relation on identi ers such that whenever X = Y , X def = ' and Y def = # then ' = #. Then = provides the desired cost e ective approximation: whenever X = Y then X and Y are indeed semantically equivalent. Moreover, = may be e ciently computed using standard xed point computation algorithms.

Fischers Protocol

From section 2 we recall that the protocol FISCHER n consists of n processes P 1 : : : P n competing for a critical section by setting and testing a shared variable V , and that the mutual exclusion property we verify is that P 1 and P 2 cannot be in their critical section at the same time, i.e: M 12 (:at(CS 1 )_ :at(CS 2 )) ^Va2S0 a] M 12 ^8 8M 12

In the remainder of this section we shall apply our compositional model{checking technique to verify the protocol. Our observation is that by rst quotienting away V; P 1 and P 2 the quotient hereby obtained simpli es to tt under our minimization heuristics. Thus no examination of the components P 3 ; : : : ; P n is required: regardless of their behaviour the mutual exclusion property M 12 will be satis ed.

In other words, state{space explosion is avoided as it is su cient to explore only a xed part of the system to prove the desired property.

Constructing the Quotient

The order by which components of a network is quotiented out may highly determine the success of our method (this resembles the importance of variable{ ordering in the BDD technology). Here, we choose to rst quotient out the variable V followed by the relevant processes P 1 and P 2 , while of course constantly minimizing the intermediate equation systems as much as possible.

Step 1: In the rst step we remove the variable V from the network and transform M 12 by quotienting it with the locations V 0 ; : : : ; V n . This will result in an equation system containing n + 1 identi ers X 0 ; : : : ; X n where X i denotes the quotient M 12 = g V i .

As the synchronization function g between V and the rest of the system is de ned as g(a; a) = a for all possible action transitions a the quotient will have exactly same conjuncts as M 12 . Further as V does in all of its locations satis es neither :at(CS 1 )nor :at(CS 2 )we get the following family of formulae X i , where i = 0; : : : ; n: X i = (:at(CS 1 )_ :at(CS 2 )) ^ = i] X i V j := j] X j ^Vj6 =i 6 = j] X i ^8 8X i : This new equation system (i.e. the top identi er X 0 ) constitutes the requirement for the remaining components P 1 ; : : : ; P n . The identi er X i expresses the requirement to the remaining system when the variable holds the value i. That is, (:at(CS 1 )_ :at(CS 2 )) should still be satis ed, and as long as the variable is only tested upon or as long as time passes X i should still hold. If the variable is set to another value j the formula de ned by X j should hold instead.

Step 2: As (:at(CS 1 )_ :at(CS 2 ))is required by all identi ers and their de nitions di er slightly the equation system cannot be simpli ed any further. Thus we proceed to transform the equation system with respect to removal of P 1 from the network. The quotient operator used to do this will be subscripted with the synchronization function f 1 . In the following we will drop the synchronization function as subscript to the quotient operator, as it is obvious which one is used. As the equation system after step 1 contains n + 1 equations and P 1 has four control locations the new equation system will contain 4 (n + 1) equations.

For each j = 0; : : : ; n we compute X j =l, where l 2 fA 1 ; B 1 ; C 1 ; CS 1 g. The three cases where j = 0; 1; 2 are treated separately, while the remaining cases are treated together. When quotienting any of the identi ers X i with A 1 , B 1 or C 1 the requirement (:at(CS 1 )_ :at(CS 2 )) disappears because :at(CS 1 ) is satis ed in all three locations. When quotienting any of the identi ers X i with CS 1 , (:at(CS 1 )_ :at(CS 2 )) remains in the de nition of the new identi er as neither :at(CS 1 ) nor :at(CS 2 ) is a satis ed by CS 1 . Due to lack of space we do not display this quotient, instead we continue the quotienting with respect to P 2 and therefore calculate M 12 =V 0 =A 1 =A 2 .

Step 3: The equation system of M 12 =V 0 =A 1 =A 2 consists of 4 4 (n+1) equations, namely the size of the product automaton of V , P 1 and P 2 . The equations can be grouped as 16 equations resulting from X 0 =P 1 =P 2 , 16 equations resulting from X 1 =P 1 =P 2 , 16 equations resulting from X 2 =P 1 =P 2 and nally 16 (n 2) equations resulting from X j =P 1 =P 2 where j = 3; : : : ; n. For a xed choice of locations, l 1 and l 2 in P 1 and P 2 the set of identi ers X j =l 1 =l 2 for j = 3; : : : ; n will describe very similar properties.

The equation system is presented as a formula graph, and part of it appears in Figure 2. Each node represents a formula identi er, and outgoing edges represents conjuncts in the de nition of an identi er. For instance, the upper most node in the graph, re ects that: X 0 =A 1 =A 2 = x 1 in (X 0 =B 1 =A 2 )^ = 0] (X 0 =A 1 =A 2 ) ^: : : : An atomic proposition (possibly a disjunction) labelling a node means that this atomic proposition appears as a conjunct in the de nition of the identi er the node represents. Hence, (:at(CS 1 )_ :at(CS 2 )) is a conjunct in the de ning equation of X 2 =CS 1 =CS 2 .

In the quotient all identi ers have a conjunct which refers to the identi er itself through the 88{modality. That is, For all Y the de nition is on the form Y = : : : ^8 8Y ^: : :. In the formula graph this would appear as self loops labelled with the 88{modality in all nodes, but in order to keep the graph simple we have omitted these loops. Further the quotient is symmetrical as P 1 and P 2 are symmetrical up to names on locations and clocks, therefore we only display half of the quotient as a formula graph.

To obtain a compact representation in Figure 2 we have used the the following abbreviations. A grey node labelled X j =l 1 =l 2 where l 1 ; l 2 are locations of P 1 and P 2 abbreviates the whole family of nodes X 3 =l 1 =l 2 ; : : : ; X n =l 1 =l 2 .

Similarly, edges labelled := j or = j really represents a whole family of edges namely one edge for each choice of j = 3; : : : ; n. E.g. the := j labelled edge from X 0 =A 1 =A 2 to X j =A 1 =A 2 in Figure 2 represents the family X 0 =A 1 =A 2 :=3 ! X 3 =A 1 =A 2 ; : : : ; X 0 =A 1 =A 2 :=n ! X n =A 1 =A 2 .

The overall structure of the formula graph for the resulting quotient is shown in Figure 3. Six typical parts of the quotient can be identi ed, these parts are labelled 1, 2, 3, 4, 5 and 6.

Part 1 of the quotient results from keeping P 2 xed in its initial location A 2 and letting P 1 and the variable V vary as much as they can. Not surprisingly this part of the quotient reduces to tt. We will later argue formally why this is actually the case.

Part 2 of the quotient corresponds to the behaviour part where rst P 1 assigns the variable, then P 2 assigns it, where after P 2 enters the critical section and hence P 1 fails to observe the variable having the value 1 and it returns to its initial state A 1 . :=j Fig. 2. Formula (sub-)Graph for M 12 =X 0 =A 1 =A 2 .

Xj/C1/A2 Xj/B1/B2 X1/C1/B2 Xj/B1/A2 =j :=j =j :=j =j :=j X2/C1/A2 X2/A1/CS2 Xj/A1/A2 =0 =0 X0/A1/A2 X0/B1/A2 X0/
Part 3 of the quotient is where P 1 and P 2 are in the critical section at the same time. The concrete manifestation of this is that the formula identi ers of this part requires (:at(CS 1 )_ :at(CS 2 )) to be satis ed by the remaining components P 3 ; : : : ; P n . It is essential to the proof of the correctness that this part of the quotient will not be required to hold for the network of processes P 3 ; : : : ; P n .

The actual proof of this relies on the use of constraint propagation: We show that from the initial clock constraint (all clocks having value 0) this dangerous part of the quotient cannot be reached. Part 4 is symmetrical to part 2 and part 5 is symmetrical to part 1. The Fig. 3. Overall structure of Formula Graph last part of the quotient, the one numbered 6, consists of the before mentioned identi ers X j =l 1 =l 2 where l 1 is a location in P 1 , l 2 is a location in P 2 and j = 3; : : : ; n. This part of the quotient is the requirement when V takes a value di erent from 0; 1 and 2.

Simpli cation

The quotient formula M 12 =V 0 =A 1 =A 2 illustrated in gure 2 is according to Theorem 3 the necessary and su cient property of the remaining components P 3 ; : : : ; P n in order that the overall system FISCHER n satis es M 12 . We may now apply our simpli cation heuristics.

To our pleasant surprise we observe the quotient formula M 12 =V 0 =A 1 =A 2 calculated above simpli es to tt when rst applying Constraint Propagation followed by Trivial Equation Elimination. Therefore we do not have to perform quotienting with respect to the remaining components in the protocol, and hence we may conclude that an increase in the number of components in the protocol only gives rise to a polynomial growth in the size of the proof. Applying constraint Propagation reveals the fact that none of the identiers X j =CS 1 =CS 2 where j = 1; : : : ; n can be reached from the initial constraint.

As none of the remaining nodes in the graph contains propositions or clock constraints Trivial Equation Elimination will immediately reduce all identi ers and especially the top identi er M 12 =V 0 =A 1 =A 2 to tt.

Constraint Propagation can be implemented on our formula graphs in the following manner: Whenever X g; ;r ! Y is an edge in the graph and we consider an implication D ) X, the constraint D may be propagated using the rewrite rules of Table 3 to a constraint on Y represented by the implication:

(frg(D ^g)) " ) Y:

(2) Thus constraint propagation in a general formula graph, where a node can have multiple outgoing edges will result in a conjunction of formulas of the type in [START_REF] Andersen | Partial Model Checking[END_REF]. What we intend to do here, however, is to direct the propagation of constraints along a speci c path in the formula graph towards speci c identi ers that we wish to prove unreachable. To this speci c purpose we introduce the notion of guided Constraint Propagation. In a guided Constraint Propagation we simply focus on a speci c path in the formula graph and disregard all other edges.

In the following we perform such a guided constraint propagation towards part 3 of the quotient by following a path through (X 0 =A 1 =A 2 ), (X 0 =A 1 =B 2 ), (X 0 =B 1 =B 2 ), (X 1 =C 1 =B 2 ), (X 1 =CS 1 =B 2 ), (X 2 =CS 1 =C 2 ), (X 2 =CS 1 =CS 2 ), see Figure 2, and discover that (X 2 =CS 1 =C 2 ) is hit by the empty constraint and thus its reference to (X 2 =CS 1 =CS 2 ) has no importance in practice.

In the propagation we jump directly to the situation where the node (X 0 =B 1 =B 2 ) has been reached by letting time pass while resetting rst the clock x 2 and then 1 . In other words we consider the implication (x 2 > x 1 ) ) (X 0 =B 1 =B 2 ): Using (2) we may propagate with respect to the edge X 0 =B 1 =B 2 x1<1; ;fx1g

! X 1 =C 1 =B 2 yielding x 2 > x 1 ) X 1 =C 1 =B 2 : Now propagating this with respect to the edge X 1 =C 1 =B 2 x1>2; ;; ! X 2 =CS 1 =B 2 yields (x 2 > x 2 x1 > 2) ) X 2 =CS 1 =B 2 : Finally propagating this constraint with respect to X 2 =CS 1 =B 2 x2<1; ;fx2g ! X 2 =CS 1 =C 2 we get x 2 in (x 2 > x 1 ^x1 > 2 ^x2 < 1) ) X 2 =CS 1 =C 2 :
Clearly the constraint (x 2 > x 1 ^x1 > 2 ^x2 < 1) is empty and hence the whole propagation simpli es to t t.

By performing this form of guided Constraint Propagation we can prove that none of the formula identi ers in the quotient requiring P 1 or P 2 not to be in the critical section are reachable from the initial time zone. Of course we can also propagate constraints to all the other parts of the quotient, but this will not reduce the quotient as all other parts really are reachable. Trivial Equation Elimination reduces all remaining identi ers to t t as they are de ned by righthand sides which after Constraint Propagation are entirely built from the following connectives: t t; g ); ^; 8 8.

Experiments

The quotient construction together with the simpli cation techniques presented in the previous section have been implemented with C++ in a prototype tool called CMC (Compositional Model-Checking) 11 . CMC enables us to compute the quotient of an L s formula with respect to a timed automaton and then to simplify the quotient using our simpli cation. In fact, CMC enables quotienting with respect to formulas of the richer logic L 14] which allows general disjunction and existential modalities (99, hai). All simpli cation techniques of L s can be applied (and have been implemented in CMC) to L with the exception that no constraint propagation has been given for general disjunction and the existential modalities.

A few new simpli cation strategies, which are quite useful in an actual veri cation, have been introduced. One of these is reduction with respect to so{called hit{zones, which essentially is an exhaustive constraint propagation 11

In the near future CMC will be integrated in and available through the tool suite Uppaal 3].

providing the automatic counterpart to the so{called guided constraint propagation used in the previous section. The idea behind this simpli cation is to precompute, for any variable, the domain (in terms of clock constraints) in which the variable will be considered during a given veri cation. Given these domains, called hit{zones, it is possible in several cases to simplify clock constraints to either `true' or `false' (and hence amenable to constant propagation). Another simpli cation which is performed by the program is to replace any variable X with the following form: X = : : : ^y < k ^8 8X by `false'.

In our experimental investigation we have compared the current version of the tool CMC with the performances of both the backward and forward reachability checker of Uppaal on an acyclic version of Fischer's protocol. During the experiment both CMC and Uppaal was installed on a machine running SunOS This paper has successfully demonstrated that the compositional proof technique of 15] may avoid the state{explosion problem. In particular, it has been shown that state{explosion is avoided in the veri cation of Fischer's protocol: the size of the correctness proof we o er grows only polynomially in the size of the number of processes in the protocol. Furthermore, this claim has been given experimental evidence by the tool CMC, which manages veri cation of 50 processes. In contrast all exiting veri cation tools will su er from state{explosion, and no tools has succeeded in verifying the protocol for more than 11 processes.

Immediate future work includes integration of the CMC implementation in the veri cation tool Uppaal which will require certain extensions as Uppaal allows integer variables as well as clocks with interval{bounded slopes. Also, the optimal order in which components are factored out needs a better understanding. This resembles the situation in BDDs, where the ordering of propositional variables strongly in uences the size of the BDD. Our ambition for future work is to get a better understanding of when and how well our technique will work.

Fig. 1 .

 1 Fig. 1. Fischers Protocol for Mutual Exclusion.

5. 5

 5 with 32MB of primary memory and 128 of swap memory. Previously the backward reachability tool of Uppaal has been demonstrated advantageous in a comparison with other veri cation tools 15] on this version of Fischer's protocol. However, as can be seen by the outcome of the present experiment in

Table 1 .

 1 De nition of satis ability.

Table 2 .

 2 De nition of Quotient '

f l

Table 4 ,

 4 Uppaal is clearly outperformed by CMC, which manages veri cation of 50 processes. tool n #-processes 2 3 4 5 6 7 8 9 10 20 30 40 50 Uppaal forwards 0.2 0.2 0.9 10.7 244.4 ? Uppaal backwards 0.1 0.2 0.3 1.2 6.2 38.5 310.6 ? CMC 0.2 0.4 0.6 0.8 1.2 1.6 2.0 2.5 3.2 14.5 40.0 88.2 172.3

Table 4 .

 4 Test results 6 Conclusion

We extend the transition relation of a timed transition system such that s 0 ! s 0 i s = s 0 .

' t t=S ] is the formula obtained by substituting all occurrences of identi ers from S in ' with the formula t t.

For the recursion{free, untimed part of the logic semantical equivalence is already NP{complete.