Julien Tierny
email: tierny@lifl.fr

Jean-Philippe Vandeborre
email: vandeborre@lifl.fr

Mohamed Daoudi
email: daoudi@lifl.fr

Topology driven 3D mesh hierarchical segmentation

. This high level information drives both the feature boundary computation as well as the feature hierarchy definition. Proposed hierarchical scheme is based on the key idea that the topology of a feature is a more important decomposition criterion than its geometry.

First, the enhanced topological skeleton of the input triangulated surface is constructed. Then it is used to delimit the core of the object and to identify junction areas. This second step results in a fine segmentation of the object. Finally, a fine to coarse strategy enables a semanticoriented hierarchical composition of features, subdividing human limbs into arms and hands for example.

. Thanks to the high level description it uses as an input, presented approach results, with low computation times, in robust and meaningful compatible hierarchical decompositions.

Introduction

Mesh segmentation consists in subdividing a polygonal surface into patches of uniform properties [START_REF] Shamir | Segmentation and shape extraction of 3D boundary meshes[END_REF][START_REF] Attene | Mesh segmentation: A comparative study[END_REF], either from a strictly geometrical point of view or from a perceptual point of view. This operation has become a necessary pre-processing tool for many applications such as shape modeling [START_REF] Funkhouser | Modeling by example[END_REF], deformation [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF], compression, texture mapping [START_REF] Zhang | Feature-based surface parametrization and texture mapping[END_REF], retrieval, etc.

According to Attene et al. [START_REF] Attene | Mesh segmentation: A comparative study[END_REF], mesh segmentation techniques can be classified into two categories. On the one hand, geometry-oriented methods aim at defining patches of homogeneous geometry. In this case, algorithms are driven by purely low level geometrical information, such as curvature [START_REF] Yamauchi | Mesh segmentation driven by gaussian curvature[END_REF]. On the other hand, semantic-oriented methods aim at distinguishing regions of perceptual interest, following high level notions such as defined in human percep-tion theory [START_REF] Biederman | Recognition-by-components: A theory of human image understanding[END_REF]. This kind of approach is particularly suited for human shape interaction applications (texture mapping, modeling, deformation, etc.), where the decomposition has to be meaningful from a human user's point of view.

In this latter case, three main difficulties arise. (i) Characterizing high level notions from low-level measurements (such as curvature or geodesic distance evaluations) remains an open issue [START_REF] Attene | Mesh segmentation: A comparative study[END_REF]. (ii) For most applications, a hierarchical decomposition scheme is expected, so as to provide a progressive understanding of the object, but it is often based on low-level measurement clustering [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF][START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF]. (iii) The compatibility of the segmentation -which means the ability to identically segment objects that are visually similar -has been recently expressed as a new challenging problem [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF].

To address these issues, Katz et al. [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF], using multidimensional scaling, introduce the high level notion of feature points -vertices located on the extremity of prominent components -to drive the segmentation of the object into core and features. Then, only patches that contain feature points are recursively sub-divided. This means that, at the finest levels of the hierarchy, the more the patches will be located near the feature points, the smaller they will be (fig. 5), which is not a natural hierarchical scheme.

Podolak et al. [START_REF] Podolak | A planar-reflective symmetry transform for 3D shapes[END_REF] propose to drive the segmentation by the notion of symmetry. They describe a clustering-based algorithm using planar-reflective symmetry transform values as an input. However, the depth of their hierarchical scheme is user-defined. Moreover, the symmetry of an object depends on its pose, which impacts the global pose robustness of the overall approach.

Li et al. [START_REF] Li | Decomposing polygon meshes for interactive applications[END_REF] propose to construct a skeleton of the object and to use this high level structural information to drive the segmentation. Lien et al. [START_REF] Lien | Simultaneous shape decomposition and skeletonization[END_REF] extend this idea constructing a hierarchical skeleton, analyzing the principal axis of its convex hull. However, this method is reported to be time consuming when dealing with non-null genus surfaces. Moreover, the convex hull of an object is dependent on its pose, which also impacts the pose robustness of the overall algorithm.

Berreti et al. [START_REF] Berreti | Partitioning of 3D meshes using Reeb graphs[END_REF] propose a method that overcomes above issue using topological skeletons [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF], defined with regard to a geodesic based mapping function computed on the surface. However, their algorithm is not hierarchical and only segments objects into core and limbs. Zhang et al. [START_REF] Zhang | Feature-based surface parametrization and texture mapping[END_REF] use topological skeletons as well but also propose to segment obtained patches along separating regions, where the area of the level sets of the mapping function suddenly varies.

In this paper, we propose to address the three previously mentioned segmentation related issues. Firstly, as the structure of an object is an important perceived shape characteristic, we propose to drive the segmentation by the topology of the object, using an enhanced version of topological skeletons [START_REF] Tierny | 3D mesh skeleton extraction using topological and geometrical analyses[END_REF]. Secondly, as this shape abstraction also encodes shape geometry, we propose a semantic-oriented hierarchical segmentation process that gathers object's features according to their topology and their geometry. Thirdly, as segmentation is topology driven, objects of compatible topology are segmented in a compatible manner.

This paper makes the following contributions. In the next section, we show how to adapt shape topological description techniques to the segmentation problem, with an appropriate surface mapping function. In the third section, we show how to use this high level information for the feature boundary definition. In the fourth section, we show how to benefit from the topological description of the shape to achieve a semantic-oriented feature hierarchy. Experimental results are presented in the fifth section, where method performance is evaluated according to criteria enumerated in latest surveys [START_REF] Attene | Mesh segmentation: A comparative study[END_REF].

Fitting enhanced topological skeletons to the segmentation problem

In the first stage of the method, given an input triangulated surface T (manifold, connected and closed), an enhanced topological skeleton is extracted. The following paragraph sums up previous work [START_REF] Tierny | 3D mesh skeleton extraction using topological and geometrical analyses[END_REF] presenting this process.

Enhanced topological skeleton extraction

First, mesh feature points (in pink in figure 1(a)) are automatically extracted, intersecting geodesic-based map ex- trema [START_REF] Tierny | 3D mesh skeleton extraction using topological and geometrical analyses[END_REF]. Then, for each vertex in the mesh, a mapping function f m is defined as the geodesic distance to the closest feature point (fig. 2(a)). Next, for each vertex v ∈ T , an upper-value approximation of f -1 m ( f m (v)), noted Γ(v), is computed along the edges of T . In particular, the connected component of Γ(v) containing v is identified and noted γ(v). Analyzing the evolution of the number of connected subsets of Γ(v) as f m evolves enables the construction of a Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF] (fig. 1(b)). At this stage, each connected component of the Reeb graph is modeled with an ordered collection of closed curves γ(v). The next step consists in identifying constrictions [START_REF] Hétroy | Constriction computation using surface curvature[END_REF] within these collections. For each connected component of the Reeb graph, the average Gaussian curvature on each curve γ(v) is computed. Then, local negative minima are identified as constrictions (fig. 1(c)). Finally, the connected components of the Reeb graph are subdivided using these constrictions as boundaries between subparts (fig. 1(d)). As a conclusion of this algorithm, the input surface is represented by an enhanced topological skeleton (fig. 1(d)), which encodes the topological and geometrical evolution of the contours of the mapping function f m .

Constriction location optimization

In this paper, we use such a high level description to drive the segmentation process. In particular, as identified constrictions are located on the narrowest parts of the surface, they are good candidates for feature boundary definition (fig. 1(c)). However, in the latter algorithm, constrictions are identified along level line approximations, which do not necessarily follow the concavity of the surface. Within the framework of segmentation, to conform to human perception theory [START_REF] Biederman | Recognition-by-components: A theory of human image understanding[END_REF], feature boundaries must be aligned with concave areas. Consequently, we force the alignment of the mapping function level lines with surface bottlenecks, by integrating surface curvature into the geodesic distance computation. Such a curvature constrained geodesic distance evaluation is computed with Dijkstra's algorithm, us-ing w(v i , v j ) as weight for every edge (v i , v j ):

w(v i , v j ) = e(v i , v j ) e + α × |I c (v i ) -I c (v j )| ∆I c (1)
where e(v i , v j ) is the euclidean distance between v i and v j , e the average edge length in T , α a predefined ratio, I c (v) the curvature index [START_REF] Koenderink | Surface shape and curvature scales[END_REF] in v (I c (v) ∈ [-1; 1]) and ∆I c the average curvature index difference between two adjacent vertices in T . The motivation of this computation is to increase the distance between two vertices when a concave region separates them. Consequently, when I c (v i ) or I c (v j ) (or both) are negative (which corresponds to concavity), α is set to α 0 and 0 otherwise. Therefore, level lines quickly visit convex areas and slowly get aligned with concave ones. This effect can be observed in figure 2, which shows on a hand model the mapping function computed in [START_REF] Tierny | 3D mesh skeleton extraction using topological and geometrical analyses[END_REF] (defined as the geodesic distance to the closest feature point), the curvature index distribution and the mapping function used in this paper (defined as the curvature constrained geodesic distance to the closest feature point, noted f m ). In particular, the reader can notice in figure 2(d) that displayed level sets are aligned with the concave parts of the finger. In these illustrations, α 0 has been set to a high value for display purpose. In practice, it has been set to 0.05 for every model. Once this curvature constrained function is computed, the Reeb graph is constructed. Then, for each connected component of the graph, a curvature estimation ζ(γ(v)) is computed for each of its curve γ(v) as follows:

ζ(γ(v)) = ∑ v i ∈γ(v) I c (v i ) × (L e 1 (v i ) + L e 2 (v i )) 2 × P (γ(v)) (2)
where P (γ(v)) is the perimeter of γ(v) and where L e 1 (v i ) and L e 2 (v i ) are the lengths of the edges adjacent to v i on γ(v). This computation gives a robust and relevant estimation of the curvature along each curve γ(v). Then, curves that locally minimize ζ(γ(v)) are identified as constrictions.

In this section, we introduced an adaptation of shape topological description techniques to the segmentation problem. We proposed a novel mapping function, based on a modified computation of geodesic distances. This computation integrates curvature so that the mapping function level lines are forced to be located along concave areas. Then, estimating the curvature on these collections of curves enables the detection of the most concave ones, that we refer as constrictions. In the next section, these constrictions are used for feature boundary computations.

Feature boundary computation

In the previous stage of the approach, a modified version of enhanced topological skeletons has been computed. Each node of this skeleton references a surface patch of the mesh. In figure 3 a distinctive color, resulting in an over-segmentation of the object. In this section, we describe a region merging algorithm based on the notions of core and junction areas, which results in a fine segmentation of the object (fig. 3(d)).

Junction areas

We classify the nodes of the skeleton into three categories, according to their degree. A node N is an extremity node if its degree (deg(N)) equals 1 (related surface patch contains a feature point, in green in fig. 3(a)). A node N is a tubular node if deg(N) = 2 (in blue in fig. 3(a)). Its surface patch has two boundaries: an outer one (γ o ), directed towards feature points, and an inner one (γ i ). A node N is a junction node if deg(N) > 2 (in red in fig. 3(a)). We define a junction area as a connected set of junction nodes. Concretely speaking, junction areas correspond to surface patches adjacent to several large features, such as the palm of the hand on a humanoid model. In figure 3(a), the set of four red nodes forms a junction area. At this stage of the algorithm, junction nodes are merged into junction areas.

Core boundary definition

We define the core of an object as a connected set of nodes, including the root node of the skeleton and its possible adjacent junction areas. The root node is displayed in figure 3(a) with a bigger radius. It corresponds to the surface patch containing the vertex which is the furthest away from feature points. Concretely speaking, the core of the object corresponds to its global inner part, like the torso of a humanoid for example. The presence of tubular nodes on prominent components prevents an undesired expansion of the core. In figure 4(d), the core expansion is stopped at the basis of the head, the arms and the legs because constrictions have been detected on these features.

An additional mechanism is proposed to optimize core boundary computation. Constrictions are often identified at the basis of prominent components, delimiting thin patches. Consequently, to prevent over-segmentation in those areas (see the fingers in figure 3(b)), we propose to merge the first tubular node of a component with the core in the following cases: In figure 3, the first tubular node of the inch (white patch in fig. 3(b)) is large but its outer boundary (with the cyan patch) is deeply concave and better suits the purpose of segmentation as a delimitation between the inch and the palm. This tubular node is merged into the core. The first condition makes, for example, the first tubular node of the middle finger (pink patch in fig. 3(b)) merge with the core.

As a conclusion, in figure 3, junction nodes first merges to form a junction area. Then this junction area merges with the root node. Finally, core boundary is optimized merging first tubular nodes if necessary (see figure 3(c)).

Topology driven feature hierarchy

In the previous stage of the approach, a fine segmentation of the object has been obtained running a node merging algorithm on the skeleton (figs. 4(a), 5(e)). Based on this fine segmentation, we propose in this section a fine to coarse strategy for the semantic-oriented composition of features.

Traditional hierarchical segmentation algorithms [10, 9, 2] base their hierarchy scheme on low level measurement clustering. In this paper, we propose to drive the hierarchy composition using a high level shape description so as to provide a semantic-oriented feature hierarchy. The key idea of our hierarchical scheme is that the topology of a feature is a more important decomposition criterion than its geometry. Consequently, we propose to subdivide features in priority along junction areas, which corresponds to skeleton's topological variations.

First, each core-adjacent connected component of the skeleton is identified as a feature. This identification results in the coarsest level of hierarchy (figs. [START_REF] Berreti | Partitioning of 3D meshes using Reeb graphs[END_REF]

(b) and 5(b)).

Then features are recursively subdivided according to the following strategy. Each feature is swept from its inner boundary to its extremity (extremity node). If a junction area is encountered, the previously swept nodes are gathered into a new feature, and the junction area is gathered with the remaining nodes into another feature. If there is no junction area in the feature, it is subdivided into the nodes it is composed of and the recursive algorithm stops.

In figure 4(a), the baby's left arm is swept from its basis until the hand palm (junction area, in red) is reached. The two tubular nodes (in blue) are gathered into a feature (in white in fig. 4(c)) and the hand into another (in pink in fig. 4(c)). At this level of hierarchy, arms are divided into hands and arms. Then, these features are subdivided into the nodes they are composed of in the last hierarchical level.

In this section, we presented a semantic-oriented hierarchical mesh decomposition algorithm. Contrary to state of the art methods, it subdivides features according to their topological characteristics, providing meaningful decompositions of limbs into arms and hands or legs and feet (fig. 5(c)) for example, which benefits shape understanding.

Experiments and results

In this section, we present experimental results obtained on manifold, connected and closed triangulated surfaces extracted from the Princeton [START_REF] Shilane | The Princeton shape benchmark[END_REF] and the INRIA [1] shape repositories.

Segmentation evaluation

To our knowledge, no ground truth 3D mesh segmentation evaluation process has been proposed in the past. However, Attene et al. [START_REF] Attene | Mesh segmentation: A comparative study[END_REF] enumerate seven criteria for evaluating a segmentation algorithm. 1. Versatility: figure 8 shows that visually meaningful segmentations are obtained on primitive shapes (with genus 2), animals, anatomic and mechanical parts. 2. Decomposition semantic: in figures 4, 5, and 8, animals and humanoids are first segmented into core and limbs. Moreover, figure 6 shows a compatible segmentation of two horse models (into core and limbs, at first level). 3. Boundary location: a curvature constrained geodesic distance computation has been proposed in section 2 so as to align constrictions along concavities. This can be observed on models with frank or coarse concavities (figs. 8(a), 8(e)). However, when no constriction is identified (fig. 8(d)), the object is segmented with coarse boundaries.

4. Hierarchy semantic: on the contrary to related work, our hierarchy scheme is driven by high level shape descriptions. Consequently, it produces a more natural categorisation of features. In figure 5(c), dinopet's limbs are decomposed into arms and hands or legs and feet while those limbs are simply subdivided in the middle in [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF] at the same level of hierarchy. Moreover, figure 6 shows that this hierarchy is compatible when dealing with objects that are visually similar. Horses are first decomposed into core and limbs, then limbs are subdivided and the head is decomposed at the last level of hierarchy in both cases. 5. Robustness: In fig. 7(b), the number of vertices has been reduced by 5. In fig. 7(c), each vertex has been moved randomly into a volume equal to 1% of that of the object, reproducing scanner noise for example. In fig. 7(d), the object has been deformed via user-interaction as described in [START_REF] Tierny | 3D mesh skeleton extraction using topological and geometrical analyses[END_REF]. After these degradations, the number of extracted patches N P is quite stable, even with strong surface noise. Moreover, the deformed object is segmented identically to the original one, which underlines the pose robustness of the algorithm. Model with N J the number of junction areas in the feature. Table 1 shows running times obtained on a 3GHz P4 PC. For information, the dinopet is segmented in 28 s. in [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF] while our method takes 3.9 s. 7. Control parameters: some low level parameters are used in our method (α 0 , ∆ f m and ζ 0 ) but they have been fixed experimentally for all models. Therefore, our algorithm is fully automatic.

As a comparison to related work, feature points are accurately extracted, even on small features such as the horse's ears, where segments are isolated in figure 9 (black and white patches), on the contrary to [START_REF] Zhang | Feature-based surface parametrization and texture mapping[END_REF][START_REF] Lien | Simultaneous shape decomposition and skeletonization[END_REF]. Also notice that limbs are more accurately subdivided with our method than in [START_REF] Zhang | Feature-based surface parametrization and texture mapping[END_REF] (where the front right leg, in red, is not subdivided). Moreover, our hierarchical scheme decomposes features in a more semantic-oriented manner than [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF] (fig. 5). Finally, running times are significantly lower than clustering based methods [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF][START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF]. On the contrary to [START_REF] Lien | Simultaneous shape decomposition and skeletonization[END_REF], the genus of the surface does not impact the computation time (cf. table 1).

Limitations

Because it is topology based, our method cannot distinguish features when they form a compact connected com-(a) Zhang [START_REF] Zhang | Feature-based surface parametrization and texture mapping[END_REF].

(b) Lien [START_REF] Lien | Simultaneous shape decomposition and skeletonization[END_REF].

(c) Our algorithm. ponent. For example, a closed fist will not be decomposed if the fingers are stuck to each other. This drawback can be observed in figure 3(d) where fingers are stuck at their basis, preventing the extraction of the first phalanx as an individual feature. The same effect can be observed in figure 5, where the upper part of the dinopet's legs cannot be distinguished from the core. To overcome this issue, extending extracted boundaries with an optimization algorithm could be an interesting direction for future work. Moreover, local re-meshing techniques could improve boundary locations on coarsely designed objects, such as the bi-torus (fig. 8(d)).

Conclusion

In this paper, we presented a fully automatic topology driven 3D mesh hierarchical segmentation algorithm. Feature boundaries and feature hierarchy are both computed in a semantic-oriented manner, subdividing character limbs into arms and hands for example. Method performance has been evaluated according to latest evaluation methodology [START_REF] Attene | Mesh segmentation: A comparative study[END_REF]. Experiments shown the rapidity of our algorithm as well as its robustness against various surface degradations.

We made the following contributions. We adapted shape topological description techniques to the segmentation problem, using a curvature constrained geodesic based mapping function and identifying constrictions along resulting level lines. We also proposed a region merging strategy that produces a fine segmentation, based on the previously computed topological decomposition. Finally, we introduced a new semantic-oriented hierarchical scheme, where features are subdivided according to their topological complexity. In the future, thanks to the progressive shape understanding provided by the algorithm, we would like to address the partial shape retrieval problem [START_REF] Funkhouser | Modeling by example[END_REF].

Figure 1 .

 1 Figure 1. Enhanced topological skeleton extraction overview [18].

Figure 2 .

 2 Figure 2. Original geodesic based function (a). Curvature index [11] (b). Curvature constrained geodesic based function (c). In (d), level sets are aligned with concave areas.

  (b), each patch has been displayed with(

Figure 3 .

 3 Figure 3. Core and junction computation results in the finest segmentation of the object.

Figure 4 .

 4 Figure 4. Hierarchical segmentation of a baby model.1. Related surface patch is small: its interval on f m is lower than a fixed threshold ∆ f m (∆ f m = 0.1); 2. Its outer boundary (γ o ) is a better candidate for core boundary definition than its inner boundary (γ i ):ζ(γ o ) < ζ 0 and ζ(γ o ) < ζ(γ i ) (ζ 0 = -0.2).In figure3, the first tubular node of the inch (white patch in fig.3(b)) is large but its outer boundary (with the cyan patch) is deeply concave and better suits the purpose of segmentation as a delimitation between the inch and the palm. This tubular node is merged into the core. The first condition makes, for example, the first tubular node of the middle finger (pink patch in fig.3(b)) merge with the core.As a conclusion, in figure3, junction nodes first merges to form a junction area. Then this junction area merges with the root node. Finally, core boundary is optimized merging first tubular nodes if necessary (see figure3(c)).

Figure 5 .

 5 Figure 5. Comparison to the hierarchical decomposition from [9]. At the second level, our algorithm segments dinopet's limbs into arms and hands or legs and feet.

(a) 1

 1 st level. (b) 2 nd level. (c) 3 rd level. (d) 1 st level. (e) 2 nd level. (f) 3 rd level.

Figure 6 .

 6 Figure 6. Compatible hierarchical segmentation on two horses. The models are identically segmented at each step while the second one has an additional feature (tail).

Figure 7 .

 7 Figure 7. Algorithm robustness against various surface degradations.

6 .Figure 8 .

 68 Figure 8. Finest level of hierarchy for various objects.

Figure 9 .

 9 Figure 9. Comparison to other skeleton driven segmentation algorithms [20, 13].

Table 1 . Computation times and number of surface patches at the coarsest and the finest levels of hierarchy.

 1 Faces 1 st level Last level Time (s.)

	Hand	52 000	7	12	107
	Baby	10 000	6	24	5.4
	Chess piece	600	3	5	0.05
	Human jaw	2 300	3	7	0.43
	Dolphin	4 200	6	8	0.89
	Bi-torus	6 000	5	7	1.3
	Bird	2 000	6	14	0.4
	Fan	450	6	8	0.06
	Mechanical	2800	3	5	0.53
	Horse 1	40 000	6	21	34
	Dinopet	9 000	7	28	3.9
	Horse 2	5 000	7	20	1.35

ton. Feature hierarchy composition takes O((N J + 1) × N)

Acknowledgements

This work is partially supported by the European Network of Excellence Delos No.

507618http://www.delos.info.


Nuno Gaspar
email: nuno.gaspar@inria.fr

Eric Madelaine
email: eric.madelaine@inria.fr

P Collet

P Merle

Fractal à la Coq

Component-based Engineering aims at providing a modular means to specify a wide range of applications. The idea is to promote a clean separation of concerns, and thus reusability, in order to ease the burden of software development and maintenance. The specification of such component models however, tends to be informal, leaving their inherent ambiguities open to interpretation.

In this paper we present our ongoing work towards a formal specification of the Fractal Component Model mechanized in the Coq Proof Assistant. An operational semantics for building component-based architectures is presented, along with its compliance with the Fractal specification.

Introduction

Building software is becoming a more and more complex task. Our era of computers and electronics demands highly sophisticated applications, where not only performance but also reliability are expected. Coping with such needs requires special care when designing software. The approach promoted by Component-based Engineering tackles these requests by leveraging the view of software as building blocks that when put together form the intended functionality. This is the rationale followed by the Fractal Component Model [START_REF] Bruneton | The fractal component model[END_REF].

Essentially, the essence of the Fractal Component Model lies around the notions of component, interface and binding. The components represent entities, generally pieces of software code, that communicate explicitly and necessarily by means of interfaces. This is in fact the main difference between Component-based programming and standard Object-oriented programming. The inherent advantage is that the dependencies between software blocks become explicit, permitting an easier software maintenance. At last, bindings convey procedures calls or messages passing between interfaces.

Components are deemed composite if they possess subcomponents. These can be shared, thus possibly yielding graph architectures and not only tree structures. Moreover, monitoring and on the fly reconfiguration of applications are other concerns tackled by Fractal.

Motivation and Related Work

Many approaches regarding the formalization of component models can be found in the literature. While it is still an utopia to see formal methods as part of the standard software engineering discipline, the maturity attained by formal tools already permit considerable achievements.

In [START_REF] Henrio | A framework for reasoning on component composition[END_REF], Henrio et al. present a framework for reasoning on component composition mechanized in Isabelle/HOL. Their focus is on a GCM-like component model 1 and show promising results for dealing with reconfiguration properties in a mechanical way.

Moreover, a formal specification of the Fractal Component Model has been proposed in the Alloy specification language [START_REF] Merle | A formal specification of the Fractal component model in Alloy[END_REF]. This work proves the consistency of a (set-theoretic) model of Fractal applications. This is perhaps the work most closely related with ours.

Their specification however, is constrained by the first-order relational logic nature of the Alloy Analyzer. In fact, they point to the use of the Coq Proof Assistant in order to overcome this limitation.

In this paper we present our work in progress towards a mechanization of the Fractal Specification in the Coq Proof Assistant. Our approach differs from the remaining in that we aim at providing a semantics that allows the building of correct-by-construction Fractal architectures.

The remainder of this paper is organized as follows. Section 2.1 shows the encoding of Fractal's basic definitions. Next, section 2.2 details a semantics for building Fractal Architectures. The kind of properties that can be proved in our development are discussed in section 2.3. Finally, section 3 concludes by pointing some directions for future work.

Mechanizing Fractal Architectures

As mentioned above, the Fractal Component Model has three core elements: components, interfaces and bindings. As such, we shall consider a simple language built around these core elements. In the following, we omit the obvious definitions for the sake of space.

Basic Definitions

Let us first demonstrate how interfaces and components are mechanized in the Coq Proof Assistant. An interface is characterized by an identifier, a type, a path identifying the component where the interface belongs to, whether it is accessible internally or externally, whether is it supposed to communicate as client or server, whether it serves functional or control purposes, and finally the language it implements. The intended meaning of each of these fields should be clear. The only complex aspect may arise from the path field. Fractal is a hierarchical component model, and since by introspection an interface is able to identify the component it belongs to, a path identifying this component is necessary. Its definition is a list of identifiers, where these identifiers indicate the components that need to be traversed in the hierarchy to reach the component holding the interface. A binding is always established between a client and a server interfaces. It can be a normal, export or import binding. These are defined by indication of the component's path which they belong to, and by the identifiers of the involved components and interfaces. For export and import bindings the identifier of one intervening component can always be inferred and thus is omitted.

Operational Semantics of Fractal Architectures

Having defined our core elements, it is now time to introduce a notion of state. A state is the structure that holds all the relevant information. For instance, for a tiny programming language it is generally a function mapping variables to their values.

1 D e f i n i t i o n s t a t e := component . In our case, a state has the same shape of a component. An empty state is therefore a component named Root without any subcomponents, interfaces and bindings.

The syntactic categories of our language for building Fractal Architectures are defined as follows.

a ::= mk component component | mk interface interface | bind binding | a; a | done
The meaning of the first three constructs should raise no doubt: making components, making interfaces and establishing bindings. We use a for representing actions2 . Allowing for multiple actions to execute as a sequence is defined by means of the standard operator ;. Last, done stands for the completed action.

The design of software architectures can be seen from a transition system point of view. One makes some action a, in some state σ, and ends up with a reduced action a' in some state σ'. This can be represented by the following manner. a, σ -→ a , σ .

Building Fractal architectures will therefore require to define these transition rules for each constructor of our language. In other words, to define a semantics.

In the following we use the . notation for projections. The above rule dictates that for creating a component one must supply a valid path, it must be itself well formed, and its identifier must be distinct from the ones in the same scope. The valid component path predicate expresses the well-formedness of a path in the hierarchy of components for the current state. Component well-formedness is inductively defined in Coq as follows. As for interfaces, they are deemed well formed provided that their pair of values ident × accessibility is unique. Moreover, each binding must be valid, that is, be of normal, import or export kind. These predicates are defined analogously to the precedent. At last, it is requested that inserting a component will not break the unicity of identifiers in the hierarchical level being inserted to. The function get scope : path → state → list component returns the list of components where the insertion is supposed to occur. As in the above rules, a valid path is expected. Further, establishing a binding also requires the valid component binding predicate to hold. As stated in the Fractal specification, a binding is either a normal, export or import binding. For instance, a binding is deemed normal if it is established between two external interfaces, and the components owning these interfaces have same enclosing component, i.e., they are in the same composite component. Moreover, a binding must be established between client to server interfaces.

Making Components

Making Interfaces

Sequence of Actions Finally, our last action serves the purpose of composing actions. Its definition is straightforward: having a composition of an action a 1 with an action a 2 , we first need to reduce a 1 so we can reach a 2 .

a 1 , σ -→ a 1 , σ a 1 ; a 2 , σ -→ a 1 ; a 2 , σ (SSeq 1 ) a 1 , σ -→ done, σ a 1 ; a 2 , σ -→ a 2 , σ (SSeq 2 )
So far, all the above rules defined a one step transition. However, there are cases where we may want to reason about multiple steps transitions. For instance, the action M k i n t e r f a c e ( ( I d e n t "A" ) : : n i l ) ( I d e n t "X" ) SomeType E x t e r n a l C l i e n t C o n t r o l C ; requires five steps to complete. Here, what we would like to do is to be able to reason on the overall execution of this action. As such, we need to formalize this notion of multiple steps. This is achieved by the reflexive transitive closure of our transition relation. Considering step : action * state → action * state → P rop, our inductive definition modelling our transition rules, we can define the notion of multiple steps in the following manner.

1 N o t a t i o n " a ' / ' s '---> * ' a ' ' / ' s ' " := ( r e f l s t e p c l o s u r e s t e p ( a , s ) ( a ' , s ' ) ) . Having this notion defined we can now reason on the overall execution of a sequence of actions. For instance, the above trivial lemma states that the execution of our SimpleArchitecture will indeed terminate.

Meeting the Specification

The rigorous definition of the possible actions that one can take allow us to formally reason on all possible outcomes. One important theorem is to prove that starting an architecture in a well formed state, when the execution completes it will end up in a state that is also well formed. It should be noted that the above theorem is an universal quantification over all possible Fractal architectures. Thus, every architecture built with our semantics rules and obeying the wellformedness hypothesis will be well formed.

Preliminary results allow us to conclude on the feasibility of this approach. At this stage however, the proof is not completed and a more precise definition of our well formed predicate is needed for a full compliance with the Fractal specification.

Final Remarks

In this paper we presented the first steps towards a mechanized specification of the Fractal Component Model. At this stage, there is still some remaining work to be included in our Coq development. Inclusion of a progression theorem and typing rules on bindings are natural next steps. Furthermore, a seamless integration with the Vercors platform [START_REF] Rabéa | Verifying safety of faulttolerant distributed components[END_REF] could be envisaged. In the long term one of our goals is to use this formalization in conjunction with model checking techniques in order to address behavioural properties of complex component architectures.

Moreover, our action language opens interesting perspectives w.r.t. reconfiguration capabilities of component-based systems. Another research direction to follow would be to see how could we best leverage this kind of approach with the concerns of autonomic reconfigurations.

1 2 | 3 t y p e -> ( * i t s t y p e * ) 4 path -> ( * path t o t h e component i t b e l o n g s * ) 5 a

 2345 I n d u c t i v e i n t e r f a c e : Type := I n t e r f a c e : i d e n t -> ( * i t s i d * ) c c e s s i b i l i t y -> ( * I n t e r n a l / E x t e r n a l * ) 6 communication -> ( * c l i e n t / s e r v e r * ) 7 f u n c t i o n a l i t y -> ( * f u n c t i o n a l / C o n t r o l * ) 8 l a n g u a g e -> ( * t h e l a n g u a g e i t implements * ) 9 i n t e r f a c e .

1 2 | 3 t y p e -> ( * i t s t y p e * ) 4 path -> ( * i t s path on t h e h i e r a r c h y * ) 5 c 2 | 3 | 4 |

 2345234 I n d u c t i v e component : Type := Component : i d e n t -> ( * i t s u n i q u e i d * ) o n t r o l L e v e l -> ( * Lowest / I n t r o s p e c t i o n / C o n f i g u r a t i o n * ) 6 l i s t component -> ( * i t s sub components * ) 7 l i s t i n t e r f a c e -> ( * i t s i n t e r f a c e s * ) 8 l i s t b i n d i n g -> ( * i t s b i n d i n g s * ) 9 component . 2 Actes de CIEL 2012 A component has an identifier, a type, a path indicating its level in the hierarchy, a control level determining what it can do, subcomponents, interfaces and bindings. Bindings are connecting components together by means of their interfaces. 1 I n d u c t i v e b i n d i n g : Type := Normal : path -> i d e n t -> i d e n t -> i d e n t -> i d e n t -> b i n d i n g Export : path -> i d e n t -> i d e n t -> i d e n t -> b i n d i n g Import : path -> i d e n t -> i d e n t -> i d e n t -> b i n d i n g .

2 3

 2 D e f i n i t i o n e m p t y s t a t e : s t a t e := 4 Component ( I d e n t " Root " ) Top n i l C o n f i g u r a t i o n n i l n i l n i l .

c=

  Component id t p cl subComps interfaces bindings valid component path p σ well formed component c ∀ c', c' ∈ (get scope p σ) → (c'.id = id) make component c, σ -→ done, add component σ c (SM akeComponent)

1 2 | 4 w|

 24 I n d u c t i v e w e l l f o r m e d c o m p o n e n t ( c : component ) : Prop := WellFormedComponent Base : 3 w e l l f o r m e d i n t e r f a c e s ( c->i n t e r f a c e s ) -> e l l f o r m e d b i n d i n g s ( c->b i n d i n g s ) ( c->subComponents ) ( c->i n t e r f a c e s ) -> 5 ( c->subComponents ) = n i l -> 6 w e l l f o r m e d c o m p o n e n t c 7 WellFormedComponent Step : 8 ( f o r a l l c ' , I n c ' ( c->subComponents ) -> w e l l f o r m e d c o m p o n e n t c ' ) l f o r m e d i n t e r f a c e s ( c->i n t e r f a c e s ) -> 11 w e l l f o r m e d b i n d i n g s ( c->b i n d i n g s ) ( c->subComponents ) ( c->i n t e r f a c e s ) -> 12 w e l l f o r m e d c o m p o n e n t c .

  i = Interface id t p a c f l valid interface path p σ c = get component with path p s ∀ i', i' ∈ (c.interfaces) → i'.id =id → i'.accessibility = a make interface i, σ -→ done, add interface σ i (SM akeInterf ace) As in the creation of components, creating an interface also requires a valid path. The function get component with path : path → state → component returns the component where the interface is supposed to be added. Inserting the interface must not break the property of unique pair of values ident × accessibility. Making Bindings valid component path (b.path) σ valid component binding b components make binding b, σ -→ done, add binding σ b (SM akeBinding)

1 2 | 3 | 4 |

 234 I n d u c t i v e v a l i d c o m p o n e n t b i n d i n g ( b : b i n d i n g ) ( s : s t a t e ) : Prop := NormalBinding : n o r m a l b i n d i n g b s -> v a l i d b i n d i n g b s E x p o r t B i n d i n g : e x p o r t b i n d i n g b s -> v a l i d b i n d i n g b s I m p o r t B i n d i n g : i m p o r t b i n d i n g b s -> v a l i d b i n d i n g b s .

1

  D e f i n i t i o n S i m p l e A r c h i t e c t u r e : a c t i o n := 2 Mk component n i l ( I d e n t "A" ) SomeType C o n f i g u r a t i o n n i l n i l ; 3 Mk component n i l ( I d e n t "B" ) SomeType C o n f i g u r a t i o n n i l n i l ;

4

 4

5 M

 5 k i n t e r f a c e ( ( I d e n t "B" ) : : n i l ) ( I d e n t "Y" ) SomeType E x t e r n a l S e r v e r C o n t r o l C ; 6 Bind ( ( I d e n t "A" ) : : n i l ) ( I d e n t "X" ) E x t e r n a l 7 ( ( I d e n t "B" ) : : n i l ) ( I d e n t "Y" ) E x t e r n a l .

2 3

 2 Lemma w e l l f o r m e d a r c h i t e c t u r e :

4 e x i s t s s , 5 S

 45 i m p l e A r c h i t e c t u r e / e m p t y s t a t e ---> * Done / s . 6 P r o o f .

1 2 match s with 3 |8

 23 D e f i n i t i o n w e l l f o r m e d ( s : s t a t e ) : Prop := Component i d t p c l l c l i l b => 4 i d = ( I d e n t " Root " ) /\ t = Top /\ p = n i l /\ c l = C o n f i g u r a t i o n /\ 5 l i = n i l /\ w e l l f o r m e d c o m p o n e n t ( Component i d t p c l l c l i l b ) Theorem v a l i d i t y : 9 f o r a l l a s s ' , a / s ---> * Done / s ' l f o r m e d s ' .

In the realm of programming languages semantics, instruction or statement are the terms usually employed. In our case we consider that action is a more adequate nomenclature.

Actes de CIEL 2012

Actes de CIEL 2012

Acknowledgements

The authors are grateful to Ludovic Henrio for numerous remarks in earlier versions of this paper and engaging discussions about possible future research directions.

The author is partially funded under contract CIFRE 2012/0109 ActiveEon/INRIA


Yanwen Chen

Yixiang Chen

Eric Madelaine

Timed-pNets: A formal communication behavior model for real-time CPS system

We propose a semantic model named timed-pNets to define hierarchical structures for CPSs as well as its communication behaviors semantic with time constraints. Logical clocks relations are introduced to describe the partial order of event occurring. After setting (time)boundaries and designing properties, we use the TimeSquare tool to simulate the system and check its properties.

Introduction

Nowadays, cyber physical systems have received much attention since the next generation of computing revolution is integration of computation, control and communication [START_REF] Lee | Cyber physical systems: Design challenges[END_REF]. In CPSs, heterogenous embedded devices are connected by wire or wireless networks to communicate among each other via sets of sensors and actuators. One typical application domain is intelligent transportation systems (ITS), in which cars and infrastructures are equipped with sensors and actuators. They communicate with each other to update physical information and accomplish remote controlling. Currently, Research and Innovative Technology Administration (RITA) [START_REF]Intelligent transportation systems[END_REF] in U.S. Department of Transportation has started research work on it to achieve a vision of national transportation by feature a connected transportation environment among vehicles, infrastructures and passengers' portable devices. They raised the importance of real-time communication among these distributed nodes since the data out of date would make big mistakes even sometime could lead to car accidence. For example, the delay of sending global traffic information to cars may result to a wrong guiding for car to choose its best way. Also the delay of information exchanging among cars may cause a car accidence especially when they cannot see each other at cross.

Our aim is to build a low-level semantic model for the system and then analyse its properties. Since each distributed device has its own clock(s) and the communication delay between devices is uncertain, synchronous and asynchronous communication behavior will be considered in our model. We propose timed-pNets, which is an extension of pNets [START_REF] Barros | Behavioural models for distributed fractal components[END_REF], to describe the communication behavior by adding logic clock relation. Time-pNets absorbs many advantages of pNets such as well hierarchical structure, flexible communication models, compact format, expressiveness, etc. By introducing clock relations which is defined in CCSL [START_REF] Mallet | CCSL: specifying clock constraints with UML/MARTE[END_REF], the asynchronous communication behavior can be well defined which leads to well checking for real-time properties. Finaly, we use the TimeSquare tool to simulate the clock relation of models and then check the correction of time logic as well as properties.

The main contributions of our work in the paper include: 1) Introduce logical clocks into pNets, 2) Propose a model language timed-pNets , 3) Use timed-pNets to model ITS, 4) Use TimeSquare to simulate clock relations and check its properties.

In section 2 we introduce the related work. In section 3 we give the definition of timed-pNets. In section 4, we propose a simple use case and describe how we use timed-pNets to model it. In section 5, we use time square to simulate the clock relations of the system and check its properties. In the last section, we give an conclusion for our current work.

Related Work

Prior research works related to building model for ITS applications include:

-Timed-automata [START_REF] Alur | A theory of timed automata[END_REF] can be used to model the behavior of real-time systems.

They provide a simple, and yet powerful, way to annotate state-transition graphs with timing constraints using finitely many real-valued clocks. Closure properties, decision problems as well as automatic verification of realtime requirements of finite-state systems are considered in timed-automata, and are supported by a number of tools, e.g. in UPPAAL [START_REF] Bengtsson | Uppaal -a Tool Suite for Automatic Verification of Real-Time Systems[END_REF]. -Globally Asynchronous Locally Synchronous (GALS) systems [START_REF] Muttersbach | Practical design of globallyasynchronous locally-synchronous systems[END_REF] combine the benefits of synchronous and asynchronous systems, in which each embedded node is modeled as a FSM (Finite State Machine) and the communication between them as buffers. This architecture provides a methodology for combining concurrent embedded systems within loosely coupled systems. -The BIP framework [START_REF] Basu | Modeling heterogeneous real-time components in BIP[END_REF], which can be used to model heterogeneous real-time components. BIP provides a powerful mechanism for structuring interaction for layered components of which synchronous and time systems are particular classes. The BIP framework produces a very fine grained formal computational model of the system functional level and executes the model semantics at runtine. In parallel, recent evolutions of the BIP framework allow the use of timed models and provide a real-time BIP Engine. -Spatio-temporal consistence language(STeC) [START_REF] Chen | A location-triggered specification language for real-time systems[END_REF], which provides a locationtriggered specification as well as operational semantics for describing distributed system with time and location constraints.

Compared to these previous works, our appraoch uses logical clocks instead of physical clocks, so that our approach is flexible enough to describe a communication delay by clock relation rather than by concrete time units. This gives us a very flexible way to specify interaction of systems with different clocks.

Timed-pNets

In this section, we define Timed-pNets, which are an extension of pNets (parameterized networks of synchronized automata), a very expressive and flexible semantic model developped by the Oasis team at INRIA for the modeling and verification of (untimed) distributed systems [START_REF] Barros | Behavioural models for distributed fractal components[END_REF].

Definition 1 (Timed-Actions). Let T be a set of discrete time variables taken from non-negative natural numbers N. B T is the set of boolean expressions (guards) over time variables and L A,T is an action set built over T , in which each action has a free time variable t ∈ T . We call a t|B ∈ L A,T , with B ∈ B T a timed action, in which t describes a time delay before the action can be executed.

We set a 0 = a that means the action a is always ready. As an example, a t|1≤t≤3 means the action a cannot be executed until t units times are passed.

Logical Clocks of timed-Actions Logical clocks [START_REF] Mallet | CCSL: specifying clock constraints with UML/MARTE[END_REF] represent a relaxed form of time where any events can be taken as a reference for counting. It can be used for specifying classical and multiform real-time requirements as well as formally specifying constraints on the behavior of a model.

Definition 2 (Logic Clocks).

A clock C a consists of an infinite set of discrete ticks of timed-action a t . We write

C a = {(a t ) 1, (a t ) 2, ..., (a t ) k, ...}(k ∈ N, a t ∈ L A,T ), in which (a t ) k denotes the k th instance of clock C a .
Clock Constraints Clock constraints are predicates built from binary relations between clock expressions. We take the syntax and semantics of clock relations from [START_REF] Mallet | CCSL: specifying clock constraints with UML/MARTE[END_REF], which is a language to express multi-clock time specifications by defining clock relations of time models for real-time systems. The clock relations include: ⊆ (subclock), (exclusion), = (coincidence), ≺ (strict precedence), ( precedence) and defined as:

a 1 ⊆ a 2 (a 1 is a subclock of a 2 ), which means each instance of a 1 must be coincident with an instant of a 2 . a 1 a 2 (a 1 excludes a 2 ), which means none of their instances coincide.

a 1 = a 2 (a 1 coincides with a 2 ), which means the action a 1 ticks if and only if the action a 2 ticks. a 1 ≺ a 2 (a 1 strictly precedes a 2 ), which means for every instant k (k ∈ N ), the k th instant of a 1 strictly precedes the k th instant of a 2 . a 1 a 2 (a 1 precedes a 2 ), which similar to the previous one. The only different is the action a 1 can tick as late as when a 2 ticks. -New relations or expressions can be define by combining these basic relations with arithmetic and boolean. For example:

a 1 -a 2 ≺ a 3
The next definition extends the classical pNets definition from [START_REF] Barros | Behavioural models for distributed fractal components[END_REF] with timedactions and clock constraints. From pNets, we retain the hierarchical structure that is essential in structuring our heterogeous systems, but also the parameterization of subnets: holes in a pNet can be instanciated by a variable number of subnets (as e.g. a number of cars in the forthcoming case-study). Then synchronisation vectors allow very flexible and expressive multi-way synchronisation mechanisms, that naturally we extend here with clock constraints. Definition 3 (Timed-pNets). A Timed-pNet is a tuple < P, A G , R G , J, C, O J , R J , -→ V >, where:

-P = {p i /p i ∈ Dom i } is a finite set of parameters -A G ⊆ L A,T
is a set of global actions -C is a set of clocks for all actions -R G is a set of relations between actions taken from each subnet -J is a countable set of argument indexes: each index j ∈ J is called a hole and is associated with a sort O j ⊆ L A,T and a set of clock constraints

R J - - → V = { - → v } is a
a tg g ∈ A G , k i1 ∈ Dom 1 , k i2 ∈ Dom 2 , !a t1 ∈ O i1 , ?a t2 ∈ O i2 , t g = max{t1, t2} * or (visibility) - → v =< ..., a t1 [k1] , ... >→ (a tg g ), in which a tg g ∈ A G , k1 ∈ Dom 1 , a t1 ∈ O i1 , t g = t1

Case Study

In this part we illustrate our approach with a simple use case called speed controlling system taken from [START_REF]Intelligent transportation systems[END_REF]. Here cars' speed are monitored by infrastructure that collects information from cars and sends brake signal back to cars under some global decision procedure. The communication protocol is described as following:

-Cars send heartbeat signals "I'm here" with parameters "(location, speed)".

-Infrastructure collects heartbeat signals from cars.

-Infrastructure sends "brake" signal to the car if it is over the speed limitation.

-A car reduces its speed when it gets the "brake" signal. Fig. 1 presents its architecture in which cars and infrastructures are distributed nodes. A cars consists of three sub components: a sensor, a controller and a brake system. The car sensor is used to detect its current location and speed and to receive control signals received from the infrastructure. The Car controller gets signals from the sensor and then call the relevant systems to execute brake operations. The local communications between the sub components of cars are synchronized, which means that the sending event and receiving event coincide. These sub components' LTS are shown in Fig. 1. The car sensor is modeled by two LTSs: one defining the periodical sending of heartbeat signals to report its location and speed, the other describing its reactions to control signals.

Now, we describe how to use timed-pNets to model the communication behavior and how to build clock relations in and between components. By lack of space we only explain the top-level synchronisation elements, and the behaviour of the car's subprocesses.

Formalisation of timed-pNets Architecture

The building of timed-pNets is hierarchical. For our exemple, we generate separately a pNet structure for the toplevel assembly, and for each of the Car and Infrastructure components. The top-level pNet has 2 holes, the first one receiving an arbitrary number of Cars, the second one a single Infrastructure. Within each "second layer" pNet local communications will be synchronous, while at toplevel communications are asynchronous. This shows in the following pNet, where clock relations on the "hb" and "ctrl" links are defined as precedence. where [START_REF] Alur | A theory of timed automata[END_REF] describes that the heartbeat signal is sent periodically; (2) indicates that events τ and !ctrl(brake), from different paths of the same LTS, are exclusive; (3) denotes that the event ?c ctrl(result) always precedes the event τ and !ctrl(brake); (4) tells us that event !ctrl(brake) precedes the event ?T s; [START_REF] Chen | A location-triggered specification language for real-time systems[END_REF] explains that the events in the i th cycle precedes those in the (i + 1) th cycle.

< P = {k : N}, A G , R G , J, C J , O J , R J , - → V > A G = {CI hb t h (k, loc, speed), CI ctrl tct (k, brake)} J = {car[k], inf rastructure} O Car = {!c hb t hb c (loc, speed), ?c ctrl tct c (brake), !call(brake), T s , ...} O Inf rastructure = {?I hb t hb I (k, loc, speed), !I ctrl t ct I (k, brake), ...} R G = {!c hb t hb c [k](loc, speed) ≺?I hb t hb I [k](loc, speed); ............!I ctrl t ct I (k, brake) ≺?c ctrl tct c [k](brake);} - → V :< O car [k], O inf

Simulation

We use TimeSquare [START_REF] Deantoni | TimeSquare: Treat your Models with Logical Time[END_REF] to simulate the clock relations and check its logic correction. The input of TimeSquare is a CCSL file including clock relations, bound requirements and properties. The tool proceeds with a symbolic simulation, and generates a trace model (one partial order satisfying the specification). Output files (text and graph) are generated to display the traces and eventually show the property violations.

In our use-case the clock relations generated from the pNets model. Then we representing the communication and computation delays by fixing delay bound in timed-action guards. We set heartbeat interval hi = hb (i + 1) -hb i. The minimum and maximum communication delay between car and infrastructure is set as (1/5)hi and (3/5)hi. For the computation delay, we assume that it takes at most (2/5)hi for each action transition so that for instance we can set (?c ctrl(result)-?T s)

(2/5)hi. For deadline, assuming each heartbeat signal should be processed before sending next heartbeat, then for each action a i, we set deadline = hb (i + 1). All these boundaries are merged into the TimeSquare CCSL input file. We expressed a "boundary liveness property" for this simulation to see if the system satisfies the real-time property under the hypothesis of boundaries. (?T s (i) ≺ hb (i + 1)) ∧ (τ i ≺ hb (i + 1)) denotes each heartbeat signal finally will be processed before deadline.

After simulation, we got the result as Fig. 2, which shows that the real-time property is not satisfied since the clock !ctrl(brake) is over the deadline. The boundary condition we set in previous part was too large. After we modify the maximum computation boundary, the simulation shows no more violations.

Conclusion

In this paper, we have added time constraints to the pNets behavior semantic model. In timed-pNets, logical time relations in lower-level (synchronous) components are derived from the corresponding label transition systems. Then Fig. 2. property checking communication delays, processing delays and required global properties are defined by the user. We illustrate our approach on a simple use case from Intelligent Transport Systems, show how our Timed-pNets models are constructed, and how timed properties are validated through simulations in the TimeSquare tool.

Fig. 1 .

 1 Fig. 1. Timed-pNets architecture with details of the car's subprocesses

  rastrcuture >→ A Car inf rastrcuture ... =<!c hb t hb c [k](loc, speed), ?I hb t hb I [k](loc, speed) >→ CI hb t h (k, loc, speed); ...... <?c ctrl tct c [k](brake), !I ctrl t ct I (k, brake) >→ CI ctrl tct (k, brake).} An interesting point is that the Infrastructure receives independent heartbeats for the Cars, that are subsequently interleaved within the Infrastructure structure. This is expressed by a clock relation on the link between Infrastructure sensors and control: ?I hb t hb I [k](loc, speed) ⊆!hb t I I (k, loc, speed), telling the each Car heartbeat clock transmitted by the Sensor is a subset of the (single) heartbeat clock received by the Control. Finally, we give an example of the set of Clock constraints that we generate from a pLTS, here for the Car Sensor sub component: R CarSensor = {hb(loc, speed) idealClockdiscretizedByrate (1); ........................(τ !ctrl(brake)) (2); ........................?c ctrl(result) ≺ (τ ∧!ctrl(brake)) (3); ........................!ctrl(brake) ≺?T s (4); ........................(?T s[i] ∨ τ [i]) ≺?c ctrl(result)[i + 1] (5);}

  set of synchronous vectors of the form:

* (binary communication) -→ v =< ..., !a t1 [ki1] , ..., ?a t2 [ki2] , ... >→ (a tg g ), in which

This work was partialy funded by the Associated Team DAESD between IN-RIA and ECNU, Shanghai; by ECNU for overside visiting funding, China 973 project(No.2011CB302802) and NSFC(No.610211004).


Alexandre Janon
email: alexandre.janon@imag.fr

Reduced-basis solutions of affinely-parametrized linear partial differential equations

Introduction and overview

In some contexts (for example, in Monte-Carlo, or quasi-Monte-Carlo sensitivity analysis), it is necessary to compute the solution u µ of some partial differential equation (PDE) depending on an input parameter vector µ (dependence can be obviously through the PDE coefficients, and also through the domain on which the equation is posed, or through initial/boundary data), for a large number (say, thousands) values of µ. When no analytical solution of the PDE is known, one has to use computationally-intensive numerical methods -such as finite elements or finite differences -to approximate the solution of the PDE. This leads to complex computer codes, which can take several hours to runfor a single value of µ. In this case, one easily sees that simply calling the numerical code for each required value of µ -one-by-one -is not practicable at all ; and that we need to use a computationally cheaper code giving a surrogate solution u µ (also called metamodel) from a value of µ. More specifically, we split our computation in two phases:

• an offline phase, where we make a reasonable number of calls to the expensive code in order to "learn" about the solutions of the PDE ;

• an online phase, where, given some value of µ, we use the information collected during the offline phase to produce a good approximation of the solution of the PDE for this specific value of µ, faster than if we had no information at all. This strategy is efficient if the fixed additional cost of data collection during the offline phase (done once for every value of µ) is dominated by a strong reduction in the marginal cost (i.e., per value of µ) when one replaces the original, expensive numerical code by the cheaper "online" code.

In this report, we describe a procedure, named reduced basis approximation (RB), specifying the offline and online phases. The key to RB approximation is to look for (during the online phase) a surrogate u µ solution of the PDE as a linear combination of functions from a linearly-independent family (called the reduced basis) that has to be wisely chosen during the offline phase. For the resulting problem to be well-posed, one has to somewhat weaken the PDE, using a Galerkin-like method. An advantage of using such a procedure is that one can quantify the loss of information made when replacing the expensive numerical solution u µ by u µ , using an a posteriori fastly-computed error bound ǫ µ so that ||u µu µ || ≤ ǫ µ for some functional norm || • ||.

The report is organised as follows: in a first part, we present the necessary assumptions on the PDE for RB to be applicable to ; for the sake of simplicity, we restrict ourselves to the case of coercive, linear and "affinely parametrized" (we will explain what this means) PDEs, but it should be noted there are extensions to RB which somewhat relax some of these assumptions ; in a second part, we describe the core of the RB approximation, that are the offline and online phases of computation of the surrogate u µ ; in a third part, we give a procedure to compute the error bound ǫ µ ; in a fourth part, we show two different ways (greedy and POD) to choose the reduced basis family ; in a fifth and final part, we address the question of using RB to approximate outputs of the solution, that are (scalar) functionals of u µ .

Notations and assumptions

Our main reference for sections 1 to 3 is [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF].

As said earlier, we denote by µ ∈ ⊂ R p our parameter vector.

Our unknown is a function u µ -depending on µ, besides the usual (e.g. time and space) variables -belonging to some function space X , satisfying a µparametrized linear partial differential equation, which can be written under a variational (a.k.a. weak) form:

a µ (u µ , v) = f µ (v) ∀v ∈ X (1)
where a µ is a µ-dependent (continuous) bilinear form on X , and f is a (continuous) linear form on X .

We further assume that our PDE is affinely parametrized, that is:

• f µ ≡ f is µ-independent;
• and a µ can be written as a sum of Q bilinear forms:

a µ (w, v) = Q q=1 Θ q (µ)a q (w, v)
where:

a q are µ-independent bilinear forms on X ; -Θ q are arbitrary functions on , capturing all dependence in µ ; since these functions are to be evaluated many times during the online phase, we shall require these functions to be fast-computed.

Reduced basis approximation 2.1 Opportunity for a RB approximation

The function space X is typically a Sobolev space; such spaces, due to their infinite-dimensional nature, are not amenable to numerical computations. In order to discretize equation ( 1) and get an approximation of the actual mathematical solution, one can set X to some finite-element space (e.g., 1 ) and use

(1) to get a linear system of equations satisfied by the coefficients of u µ in a basis of X . In the sequel, X will always stand for such a space. This means that the reduced basis solution we are to compute is an approximation to the discrete numerical solution, and thus we rely on (and assume) the fact that the latter is a good approximation to the analytical solution. A good treatment of finite element discretizations can be found in [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF].

The computational cost required to compute this approximation is an increasing function of the dimension of the chosen discrete function space, and unfortunately, it has to be very large for the approximation to be accurate. One explanation for this fact is that those spaces are "too general", that is, they contain too many different functions and so we need to deal with a large number of degrees of freedom for our unknown.

On the other hand, we expect that two solutions of (1), for two different values of µ, will exhibit common "features". We could "encode" those features into a space X , spanned by some (µ-independent) functions ζ 1 , ζ 2 , . . . , ζ N ∈ X , and replace X by X in [START_REF] Bui-Thanh | Goaloriented, model-constrained optimization for reduction of large-scale systems[END_REF]. The family {ζ 1 , . . . , ζ N } is called the reduced basis. For a sufficiently regular-in-µ problem, we expect that taking N << dim X will still give -for a well-choosed reduced basis -an accurate enough solution, yet providing substantial computational economy.

Offline and online procedures

We suppose, for now, that we are given a reduced basis {ζ 1 , . . . , ζ N } of linearly independent functions in X (discussion about choice of this basis can be found in section 4). We set X = Span{ζ 1 , . . . , ζ N }.

Given µ ∈ , we look for u µ as a function belonging to X , thus :

u µ = N n=1 c n (µ)ζ n (2)
where the unknowns are c n (µ) ∈ R, n = 1, . . . , N .

These unknowns are to be found by imposing that u µ satisfy a relaxed version of (1); instead of verifying (1) for every v ∈ X , u µ satisfies (1) for every v ∈ X :

a µ ( u µ , v) = f (v) ∀v ∈ X (3)
Since X ⊂ X , (3) is weaker than (1). This standard technique is known as

Galerkin projection.

Substituting (2) into (3), and setting v = ζ p (p = 1, . . . , N ) gives that the c n (µ)'s satisfy the following system of N linear equations:

N n=1   Q q=1 Θ q (µ)a q (ζ n , ζ p )   c n (µ) = f (ζ p ) ∀p = 1, . . . , N (4)
This justifies the following offline/online procedure for the computation of u µ :

• in the offline phase:

we choose a reduced basis {ζ 1 , . . . , ζ N } (more details about this choice later)

we assemble, and store, Q matrices of size N × N :

a q = a q (ζ n , ζ p ) p,n=1,...,N
and the N -sized vector:

f = f (ζ p ) p=1...N
• in the online phase, being given a value of µ:

we assemble the N × N matrix:

a(µ) = Q q=1
Θ q (µ)a q -our vector of unknowns c(µ) = c n (µ) n=1...N is given by :

c(µ) = a(µ) -1 f
where the matrix inversion solves the system (4),

and so we get:

u µ = N n=1 c n (µ)ζ n
The complexity of the online phase is as follows: assembling of a(µ) takes O(QN 2 ) operations. Its inversion takes O(N 3 ) operations (one should note that a(µ) has no reason to be a sparse matrix -it is not, indeed). So the coefficients of u µ with respect to the reduced basis can be computed in less

than O(QN 2 + N 3 ) operations.
It is important to note that this complexity is independent of dim X . Thus online computation time is independent of the "quality" of the underlying reference numerical approximation, which can be as good as wanted without impacting marginal cost. We notice that we neglected the cost of evaluations of Θ q (µ), hence the hypothesis of fast-evaluation of Θ q we have made. We also notice that online complexity (in computational time and in storage requirements) is proportional to Q, the number of terms in the affine decomposition of a µ ; for our procedure to be efficient, Q should therefore be kept orders of magnitude smaller than dim X .

3 Error bound

Derivation of the error bound

In this section, we design an error bound ǫ µ so that:

||u µ -u µ || ≤ ǫ µ for all µ ∈ .
We should note once again that, thanks to our assumption about the quality of the underlying approximation, we measure only the error between the reference, expensive numerical solution u and our surrogate u ; error between u and the analytical solution is neglected.

We endow X with an Hilbert structure: an inner product < •, • > on X × X yielding a norm || • || on X . We need a crucial assumption about the bilinear form a µ describing our PDE: the α(µ) defined by:

α(µ) = inf w∈X sup v∈X a µ (w, v)
||w||||v|| has to be strictly positive: α(µ) > 0 ∀µ ∈ .

The quantity α(µ) is known as Babuška's "inf-sup" stability constant. The infsup stability hypothesis α(µ) > 0 is a sufficient condition for well-posedness of the underlying variational problem (1) (Babuška-Brezzi condition).

In the sequel, we suppose that a µ is symmetric. The inf-sup stability hypothesis is then equivalent to the coercivity of a µ (and the well-posedness of the variational problem comes from Lax-Milgram), and we have:

α(µ) = inf w∈X a µ (w, w) ||w|| 2
Let's move on to the error bound. First we remember that:

a µ (u µ , v) = f µ (v) ∀v ∈ X
Subtracting a µ ( u µ , v) from both sides of this equation, we get, thanks to bilinearity of a µ :

a µ (u µ -u µ , v) = f µ (v) -a µ ( u µ , v) ∀v ∈ X So: a µ (e µ , v) = r µ (v) ∀v ∈ X (5)
if we set:

e µ = u µ -u µ r µ (v) = f (v) -a µ ( u µ , v)
(e µ is known as the error, and r µ as the residual form).

Introducing α(µ), we have the following inequality:

a µ (e µ , e µ ) ≥ α(µ)||e µ || 2
We also have:

||r µ (v)|| ≤ ||r µ || X ′ ||v||
where

||r µ || X ′ = sup v∈X ||v||=1 r µ (v)
is the dual norm of the residual.

Thanks to these inequalities, and ( 5) with v = e µ , we get:

α(µ)||e µ || 2 ≤ ||r µ || X ′ ||e µ ||
and so, thanks to our stability assumption:

||e µ || = ||u µ -u µ || ≤ ||r µ || X ′ α(µ) (6)
So we have:

ǫ µ = ||r µ || X ′ α(µ)
In order to practically exploit the error bound [START_REF] Nocedal | Numerical optimization[END_REF], it remains to explain how to efficiently work out ||r µ || and α(µ). This question is addressed in the two following sections.

Dual norm of the residual ||r µ || X ′

Riesz representation theorem ensures that there exist y µ ∈ X so that

< y µ , v >= r µ (v) ∀v ∈ X (7)
Moreover, computing y µ is sufficient because of the additional property: 7) can be rewritten, thanks to the definition of r µ and the affine decomposition of a µ :

||r µ || X ′ = || y µ || Equation (
< y µ , v >= f (v) - Q q=1 Θ q (µ)a q ( u µ , v) ∀v ∈ X
and so, making use of the reduced basis decomposition of u µ (equation ( 2)):

< y µ , v >= f (v) - Q q=1 N n=1 c n (µ)a q (ζ n , v) ∀v ∈ X
This equation can be viewed as a linear variational problem whose unknown is y µ ∈ X . By Duhamel superposition principle, y µ can be written as:

y µ = γ - Q q=1 N n=1 Θ q (µ)c n (µ)Γ qn (8)
where γ ∈ X , and Γ qn ∈ X (q = 1, . . . , Q, n = 1, . . . , N ) satisfy:

< γ, v >= f (v) ∀v ∈ X (9) < Γ qn , v >= a q (ζ n , v) ∀v ∈ X (10)
Law of cosines applied to [START_REF] Volkwein | Proper orthogonal decomposition and singular value decomposition[END_REF] gives:

|| y µ || 2 = ||γ|| 2 -2 Q q=1 N n=1 Θ q (µ)c n (µ) < γ, Γ qn > + Q q=1 Q q ′ =1 N n=1 N n ′ =1 Θ q (µ)c n (µ)Θ q ′ (µ)c n ′ (µ) < Γ qn , Γ q ′ n ′ > (11)
Equation (11) will be our key to an offline/online procedure for computing

||r µ || X ′ .
Each of those phases should be done after the phases of computation of u(µ) described in section 2.2. The procedures read:

• in the offline phase:

we solve the QN + 1 linear variational problems (9) and (10) for γ and Γ qn (q = 1, . . . , Q, n = 1, . . . , N ) ;

we compute and store 1 + QN + (QN ) 2 inner products:

< γ, γ >, < γ, Γ qn >, < Γ qn , Γ q ′ n ′ > (q, q ′ = 1, . . . , Q, n, n ′ = 1, . . . , N )
• in the online phase, we compute ||r µ || X ′ = || y µ || using equation (11).

We note that the vectors γ and Γ qn computed during offline phase are not saved for the online phase, we store only 1 + QN + (QN ) 2 scalars. Thus in the online phase, the actual vector y µ is never formed. Doing so would make complexity of the online phase depend on dim X ; this is something we would like to avoid.

Stability constant α(µ)

We now turn to the problem about estimating α(µ).

Recall that

α(µ) = inf v∈X ||v||=1 a µ (v, v)
α(µ) is thus the optimal value of a minimization problem posed on X . The objective function in this problem is a quadratic form, we seek to minimize on the unit sphere. It is a well-known fact that α(µ) is the smallest eigenvalue of a (symmetric) matrix associated to the quadratic form a µ (•, •). Indeed, if for all v and w in X :

a µ (v, w) = v T M µ w and if < v, w >= v T Ωw for symmetric matrices
M µ and Ω, then α(µ) is the smallest λ ∈ R so that equation :

M µ v = λΩv admits a nonzero solution v.
Due to the high-dimensionality of X , it would be too expensive to solve this (generalized) eigenproblem online, for each value of µ. However, we expect that affine decomposition of a µ will provide us with an efficient offline/online procedure for estimating α(µ). By estimating here we mean finding a lower bound for α(µ), because the form of the error estimator [START_REF] Nocedal | Numerical optimization[END_REF] shows us we only need such a lower bound:

α(µ) ≤ α(µ)
because we have therefore:

||u -u|| ≤ ||r|| X ′ α(µ)

An easy case

We should first address the case where Θ q (µ) ≥ 0 ∀q = 1, . . . , Q, ∀µ ∈ .

If this condition is met, then for any µ ∈ , µ ∈ :

Q q=1 Θ q (µ)a q (u, v) ≥ Q q=1 Θ q (µ) Θ q (µ) Θ q (µ)a q (u, v) ≥ inf q ′ =1,...,Q Θ q (µ) Θ q (µ) Q q=1 Θ q (µ)a q (u, v) ≥ inf q ′ =1,...,Q Θ q (µ) Θ q (µ) α(µ)||u||||v||
This leads us to the following offline/online procedure for α(µ):

• Offline:

choose a "nominal" value µ ∈ for which Θ q (µ) = 0 ∀q = 1, . . . , Q ;

compute and store α(µ) by solving an eigenproblem on X .

• Online, given a value of µ, take

α(µ) = inf q=1,...,Q Θ q (µ) Θ q (µ) α(µ)
This procedure is simple and efficient, with a dim X -independent online stage.

However, if our positivity hypothesis on Θ q is violated, we have to resort to a more complex method, which we describe below.

The general case: successive constraints method (SCM)

The successive constraints method (SCM) is presented in [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF].

α(µ) can be rewritten

α(µ) = inf w∈X ,||w||=1
a µ (w, w)

We further note that:

α(µ) = inf w∈X ,||w||=1    Q q=1 Θ q (µ)a q (w, w)    = inf y∈    Q q=1 Θ q (µ) y q    (12)
where:

= ( y 1 , . . . , y q ) ∈ R Q st. ∃w ∈ X , ||w|| = 1 and y q = a q (w, w)∀q = 1, . . . , Q

To obtain a lower bound on α(µ), we should replace by a ⊇ ; in other words, we should relax the constraints on the optimization problem (12).

To define , we begin by noticing that is bounded: let:

σ - q = inf w∈X ,||w||=1 a q (w, w), σ + q = sup w∈X ,||w|||=1 a q (w, w), = Q q=1 [σ - q ; σ + q ].
Now let y = ( y 1 , . . . , y Q ) ∈ and let w satisfying the properties in the definition of . For every µ ′ ∈ , we have:

q Θ q (µ ′ ) y q = q Θ q (µ ′ )a q (w, w) ≥ inf w ′ ∈X ,||w ′ ||=1 q Θ q (µ ′ )a q (w ′ , w ′ ) = α(µ ′ )
So far we have:

⊂    y ∈ | Q q=1 Θ q (µ ′ ) y q ≥ α(µ ′ ) ∀µ ′ ∈   
The right-hand side of this inclusion would be a good candidate for ; however it is defined by an infinite number of µ-independent constraints (each for one µ ′ ∈ ). We should relax one more time these constraints in order to get a finite number of (now, µ-dependent) constraints. We take:

=    y ∈ | Q q=1 Θ q (µ ′ ) y q ≥ α(µ ′ ) ∀µ ′ ∈ M (µ, )   
Where M ∈ N, is a finite subset of and S M (µ, ) is the set of the M points in that are closest (with respect to some metric) to µ.

So our SCM offline/online procedure is:

• offline:

choose M ; choose (we postpone discussion about the choice of until the next subsection) ;

for each µ ′ ∈ , compute and store α(µ ′ ) by solving an eigenproblem on X ; compute and store σ - q and σ + q the same way, for each q = 1, . . . , Q.

• online: solve optimization problem:

α(µ) = inf y=( y 1 ,..., y Q )∈    Q q=1 Θ q (µ) y q   
The online optimization problem is what is called a linear programming problem: the goal is to minimize a linear function under a finite number of linear inequalities constraints. In our case there are Q variables and Q+M constraints.

There exist efficient algorithms, such as the simplex algorithm (see [START_REF] Nocedal | Numerical optimization[END_REF] for instance), which solve such optimization problems under (on average) polynomial complexity with respect to the number of variables and number of constraints, even if they can be exponential in the worst cases. The key here, once again, is that online complexity is still independent of dim X .

A last remark we can do on the algorithm is about the trade-off in the choice of M : whatever M is, we always get a certified bound on α(µ), but increasing M will improve sharpness of this bound, at the expense of an increase in online computation time.

Choice of

We will use a so-called greedy procedure for choosing the finite subset ⊂

. Roughly speaking, a decision algorithm is said to be greedy if he makes the best possible choice at each step. This strategy does not always yield a globally optimal decision but often gives a good solution for cases in direct global optimization would be unpractical.

For our greedy choice of , we will start with = {µ 1 }, with an arbitrarily chosen µ 1 ∈ . Then, at each step of the algorithm, we are going to add a point in that gives the best possible expected improvement in the precision of our bound α(µ).

Thus we get the following procedure for choice of :

1. choose M ;

2. initialize = {µ 1 } with some arbitrary µ 1 ∈ ;

3. choose a finite-sized sample set Ξ ⊂ ;

4. repeat:

• add

µ * = argmax µ∈Ξ α(µ) -α(µ)
to (here α(µ) stands for the SCM lower bound computed using the "current" )

The repeat loop can be stopped either when # has reached a maximal value, or when max µ∈Ξ α(µ)α(µ) is less than a desired precision.

The greedy algorithm written below would require computation of α(µ) for This cheap upper bound could be built using the following: let us define:

#Ξ × (1 + 2 + . . . + (# - 1
up = { y * (µ k ), k = 1, . . . , K} where {µ 1 , . . . , µ K } = and y * (µ k ) = arginf y∈    Q q=1 Θ q (µ k ) y q    (k = 1, . . . , K)
We can define:

α up (µ) = inf y∈ up    Q q=1 Θ q (µ) y q    which satisfies α up (µ) ≥ α(µ) because up ⊂ .
Using α up (µ) as a surrogate for α(µ) reduces fixed complexity of the greedy algorithm to K = # eigenproblems on X (for computation of the y * (µ k )) and #Ξ trivial minimizations per step over reasonably-sized, finite, up .

Choice of a reduced basis

In this section, we present two methods for choosing our reduced basis {ζ 1 , . . . , ζ N }.

These two methods have different offline complexities and lead to different bases. Both nicely fit into our reduced basis framework, and can benefit from the same procedure for error estimation, detailed in section 3.

Greedy procedure

The greedy procedure is treated in [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF].

The following algorithm is based on the same heuristic than the one described earlier for choosing in the SCM: at each step of the algorithm, we put one new function in the reduced basis; this function is the solution u µ * of the PDE for a specific value of the parameter choosed so as to give the best expected improvement in the online error bound, i.e. the value of the parameter for which we make the greatest (estimated) error when using "current" reduced basis.

The greedy algorithm for choice of the reduced basis is the following: given N ∈ N, the desired final size of the reduced basis,

• Choose a finite-sized, random, large Ξ ⊂ sample of parameters.

• Choose µ * 1 ∈ at random, and set

ζ 1 = u µ * 1 .
• Repeat, for n from 2 to N :

-For each µ ∈ Ξ, compute the online error bound for u µ when using {ζ 1 , . . . , ζ n-1 } as reduced basis.

-Find µ * n ∈ Ξ with the greatest online error bound.

-Set ζ n = u µ * n .
This procedure requires computation of N actual solutions, and #Ξ × N online error bounds. It relies crucially on the existence of a cheap online error bound.

Of course, instead of fixing a target size for the reduced basis, one can also keep iterating and adding functions to the reduced basis until the greatest online error bound (on the Ξ sample of parameters) gets below some prescribed value.

in the introduction, is large. We enumerate the snapshots parameters: Ξ = {µ 1 , . . . , µ N snap }.

Each snapshot u(i) = u µ i (i = 1, . . . , N snap ) has an expansion in the basis of X :

u(i) = j=1 u(i) j φ j
and we can write down M , the matrix of the u(i) j coefficients:

M =        u(1) 1 u(2) 1 . . . u(N snap ) 1 u(1) 2 u(2) 2 . . . u(N snap ) 2 . . . . . . . . . . . . u(1) u(2) . . . u(N snap )       
Denote by Ω, the (symmetric, positive definite) matrix of the scalar product

<, > in X .
We can write down our objective function as:

µ∈Ξ < u µ , p > 2 = p T Ω T M M T Ωp
The condition of optimality for (13) is then given by the Lagrange multipliers theorem: there has to exist λ ∈ R, so that:

Ω T M M T Ωζ 1 = λΩζ 1 ( 14
)
and λ is the optimal value of the problem (as can be seen by left-multiplying (14) by ζ T 1 , and using ||ζ 1 || 2 = 1) -which should be made as large as possible. Now left-multiply (14) by Ω -T = Ω -1 :

M M T Ωζ 1 = λζ 1 ( 15
)
so that ζ 1 is the unit eigenvector associated with the greatest eigenvalue of M M T Ω.

So we have reduced our problem as an eigenproblem of dimension , because

M M T Ω is of size × . Since can be large, this could turn out to be highly unpractical. Fortunately, we are to see that we can solve the same problem using eigenproblem of size N snap × N snap .

The key is to look for ζ 1 as a linear combination of the snapshots, instead of the "generic" basis items {φ i }.

We write:

ζ 1 = N snap i=1 z 1i u(i) = M z 1
where

z 1 = (z 1 1, . . . , z 1,N snap ) T . Setting ζ 1 = M z 1 into (15) yields: M M T ΩM z 1 = λM z 1
and a sufficient condition for this to hold is: Proof: To see 1., first note that if we had M z = 0, then M T ΩM z = 0 and so z could not be an eigenvector associated with a nonzero eigenvalue of M T ΩM ;

M T ΩM z 1 = λz 1 so that our new unknown z 1 ∈ R N snap is
thus M z = 0. Now M M T Ω(M z) = M (M T ΩM z) = M λz = λ(M z) so that M z
is an eigenvector of M M T Ω associated with λ. Point 2. is 1. mutatis mutandis. Point 3.: it is clear from 1. and 2. that the two matrices have the same nonzero eigenvalues. To get the statement about the multiplicities, suppose that Span{z 1 , . . . , z k } is the eigenspace of M T ΩM associated with eigenvalue λ. Since M T ΩM is self-adjoint with respect to the Euclidean inner product (i.e., is symmetric), we can suppose that {z 1 , . . . , z k } is orthogonal (for the Euclidean inner product). Now, from 2., Span{M z 1 , . . . , M z k } is contained in the eigenspace of M M T Ω associated with λ. Moreover, for every i = j:

< M z i , M z j >= z T j M T ΩM z i = λz T j z i = 0 so that {M z 1 , . . . , M z k } is <, >-orthogonal
, and therefore linearly independent, and so λ has multiplicity k in M M T Ω.

Second mode

Now let's say we are looking for a set of two functions as our reduced basis:

{ζ 1 , ζ 2 }. As earlier, we want (ζ 1 , ζ 2 ) to minimize: µ∈Ξ u µ -(< u µ , ζ 1 > ζ 1 + < u µ , ζ 2 > ζ 2 ) 2 among the (ζ 1 , ζ 2 ) ∈ X × X satisfying: ||ζ 1 || = ||ζ 2 || = 1 and < ζ 1 , ζ 2 >= 0.
In other words, we aim to minimize the overall error (over our sample of parameters Ξ) made when orthogonally projecting u µ onto Span{ζ 1 , ζ 2 }, with

{ζ 1 , ζ 2 } an orthonormal family.
Mimicking what we have done earlier, we rewrite our objective function as: Optimization in ζ 2 is the following:

µ∈Ξ u µ -(< u µ , ζ 1 > ζ 1 + < u µ , p > p) 2 = µ∈Ξ ||u µ || 2 -2 µ∈Ξ < u µ , ζ 1 > 2 -2 µ∈Ξ < u µ , p > 2 + µ∈Ξ < u µ , ζ 1 > ||ζ 1 || 2 + µ∈Ξ < u µ , p > 2 ||p|| 2 + µ∈Ξ < u µ , ζ 1 >< u µ , ζ 2 >< ζ 1 , ζ 2 > = Constant - µ∈Ξ < u µ , ζ 1 > 2 - µ∈Ξ < u µ , ζ 2 > 2
ζ 2 = argmax p∈X ,||p||=1,<p,ζ 1 >=0 µ∈Ξ < u µ , p > 2
The objective function has exactly the same matrix expression as the one of the problem defining the first POD mode:

µ∈Ξ < u µ , p > 2 = p T Ω T M M T Ωp
The Lagrange multiplier theorem now gives that there exist λ 2 ∈ R so that:

Ω T M M T Ωζ 2 = λ 2 Ωζ 2 (16)
As before, left multiplying by ζ T 2 gives that λ 2 is the optimal value of the problem, which should be made as large as possible.

Left-multiplying (16) by Ω -T we get that ( 16) is equivalent to:

M M T Ωζ 2 = λζ 2 so λ 2 is an eigenvalue of M M T Ω.
Now we have, for every x, y ∈ X :

< M M T Ωx, y > = y T ΩM M T Ωx = y T M M T Ω T Ωx = < x, M M T Ω y >
so that M M T Ω is self-adjoint with respect to <, >. So we can find a unit eigenvector ζ 2 , orthogonal to ζ 1 associated with the second greatest eigenvalue (note that this eigenvalue is equal to λ 1 if λ 1 has multiplicity > 1).

As earlier, we can use the fact that M M T Ω and M T ΩM have the same nonzero eigenvalues, to get ζ 2 = M z 2 where z 2 is an eigenvector of the smaller matrix M T ΩM . The fact that multiplicities of nonzero eigenvalues in M M T Ω and M T ΩM are the same ensures that this procedure is correct even if λ 1 = λ 2 .

Following modes, and summary for POD

We can readily repeat what we have done for ζ 2 to prove:

Let Ξ = {µ 1 , . . . , µ N snap } be a finite part of . Let M be the matrix whose columns are the coefficients of u µ 1 , . . . , u µ N snap with respect to some basis of X .

Let Ω be the matrix of the inner product in X .

Then, the family (ζ 1 , . . . , ζ N ) of vectors in X that minimizes:

µ∈Ξ u(µ) -< u µ , ζ 1 > ζ 1 + < u µ , ζ 2 > ζ 2 + . . . + < u µ , ζ N > ζ N 2 under the constraints < ζ i , ζ j >= 1 if i = j 0 else.
20 is defined by:

         ζ 1 = 1 ||M z 1 || M z 1 . . . ζ N = 1 ||M z n || M z n
where z k (k = 1, . . . , N ) is a nonzero eigenvector, associated with k-th largest eigenvalue (counting repeated nonsimple eigenvalues) of M T ΩM .

The statement above holds for every N ∈ N * so that the N -th largest eigenvalue of M T ΩM is nonzero.

The functions (ζ 1 , . . . , ζ N ) defined above are called the first N POD modes associated with the snapshot ensemble Ξ.

A note about the (offline) complexity of the POD: for a N -sized reduced basis, it requires computation of N snap "expensive" solutions u µ and extraction of the N leading eigenelements of an N snap × N snap symmetric matrix.

Case of a scalar output

One is often not interested in the full solution u µ , but rather in an "output" calculated from u µ :

s µ = η(u µ )
where η : X → R is a (µ-independent) functional.

We present here some approaches for computation of a surrogate output s µ from u µ .

Primal approach

A trivial approach to the question is to compute u µ , and then, online:

s µ = η( u µ ) (17)
Done naively, this procedure might be dim X -dependent, since evaluation of η may require access to every component of u µ in the original X space.

A way of overcoming this problem, if η is a linear functional, is to write that:

s µ = η N n=1 c n (µ)ζ n = N n=1 c n (µ)η(ζ n )
21 and the offline/online procedure becomes:

• offline: compute and store η(ζ n ), n = 1, . . . , N

• online: find c n (µ) (n = 1, . . . , N ) satisfying (2) using procedure described in 2.2 and compute

s µ = N n=1 c n (µ)η(ζ n ) (18)
And now the online phase is still dim X -independent.

For a quadratic output functional, one proceeds the same way using:

s µ = N n=1 N n ′ =1 c n (µ)c n ′ (µ)η b (ζ n , ζ n ′ )
where η b is the symmetric bilinear form associated with η.

This time, one has to compute and store (offline

) η b (ζ n , ζ n ′ ), for n, n ′ = 1, . . . , N .
The same procedure can easily be generalized for polynomial outputs.

Error bound for the output

A straightforward majoration for |s µs µ | can be found under a Lipschitz hypothesis for η: let L > 0 be so that, whatever v, w ∈ X :

|η(v) -η(w)| ≤ L||v -w||
(for example, if η is a continuous linear functional, L can be chosen as the dual norm of η).

Then we have:

|s µ -s µ | ≤ Lǫ µ (19)

Primal/dual approach

For linear outputs, error bound (19) can be improved, at the expense of, in some cases, an increase in computation time and program complexity. The idea, present in [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF], is to introduce the dual problem of (1) associated with η:

a µ (v, ψ µ ) = -η(v) ∀v ∈ X ( 20
)
whose unknown is ψ µ ∈ X .

we have:

r µ ( ψ µ ) = N du p=1 q n (µ)r µ (ζ du p ) (23)
and, for each p = 1, . . . , N du , we have, from definition of r µ , affine expansion of a µ and reduced basis expansion (2) of u µ :

r µ (ζ du p ) = f (ζ du p ) - Q q=1 Θ q (µ) N n=1 c n (µ)a q (ζ n , ζ du p ) (24)
So our full offline/online procedure for primal/dual evaluation and error bound on s µ reads:

• offline:

perform offline-stage RB approximation (part 2.2) and error estimation (part 3) for the two problems ( 1) and ( 20

Output-aware choice of basis: goal-oriented POD

When one is interested in a special output of u µ , it could be wise to choose an adapted reduced basis for this output. For instance, if the output depends strongly on the values of u µ (x) for x in some region R 1 contained in the domain of calculus, and much less strongly for x in some region R 2 , then a well-chosen reduced basis should lead to a precise u µ approximation in R 1 , which could be much less precise in R 2 .

The goal-oriented POD method, presented in [START_REF] Bui-Thanh | Goaloriented, model-constrained optimization for reduction of large-scale systems[END_REF] in the context of linear-timeinvariant (LTI) dynamical systems and balanced truncation, is a variation on POD taking the output functional in account. We believe it is possible to adapt the optimization framework described in [START_REF] Bui-Thanh | Goaloriented, model-constrained optimization for reduction of large-scale systems[END_REF] to solve (a weaker variant) of the above problem for (ζ 1 , . . . , ζ N ). However this optimization problem seems rather difficult (non-convex, with local non-global minima), thus requiring much more computational burden than classic POD.

Conclusion

We have presented a PDE-based, deterministic approach for model reduction called reduced basis approximation. The advantages of this approach are: a certified, computationally-efficient bound for reduction error, and a wide flexibility in the reduced basis chosen so as to enable output-aware choices of basis. Its main drawback is the need to make assumptions on the model that is being reduced. However work have been done to extend the reduced basis, and its associate procedures (error estimation, basis selection) to non-affinely

  )) values of µ. If one choose a large Ξ set, the cost of the choice of might be too prohibitive. Rather, one could consider a cheap upper bound α up (µ) for α(µ) and replace α(µ)α(µ) by the relative surrogate sharpness estimator:α up (µ)α(µ) α up (µ) .

  an eigenvector of M T ΩM (which of size N snap × N snap ) associated with its greatest eigenvalue. The absence of optimality loss when looking for ζ 1 in Range(M ) (in other words, ζ 1 = M z 1 with z 1 as defined below satisfies (13)) follows clearly from these facts: 1. if z is a nonzero eigenvector of M T ΩM (associated with λ = 0), then ζ = M z is a nonzero eigenvector of M M T Ω associated with λ ; 2. if ζ is a nonzero eigenvector of M M T Ω (associated with λ = 0), then z = M T Ωζ is a nonzero eigenvector of M T ΩM associated with λ; 3. M T ΩM and M M T Ω have the same nonzero eigenvalues, each with the same multiplicity in one and the other.

because of the constraints on ζ 1 and ζ 2 .

 2 The objective function, apart from the terms that are independent of ζ 1 andζ 2 ,is a sum of a function dependent on ζ 1 only and a function dependent on ζ 2 only. Thus the joint optimization on (ζ 1 , ζ 2 ) is equivalent to successive optimization on ζ 1 , then on ζ 2 . Optimization on ζ 1 has been done previously.

) ;-

  compute and store η(ζ n ) for n = 1, . . . , N ; -compute and store f (ζ du p ) (p = 1, . . . , N du ) and a q (ζ n , ζ du p ) (q = 1, . . . , Q; n = 1, . . . , N ; p = 1, . . . , N du ) ; • online: -perform online-stage RB approximation and error estimation for the two problems (1) and (20); this gives respective expansions {c n (µ)} n=1,...,N and {q p (µ)} p=1,...,N du of u µ and ψ µ in, respectively primal reduced basis {ζ 1 , . . . , ζ N } and dual reduced basis {ζ du 1 , . . . , ζ N du }, together with α(µ), ǫ µ (error bound for primal problem) and ǫ µ (error bound for dual problem) ; -compute η( µ) as left-hand side of (18) ; -compute r µ ( ψ µ ) using (23) and (24) ; -output s µ using (21) ; -give error bound (22).

Θ

  Goal-oriented POD proceeds by finding (ζ 1 , . . . , ζ N ) ∈ X minimizing the sum of squared differences (for different values of a parameter in a finite sample #Ξ) between the actual output and the output evaluated on the RB approximation taking {ζ 1 , . . . , ζ N } so as to minimizeinf µ∈Ξ η(u µ )η( u µ ) q (µ)a q (ζ n , ζ p )   c n (µ) = f (ζ p ) ∀µ ∈ Ξ, ∀p = 1, . . . , N < ζ n , ζ m >=    1 if n = m 0 else∀n, m = 1, . . . , N

parametrized PDEs ( [2]) or non-linear, time-dependent PDEs (see [4] for instance).

Future work should include further extensions to other PDEs with different parameters, and more output-aware basis choice procedures.

Proper Orthogonal Decomposition (POD)-driven procedure

To our knowledge, the following presentation of the Proper Orthogonal Decomposition is new. We refer to [START_REF] Volkwein | Proper orthogonal decomposition and singular value decomposition[END_REF] for a different presentation, based on singular value decomposition (SVD).

First mode

Suppose, for a soft start, that we want to find a reduced basis containing only one item: {ζ 1 }. Without loss of generality we can always suppose that ||ζ 1 || = 1. It seems reasonable to make the choice which makes the least summed squares error when projecting orthogonally the solutions u µ (for all µ ∈ ) onto the reduced space Span {ζ 1 }. Said with formulas, we seek ζ 1 satisfying:

The function ζ 1 will be called our first POD mode.

Unfortunately, the integral appearing in the objective function below is, in general, analytically untractable. We can approximate it using a finite sample of parameters Ξ ⊂ ; the integral then gets replaced by a finite sum over Ξ:

The set of solutions {u µ , µ ∈ Ξ} is called the snapshots ensemble.

The objective function can then be rewritten:

And so:

The objective function in (13) is a quadratic form in p, that has to be maximized under a "unity norm" constraint. We know this can be reduced to an eigenproblem.

More specifically, let us introduce a basis {φ 1 , . }, find an approximation ψ µ ∈ X du to ψ µ satisfying:

with error bound:

After having found reduced basis solutions for the primal and dual problems, the following surrogate output can be taken:

The error bound on the output now reads: Proof of (22): We omit the µ subscripts.

= |r(ψψ)| by the very definition of r ≤ ||r|| X ′ ||ψ -ψ|| ≤ ||r|| X ′ ǫ du = αǫǫ du from the definition of ǫ.

We shall now pass to the offline/online procedure: we first require computation of r µ ( ψ µ ) in (21).

Expanding ψ µ on the dual reduced basis: