HOW TO SAMPLE TO CALIBRATE A SPATIAL MODEL: APPLICATION FOR MODELLING VINE WATER STATUS

Keywords: sampling, vine water status, ancillary information sources, modelling Mots -Clés : échantillonnage, état hydrique de la vigne, variables auxiliaires, modèle spatiale

Several studies have shown that vineyards are spatially variable at within field scale. Characterizing this variability could be used to improve grape and wine production, by adapting management practices and harvesting to the vineyard variability. The goal of this work is to calibrate a spatial model by coupling high quality, high cost punctual vine measurements with low cost medium-high density ancillary data sources. For operational constraints, the number of measurements required to calibrate the model have to be as small as possible. The approach we propose is tested on a particular but significant example: the spatial variability of vine water status. Exposed leaf area (ELA), trunk circumference (TC) and Pre-dawn leaf water potential (PLWP) were measured in 49 sites in a non-irrigated Syrah block. A sampling method was applied to rank these sites according to their suitability to represent the statistical variability on ELA and TC. The first selected sites were used for the calibration of a model relating ELA and TC with PLWP. The model was validated with the remaining sites. The proposed method allowed the model calibration using only 3 sites. From 3 to 9 sites, the model performance improved for each additional site used for calibration from a Standard Error of Prediction (SEP) of 0.12 MPa to 0.10 MPa. The block statistical variability was well represented with the first 9 chosen sites. This study proves the suitability of the approach to calibrate a spatial model describing plant water status spatial variability. The obtained site selection showed to be well suited to represent the within field variability and optimize site number and location for the calibration of the model. The work deals with an important operational constraint for the commercial application of spatial models: reducing the number of measurements needed for their calibration. This is a first, essential step for vine growers to be able to spatially optimize vine management and harvest.

INTRODUCTION

It has been extensively shown that vineyards are spatially variable at within block scale [1]. Differential management of this variability could allow optimization of vineyard management practices, and consequently to differentiate and improve grape and wine quality [2]. There is, nevertheless, an important limitation for identifying spatial variability in a vineyard: the cost or impracticality of taking direct measurements with enough spatial resolution. Measurements of soil properties, plant water status, nutrient or health status, quality attributes or even measurements of biomass are usually too costly and time consuming to be directly measured with enough spatial resolution. One extended solution is to use ancillary information, easier and more affordable to measure, in order to estimate high quality, high cost related attributes [3]. Between those high quality vine measurements, plant water status has been extensively indicated as a crucial parameter greatly affecting grape and wine quality and production [4][5][6]. In this regard, Acevedo-Opazo et al. [7], proposed a conceptual spatial model based on ancillary information sources to estimate vine water status spatial variability at within field scale. Nevertheless, to be operational, even if a spatial model has been developed an important constraint for its practical application is its calibration for a particular field or conditions. If the variable to estimate is expensive and/or too time-consuming, the model could only be used operationally if it can be calibrated with a small number of measurements. In some cases, the affordable sample size could be as low as 3 to 5 measured locations per field. Aiming to overcome this constraint, the goal of this work is to propose a sampling method to optimize the location of high cost vine measurements in order to calibrate a regression spatial model with a very small number of measured sites (3)(4)(5)(6)(7)[START_REF] Coulouma | Proximal soil sensing[END_REF][START_REF] Scholander | [END_REF][START_REF] Murisier | [END_REF]. The spatial model proposed by Acevedo-Opazo et al. [7] to estimate vine water status spatial variability at within field scale was used as a case study.

MATERIALS AND METHODS

Study Case: Spatial model for the estimation of plant water status

According to the work by Acevedo-Opazo et al. [7], a linear regression model relating ancillary data sources and plant water status is considered in the current work. Based on the results of this study, growth and vigour vine measurements were used as ancillary information sources for estimating plant water status.

Site description

The experiment was carried out on a 1.2 ha unirrigated Syrah field , established in 1990 with 1 m spacing between vines and 2.5 m between rows, with the vines trained in a vertical shoot positioning system. The soil is considerably variable in soil water availability and vine growth [START_REF] Coulouma | Proximal soil sensing[END_REF]

Measurements

The database presented by Acevedo-Opazo et al., [7] comprised plant water status measurements as well as Exposed Leaf Area (ELA), and Trunk Circumference (TC) was used in this work. All variables were measured on 49 sites, arranged in a regular grid. Each site was identified with a georeferenced position in the field with a stand-alone GPS receiver (Garmin Etrex).Vine water status was measured as predawn leaf water potential (PLWP) between 3:00 and 5:00 a.m. using a pressure chamber [START_REF] Scholander | [END_REF] on a date corresponding to a significant water restriction. Three adjacent vines were measured at each sampling site and averaged to obtain a site value (si). Exposed Leaf Area (ELA) was measured manually using the method of Murisier and Zufferey [START_REF] Murisier | [END_REF] in order to assess the seasonal vegetative growth. Trunk circumference (TC) was measured 10 cm above the graft of each target vine, as an integrative attribute taking into account the average vigor of the vines from the date of plantation.

Sampling method

Kennard and Stone method (K&S) [11], based on experimental design theory and aimed for linear regression models is proposed for the calibration of the spatial model. The approach of the method is to select the sites to be used for calibrating the model according to their values on the ancillary variables (i.e ELA and TC). The criterion for the site selection is to maximize the variance gathered by the sampled sites according to their values on the ancillary variables. The method provides a rank of sites. This rank describes the order in which the different sites need to be sampled in order to maximize the variance for each sample size, this is, in order to optimize the calibration of the model for any given affordable number of PLWP measurements.

Mapping

Mapping was done with the 3D field software (Version 2.9.0.0, Copyright 1998-2007, Vladimir Galouchko, Russia). Inverse distance weighting was used for interpolation because of the limited number of data (n=49) in the field.

RESULTS AND DISCUSSION

Spatial model

The spatial model was able to explain most of the variation on PLWP in the parcel (R2=0.85) when using all available 49 measurements for the calibration of the model. The Standard Error of Calibration (SEC) was 0.1 MPa. This value represents best expected performance of the chosen spatial model to explain the within block spatial variability of PWLP.

Results of the sampling method

The proposed method allowed the model calibration using only PLWP measurements taken on 3 sites. From 3 to 9 sites, the model performance improved for each additional site used for calibration from a Standard Error of Prediction (SEP) of 0.12 MPa to 0.10 MPa. The SEP represents the error when applying the calibrated model to estimate plant water status in sites of the parcel where only auxiliary measurements are available. The standard deviation (STD) of PLWP within the block was 0.258 MPa. The STD value accounts for the error produced when the average PLWP of the block is considered as a general estimation of PLWP for every site. Thus, the spatial model calibrated with only 3 sites provides in this case a significant improvement over using the average of the block as a general estimation of PLWP for every site.Figure 1 a) shows plant water status (PWLP) zones estimated in the parcel through the spatial model calibrated with the 3 sites selected by K&S. When compared to the zones obtained with 49 PLWP measured values, shown in Figure 1 b), similar spatial structure can be observed in both cases.

CONCLUSIONS

This work dealt with an important operational constraint for commercial application of spatial models in viticulture which is:" how to reduce the number of measurements needed for the calibration of a spatial model?"

The case study here focused on a linear regression spatial model based on ancillary information sources. The spatial model was used to estimate the spatial variability of vine water status at a within field scale. This case study is relevant for the wine industry because of its commercial implication in terms of vineyard water status management and its direct impact on yield and wine quality.

A sampling algorithm based on experimental design was tested for the selection of sites to be used in the calibration of the model. The method provides a rank of sites based on one or more ancillary information sources (AIS), i.e, plant based measurements such as trunk circumference and exposed leaf area. Measurements of plant water status taken at these selected sites can then be used to calibrate a model that relates plant water status with relevant AIS. The resulting model can be applied to the rest of the field at each location where AIS are available. Once the spatial structure of the field has been highlighted by AIS, then, the information provided by selected sites can be used to estimate the spatial variability of plant water status within the field. The current work proves that such a spatial model can be calibrated using a very reduced number of measurements (3 measurements of plant water potential taken at the time when plant water status of the field needs to be characterized). The calibrated model allows the delineation of different homogeneous within-field zones. This information could be used as a decision tool for improving the management of field operations like irrigation. This method can be extended to other variables provided the availability of relevant AIS with a convenient resolution adapted. In general terms, this approach could be used to better characterize the spatial variability of key production parameters (such as yield, water status of quality parameters like sugar content) at spatial resolutions that are cost prohibitive to generate with targeted measurements. Thus, this work represents a first step towards the optimization of vineyard management practice for routine practices in a commercial context.

Figure 1 :

 1 Figure 1 : Observed (a) and estimated (b) PLWP values using the spatial model calibrated with 3 sites, selected through K&S. the Geographical location of the selected sites is marked in the figure.

Author-produced version of the article published in IXth International Terroir Congress, Dijon, 25-29 Juin 2012


Mickaël Fabrègue
email: mickael.fabregue@teledetection.fr

Agnès Braud
email: agnes.braud@unistra.fr

Sandra Bringay
email: bringay@lirmm.fr

Florence Le Ber
email: florence.leber@engees.unistra.fr

Maguelonne Teisseire
email: maguelonne.teisseire@teledetection.fr

Including spatial relations and scales within sequential pattern extraction

Keywords: Data Mining, Sequential patterns, Spatio-temporal, Aquatic ecosystem

Georeferenced databases contain a huge volume of temporal and spatial data. They are notably used in environmental analysis. Several works address the problem of mining those data, but none are able to take into account the richness of the data and especially their spatial and temporal dimensions. In this paper, we focus on the extraction of a new kind of spatio-temporal pattern, considering the relationship between spatial objects and geographical scales. We propose an algorithm, STR PrefixGrowth, which can be applied to a huge amount of data. The proposed method is evaluated on hydrological data collected on the Saône watershed during the last 19 years. Our experiments emphasize the contribution of our approach toward the existing methods.

Introduction

Due to the recent explosion of mobile technologies and georeferenced data, a new kind of data has emerged: spatio-temporal data. Each data is associated with a given spatial reference (i.e. a localisation) and a temporal information (i.e. a timestamp). New needs for the monitoring of these data in time and space have appeared, for example to study the spread of information in social networks [START_REF] Lin | Inferring the diffusion and evolution of topics in social communities[END_REF], in epidemic surveys [START_REF] Gubler | Epidemic dengue/dengue hemorrhagic fever as a public health, social and economic problem in the 21st century[END_REF] or for hydrological monitoring as presented in this article. In these domains the volume of data is huge, and commonly contains heterogeneous information. Often several levels of spatial division describe the geographical aspect, based on an inclusion property and relationships between geographical objects. An area can be included in another area (e.g. Europe is divided into countries Spain, France, Germany, etc.). Moreover, the geographical objects are linked by spatial relations. For example, an area is close to another area, is located to the north or to the east of another area (e.g. USA and Canada are two adjacent areas and USA is south of Canada). In this article we will focus on data mining methods which consider the temporal dimension and also spatial relationship between spatial objects. The objective is to provide a method for extracting spatio-temporal patterns to highlight common behavior in large volumes of data. The method is applied to the environmental domain, and more specifically for studying aquatic ecosystems. The dataset is a collection of samples on the hydrological catchment of the Saône, from part of the Fresqueau project. This project aims to provide operational tools to study the state of aquatic systems. It falls within the scope of the Water European Framework Directive which aims to correct the state of the aquatic systems and catchments in 2015.

Related work

Pattern extraction has been the subject of lot of research in the field of data mining. Pattern discovery highlights a recurring information in data, characterizing a frequent behavior. This knowledge can be represented by various types of patterns. Several authors have proposed new methods that consider both time and space.

In [START_REF] Wang | LNCS 3453 -Mining Generalized Spatio-Temporal Patterns[END_REF], data is represented in a set of spatial grids where items (events) appear at different coordinates. Each grid describes the state of the problem at a specific timestamps t. For each date and absolute position, an itemset (a set of events), is generated. For each absolute position, a sequence of itemset is built by considering all the timestamps. Then, sequential patterns are extracted from such sequences by considering an absolute position as the reference point. An example of a pattern obtained by such a method is xpRain(0,0)qpHumidity(0,1)qy meaning that it frequently rains at coordinates (0,0) and, later in time, humidity exists at coordinates (0,1). This kind of pattern has the disadvantage of being sensitive to the choice of the reference point. Furthermore the space is reduced to a grid representation.

In [START_REF] Huang | A Framework for Mining Sequential Patterns from Spatio-Temporal Event Data Sets[END_REF], the authors proposed the concept of close events in time and space. A spatio-temporal window is defined by both a temporal and a spatial interval.

Patterns are association rules such as xRainñHumidityy, meaning that in close areas at close timestamps, the rain is frequently followed by humidity. These rules do not take in consideration potential relationships between spatial objects nor different geographical scales.

The extraction of spatio-temporals patterns with neighborhood relationships between geographical objects is proposed in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF]. Patterns have the shape xpHumidity .[Rain Wind]qpHumidity Rainqy. Neighborhood relationships are denoted by a neighborhood operator . and a grouping operator []. Consider, for example, a city in which the previous pattern is found. This pattern means that humidity appeared at a timestamp and at the same time rain and wind appeared in a nearby town (according to an euclidean distance or defined by user). Later, humidity and rain appeared in the city. This spatial relationship is simple and it is not possible to specialize it, nor to have several levels of granularities (i.e. several geographical scales). Actually, it is limited to one kind of relationship: spatial proximity.

A technique for granularity management of space is provided in [START_REF] Tsoukatos | Efficient mining of spatiotemporal patterns[END_REF]. As in [START_REF] Wang | LNCS 3453 -Mining Generalized Spatio-Temporal Patterns[END_REF], a grid of events represents spatiality and a set of spatial grids represent temporality. The user has to choose a level of granularity that will merge a set of adjacent cases in the grid. The higher the level of granularity, the bigger the set of merged cases. This technique aims to generalize data in a spatial way. To extract patterns, it is necessary to choose a granularity value and furthermore a grid representation to describe spatiality. Extracted patterns have the shape of classical sequential patterns such as xpSunqpWindqpSun,Humidityqy, meaning that frequently Sun is followed by the event W ind, followed itself by the events Sun and Humidity according to a specific level of granularity.

All these methods do not effectively consider complex data with geographic objects linked together and at different scales. The approach in this paper aims to take into account all these notions: 1) by considering the temporal and spatial dimensions 2) by generalizing the problem to a more complex spatial relationships between geographical objects 3) by including all possible granularities during the extraction process.

In section 3.1, we introduce some preliminary definitions on which our method is based. Then sections 3.2 and 3.3 present a formal framework to take into account relationships between objects and different spatial granularities. The developed algorithm is presented in section 4. In section 5, the method is applied to a real dataset and the obtained results are presented. We discuss on the prospects of the scope of this proposal in section 6.

Spatio-temporal patterns

Our approach extends the notion of sequential patterns introduced in [START_REF] Agrawal | Mining sequential patterns[END_REF] and takes into account the temporal and spatial sequential patterns defined in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF].

Preliminaries

Sequential patterns are extracted from a set of data sequences. For each coordinate or geographical object, a sequence of events is built. First, we consider the database DB presented in the table 1 which shows the set of events which appeared in three different cities in the south of France. For each city, a sequence is generated (see. the table 2) Definition 1 (Sequence) Let I = tI 1 , I 2 , . . . , I m u be the set of items (events). An itemset is a non empty set of items denoted pI 1 , I 2 , . . . , I k q where I j I (it is a non ordered representation). A sequence s is a non empty ordered list of itemsets noted xIS 1 IS 2 . . . IS p y where IS j IS, with IS the set of itemsets.

Each sequence being composed of itemset, an item can appears several times in a same sequence. 2 where each represents an event sequence for a city, we note that the sequence S xpSunq pHeat=Highqy is supported by sequences S N îmes and S M ontpellier . Then S ¨SN îmes and S SMontpellier .

A sequential pattern is a frequent sub-sequence characterized by a support which is the number of occurrences of a pattern in S, a set of sequences. The extraction of those patterns is determined by a minimum support parameter denoted θ. This means that only patterns with a support value greater than θ will be extracted. Let M be the set of extracted sequential patterns, then dM M and SupportpM q ¥ θ.

Definition 3 (Sequential pattern support) A sequence S S supports a sequential pattern M when M ¨S. The support of a pattern M is the number of sequences in S in which M is included (supported). Let S I be the set of sequences that support m, then S I tS i S such as M ¨Si u and SupportpM q |S I |.

Example 2 Consider table 2, we note that sequence S xpSunq pHeat=Highqy is supported by sequences S N îmes and S M ontpellier . Therefore SupportpSq 2.

Although sequential patterns fit the temporal aspect well, they are not able to consider the spatiality nor potential relationships between geographical objects. To take into consideration these two aspects, we now present a hierarchical approach based on dimensions.

Relationships between spatial objects

Relationships between spatial objects are potential links that exist between geographic points or objects. For instance, in an epidemiological context, spatial objects are cities or areas. Therefore, several links should be considered, as migration flows or obstructions presence like a forest or a mountain. These links can also be specialized, e.g a forest can be a fir forest or an oak forest. Therefore important to take into account these links but also their potential specializations. To achieve this goal, we use a set of spatio-relational dimensions D R described by sets of values with an associated hierarchy for each dimension.

Definition 4 (Spatio-relational dimension) A dimension D D R is defined by a domain of values X j such as dompDq tX 1 , X 2 , ..., X n u.

Example 3 Let D Orientation be the spatio-relational dimension representing orientation according to the points of the compass with the addition of North-West, North-East, South-West and South-East: dompD Orientation q tNorth, W est, South, East, N orth-W est, N orth-East, South-W est, South-Eastu. Then the following hierarchy is constructed:

After defining the concept of a spatio-relational dimension, we introduce the notion of a hierarchy on a dimension. The ojective is to easily consider the more specific relations of an existing relation.

Definition 5 (Hierarchical representation of a spatio-relational dimension) Let D D R be a spatio-relational dimension with dompDq tX 1 , X 2 , ..., X n u and let H H R be the hierarchy representation associated to this dimension, then H is a semi-lattice or an oriented tree and for all node N H, labelpN q dompDq.

Orientation

Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q
To navigate in this hierarchy, we have to define some navigation operations as defined in [START_REF] Plantevit | Mining multidimensional and multilevel sequential patterns[END_REF]. These operations represent the notions of direct and global generalization and specialization. Definition 6 (Direct and global specialization) Let down R pX i q be the operation which allows access to direct specializations of the relation X i and downAll R pX i q the operation which allows access to all specializations of relation X i . The direct specializations of X i are X j such that there is a descending edge from X i to X j in the hierarchy, and the global specializations of X i are X k such that there is a descending path X i to X k .

Example 4 Lets consider the dimension D Orientation :

-down R pWestq tNorth-West, South-W estu, -downAll R pOrientationq dompd Orientation q.

Definition 7 (Direct and global generalization) Let up R pX i q be the operation which allows access to direct generalizations of the relation X i and upAll R pX i q the operation which allows access to all generalizations of the relation X i . The direct generalizations of X i are X j such as there is an ascending edge from X i to X j in the hierarchy, and X i global generalizations are all X k such as there is an ascending path from X i to X k .

Example 5 Using the same example:

-up R pNorth-Eastq tNorth, Eastu, -upAll R pNorth-Westq tNorth, W est, Orientationu.
This hierarchy offers the possibility to extract information at different levels. For instance, it can take into consideration the presence of an event at the north, but also drill down the hierarchy to find more specific relations. From the definition of the hierarchy on spatio-relational dimensions and their operations, patterns are extracted by considering relations between spatial objects. To add this new information to patterns, we use the link operator . when an item is found in a zone linked by a relation, as in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF]. When multiple items are considered by the operator ., the n-ary group operator [] is used.

Definition 8 (Related itemset)

Let IS and IS I be two itemsets which describe two differents zones Z and Z I at the same timestamps, if there exists a link δ in the spatio-relational hierarchy between Z and Z I , then they constitute a related itemset noted IS R pIS .δrIS I sq which means the itemset IS is found in Z and at the same time the itemset IS I appears in a closed zone in δ relation with the first zone.

Example 6 Taking two cities C 1 and C 2 , humidity appears in C 1 , rain and wind in C 2 at the same timestamps t. Furthermore the hierarchy highlights the fact that C 2 is at south of C 1 then the related itemset IS R pHumidity . South [Rain Wind]q is found in C 1 .

We have now to define the inclusion of a related pattern in another related sequential pattern. This inclusion is very close to the classic sequential pattern inclusion, the difference concerns the inclusion between itemsets.

Definition 9 (Inclusion of a related itemset)

A related itemset IS R IS i .δrIS j s is included in another related itemset IS I R IS I

i .δ I rIS I j s, if and only if, IS i IS I i , IS j IS I j and δ δ I or δ upAllpδ I q (i.e. δ I is equal to δ or δ I is a specialization of δ).

Example 7 Let DtOrientationu be the spatio-relational dimension of two itemsets IS 1 and IS 2 such that IS 1 . South rHumidity,Winds and IS 2 . South¡East r Humidity, Rain, Winds. We can note that all items in IS 1 are included in IS 2 and the relationship IS 1 is more general than IS 2 in the hierarchy. Therefore

IS 1 ¨IS 2 .
The obtained sequential patterns are composed of related itemsets and form a new kind of pattern, i.e. related sequential patterns.

Definition 10 (Related sequential pattern) Let IS be the set of itemsets and IS R the set of related itemsets, a related sequential pattern M R is a non-empty ordered list of itemsets and related itemsets denoted xIS 1 , IS 2 , . . . , IS p y where IS j IS IS R with a support value SupportpM R q.

In this section, we introduced a new kind a sequential pattern which considers existing links between geographical objects. These relations are organized in a hierarchy to efficiently consider specializations and generalizations. But in the context of spatial segmentation, it is also important to take into account the spatial granularity which exist between the zones in patterns in order to provide the experts with more precise patterns. The next section presents this new feature.

Geographical granularities in patterns

Different geographical granularities describe a division of space, itself divided into sub-divisions. This segmentation can have different shapes according to the context of the problem to solve. For example, let us consider a division of the Earth with respect to a geopolitical point of view. Space is divided according to continental frontiers or country boundaries. With respect to a climatic point of view, this division is different: hot climate areas, temperate areas, etc. In addition, areas are further divided into smaller regions. It is therefore necessary to not only take into consideration areas, but also their sub-divisions.

Definition 11 (Area dimension) An area dimension D D S is defined by a domain of values X j such as dompDq tX 1 , X 2 , ..., X n u. dompd Country q tAustria, Belgium, Bulgaria, ..., Sweden, U nited Kingdomu.

Considering granularity relies on the construction of a hierarchy based on an inclusion relation on such geographical divisions. To illustrate this, let us take as an example a division of the Earth by considering multiple granularities, continents and countries. The following hierarchy describes this division.

Example 9 Earth geopolitical division hierarchy

W orld

Asia

Af rica Europe ...

Q Q Q F rance Belgium ... Q Q Q Ain Aisne ... Q Q Q
This hierarchical representation of spatial granularities is close to hierarchies representing links between spatial objects (previous section). However, this hierarchy is not based on a generalization/specialization notion but on an inclusion order. For instance, a country is not a continent specialization but a sub-division.

We have to redefine navigation operations in this type of hierarchy H H S , with H S the set of hierarchies on granularity dimension. Definition 12 (Direct and global content) Let down S pX i q be an operation that allows access to the direct content of granularity X i and downAll R pX i q the operation which allows access to all the content of granularity X i . The direct content of X i is X j such that there is a descending edge from X i to X j in the hierarchy, the global content of X i is X k such that there is a descending path from X i to X k .

Example 10 Let us take the example of Earth division:

-down S pEuropeq tBelgium, F rance...u, -downAll S pEuropeq tBelgium, F rance..., Ain, Aisne...u.

Definition 13 (Direct and global containers) Let up S pX i q be an operation that allows access to the granularity that directly contains X i and upAll S pX i q the operation which allows access to all granularities that contain X i . The direct containers of X i are all X j such as there is a ascending edge from X i to X j in the hierarchy, and the global containers of X i are all X k such that there is an ascending path from X i to X k .

Example 11 Let us take the previous example:

-up S pBelgiumq tEuropeu, -upAll S pAinq tFrance, Europe, W orldu.

We use this hierarchy to add the notion of spatial inclusion into patterns

To extract these patterns, the algorithm navigates throughout granularity hierarchies and checks if a pattern is frequent at a more specific level of the hierarchy. If it is indeed the case, the pattern becomes spatio-temporal because its frequency depends on a specific spatial area.

Definition 14 (Spatio-temporal pattern)

Let v w be an operator of spatiality and M a classic or related sequential pattern, X k D the value of a granularity dimension D, a minimal support θ and S I the set of sequences S i such that |M ¨Si | at the granularity value X k . If |S I | ¡ θ then a spatio-temporal pattern M I is created, such that M I vX k wM.

Example 12 Let M xpHumidity. N orth rRainWindsqpHumidityRainqy be the relational pattern, with θ 10% and SupportpM q 50%. The pattern M has a frequency equal to 50% over the Earth but has a frequency equal to 15% if we just consider European cities. A spatio-temporal pattern M I is created such that M I vEuropewxpHumidity. N orth rRainWindsqpHumidityRainqy and SupportpM I q 15%

The previously presented definitions allow for taking into account spatial relationships and also geographical granularities. An adapted algorithm has been implemented to extract related spatio-temporal patterns at differents scales. This algorithm is presented in the next section.

STR PrefixGrowth algorithm

To extract patterns, we used the PrefixSpan [START_REF] Pei | Mining sequential patterns by pattern-growth: The prefixspan approach[END_REF] extraction algorithm, as was also used in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF]. This is currently one of the most efficient algorithms for extracting sequential patterns, both in terms of computation time and in terms of memory consumption. Sequential patterns are extracted from common prefixes.

For instance, xpaqy, xpaqpaqy, xpaqpabqy and xpaqpabcqy are prefixes of sequence xpaqpabcqpacqpdqpcfqy.

If a prefix is present in a number of sequences greater than a minimum support value θ, then this prefix is considered as frequent. When a frequent prefix is found, the database is divided recursively. When we look for frequent patterns, it is not necessary to keep the entire database and therefore data (i.e sequences) that do not support the current pattern are not preserved in the projected database. The reason is that these sequences will not support patterns of greater length because of the antimonotonic property of support. The efficiency of this algorithm is due to (1) the non-generation of candidate patterns thanks to research of frequent prefixes, and (2) the projection of the database into smaller databases to accelerate the exploration by removing sequences no longer needed. Algorithm 1: ST R P ref ixGrowthpα, θ, DB| α , D R , D S q input : α a pattern, θ a support minimum, DB| α projected database according to pattern α , D R a set of spatio-relational hierarchies, D S a set of granularity hierarchies output: SP set of patterns extracted in this function call (i.e current recursion)

I θ Ð getListOccurences(θ, DB| α , D R ); SP Ð ∅; foreach i in I θ do β=append(α,i); SP Ð SP β; SP Ð SP prefixGrowth ST M pβ, θ, DB| β , D R , D S q; SP Ð SP exploreSpatialHierarchypβ, θ, DB| β , D S q;
end Algorithm 2: getListOccurencespθ, DB| α , D R q input : θ a minimum support, DB|α projected database according to pattern α, DR a set of spatio-relationnal hierarchies output: I θ the list of frequent occurences in DB|α input : α a pattern , θ a minimum support, DB|α projected database according to pattern α , DS a set of granularity hierarchies output: SP the set of extracted patterns

I θ Ð I θ Y searchIExtend(θ,DB|α); I θ Ð I θ Y searchSExtend(θ,
SP Ð ∅; foreach dimi in DS do foreach s in dimi do if isFrequent(α, θ, DB|α,s) then /* check if a pattern is frequent in the current granularity */ SP Ð SPY spatialPattern(α, s); end end
Our general approach is described by the recursive algorithm 1, called ST RPref ixGrowth for Spatio Temporal and Relational PrefixGrowth . This method first determines the list of frequent occurrences in the database projected according to α and depending on the minimum support θ. A frequent occurrence (e.g. a frequent item) means that a pattern of greater length is found. In the function getListOccurencespq, we explore the relationship hierarchies. Two operations are used, the searchIExtendpq and searchSExtendpq, representing the two ways to extend a pattern, the I-Extension and the S-Extension. The I-extension adds an item to the last itemset of a sequence and the S-Extension adds a new item to a new itemset at the end of a sequence, at a further timestamp. For example let us take the pattern m xpaqpbqy and a frequent occurrence representing the item c. If c is an I-extension and m I an extended pattern, then m xpaqpbcqy. If c is an S-extension and m P an extended pattern, then m P xpaqpbqpcqy. For each relationship hierarchy, searchIExtendpq and searchSExtendpq operations are used to find occurrences of relations on every level of hierarchies. Frequent relations are then considered as occurrences. Relations between sequences are managed as individual items, they are returned along with occurrences of classics items. This function is provided by algorithm 2.

Occurrences, or frequent items, will be used to extend the pattern α with the function appendpq, which considers that an item is an intension or an extension. Then, for each extended pattern β, we project the database according to this pattern and we call P ref ixGrowth ST M to continue the recursive search of patterns. Finally, each pattern is given as a parameter of the function exploreSpatialHierarchypq that explores the spatial dimensions at all levels of granularity (algorithm 3) to find new patterns (section 3.3). For each spatial dimension, it checks if a pattern is frequent at each granularity of the hierarchies. If it does, we add the spatial pattern to the set of patterns.

The PrefixSpan complexity in the worst case is Θpp2 ¤Iq L q with I the number of items and L the length of the longest sequence in the database DB. Let H R be the number of hierarchies of spatial-relations, let R be the maximal number of relations per hierarchy of spatial relations, let H S be the number of spatial hierarchies and let S the maximal number of spatial areas per spatial hierarchy, the complexity of the STR PrefixGrowt algorithm is ΘpH S ¤ S ¤ p2 ¤ N ¤ H R ¤ Rq L q. This algorithm is pseudo-polynomial, i.e. is linear according to the number of extracted patterns. The worst case corresponds to the maximal number of patterns which could be extracted in a specific dataset.

To test and validate our method, we have applied this algorithm to a real dataset and we have compared it to existing methods. These results are presented in the following section.

Mining hydrological data

The dataset has been supplied by the RMC agency in the context of the Fresqueau project. It describes the biological and physicochemical information of streams in the Saône watershed, in the east of France. The data have been collected at different timestamps on 771 sites. The information contains different kinds of characteristics as biological indicators, pH, levels of nitrates or phosphates... For each site, a set of collected data for a specific timestamp is an itemset and those itemsets are ordered according to the time to generate a sequence.

Moreover, to apply our approach, we have selected some characteristics 1) to explicit the links between river sites and 2) to consider different geographical scales.

Hierarchies

Those data are described by several dimensions with their associated hierarchies to consider granularities and links between stations. They are presented as follows:

Stream orientation: this allows us to know whether a site is located downstream or upstream from another site. This is a simple hierarchy, one level deep.

Example 13 Stream orientation hierarchy

Stream orientation

Downstream

U pstream

Q Q Q
Hydrographic zones: France is divided into general watersheds and into three more specific partitions. Each level is a sub-division of the previous level. Watersheds are the most general level, itself divided into hydrographic areas. Then, there are sectors divided into sub-sectors. This hierarchy therefore has 4 levels.

Example 14 Hydrographic zone hierarchy

N ational level

Saone watershed level

Sectors level

Sub ¡ sectors level

F rance U 0 1 2 ... ! ! ! ! J J a a a a 1 0 ... Q Q Q
Each site is upstream or downstream to a neighboring site and is associated with a hydrographic zone. Stream's orientation is used as a spatio-related dimension between site (section 3.2). Hydrographic zones are used to take into account geographical granularity (section 3.3).

Experimentation

Before extracting patterns, we have to discretize data. An arbitrary discretization with 5 intervals is selected for each type of information. A description of the information that appears in patterns given in table 4 is presented in the following: ibgn: Normalized global biological index (IBGN) is a tool used to evaluate biological quality in a watershed. This biological index has a between 0 and 20 depending on the presence of some bioindicators (invertebrates).

ibgn note: is a score ranging from 0 to 5 and is based on the IBGN value.

var taxo: this data describes taxonomic variety. This is a metric corresponding to number of taxa (freshwater macroinvertebrates) collected during a sampling and is used in the IBGN computation.

We compare our approach to classical sequential pattern extraction methods (MS) and spatio-temporal patterns obtained with the approach in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF] (M ST ). Both methods are close to ours, called M ST R for Spatio-Temporal and Related. In table 3, we vary support minimum values to observe the evolution of the number of sequential patterns according to the differents methods. Table 4 Exploring hierarchical granularities and spatial relations allows the extraction of more specific and expressive patterns, not obtainable with existing methods. For instance, the pattern p xp.Orientribgn 11-15sqpvar taxo 31-40qy 4) means that frequently an IBGN value between 11 and 15 is frequently found in a neighboring site (i.e upstream or downstream) associated with a later taxonomic variety between 31 and 40. The pattern p I xp.Downstreamribgn 11- 15sqpvar taxo 31-40qy is a specialization of p and frequently finds the IBGN value between 11 and 15 in a downstream site. The pattern p P vU2wxp.Orientribgn 11-15sqpvar taxo 31-40qy means that the pattern p is frequent in the sector U2, a more specific geographic area. These patterns cannot be obtained with classical sequential patterns, e.g. xpvar taxo 31-40qy, nor the method presented in [START_REF] Salas | The pattern next door: Towards spatio-sequential pattern discovery[END_REF] best scale can be different from another. Our approach allows the presence of different hierarchical levels in the results. Finally, our approach deals with several issues: 1) considering spatial and temporal dimensions, 2) managing relations between geographical objects, and 3) exploring all granularities.

Conclusion

The method proposed in this paper tackles on mining georeferenced data and is able to consider efficiently the spatial and temporal dimensions. Our approach differs from solutions proposed in the literature, by considering both spatial relationships and granularities in a new way. The obtained patterns are semantically richer nevertheless this type of extraction leads to the exploration of a huge search space with an important amount of patterns. In the future, we wish to adapt some interestingness measures [START_REF] Tan | Selecting the right interestingness measure for association patterns[END_REF][START_REF] Geng | Interestingness measures for data mining: A survey[END_REF] to these kinds of patterns to 1) filter the patterns according to experts' needs and 2) push it in the pattern extraction process. We aim at improving the extraction time by reducing the search space, and also provide experts with the minimal and most interesting set of spatio-temporal and related patterns. An another prospect is to define some tools to help expert's navigation in results by considering ergonomic and visualization aspect.

Example 8

 8 Let D Country D S be an area dimension describing Europe's division into countries.

  DB|α); foreach dimi in DR do /* For each dimension in DR */ I θ Ð I θ Y searchIntend(θ,DB|α,dimi); I θ Ð I θ Y searchExtend(θ,DB|α,dimi); end Algorithm 3: exploreSpatialHierarchypα, θ, DB| α , D S q

Table 2 .

 2 Sequences of cityDefinition 2 (Sub-sequence) A sequence A xIS 1 IS 2 . . . IS p y is a sub-sequence of another sequence B xIS I ¨B) if p ¤ m and if there exists integers j 1 j 2 . . . j k . . . j p such as IS 1 IS j1 , IS 2 IS j2 , . . . , IS p

		City	Month Items
		Nîmes	2011/01 Humidity=Low, Sun
		Montpellier 2011/02 Sun
		Nîmes	2011/03 Heat=High
		Montpellier 2011/03 Humidity=Low, Heat=High
		Nîmes	2011/04 Heat=Low, Wind
		Orange	2011/04 Rain
		Orange	2011/06 Rain, Wind
			Table 1. Database
	City	Sequence
	Nîmes Montpellier x(Sun)(Humidity=Low Heat=High)y x(Humidity=Low Sun)(Heat=High)(Heat=Low Wind)y Orange x(Rain)(Rain Wind)y
	Extracting knowledge from sequences search frequent sub-sequences, named
	as sequential patterns. Several algorithms have been proposed for sequential
	pattern mining [7-12].
			1 IS I 2 . . . IS I
	IS jp .
	Example 1 Consider the sequences presented in table

m y (A

Table 3 .

 3 presents an example of patterns extracted for each method. Number of extracted patterns according to minimum support

	MS MST MST R
	0.5 1	4	4
	0.4 4	12	12
	0.3 22	60	64
	0.2 75 186	233
	0.1 180 445	1882

Table 4 .

 4 , e.g. xp.ribgn 11-15sqpvar taxo 31-40qy. Experts often have difficulties to determine the best scale to obtain the best observations, and for each parameter, the Downstreamribgn 11-15sqpvar taxo 31-40qy 0.051 vU1wxp.Orientribgn 11-15sqpvar taxo 31-40qy 0.054 vU2wxp.Orientribgn 11-15sqpvar taxo 31-40qy 0.073 xp.Orientribgn note 3sqpvar taxo 31-40qy 0.056 vU1wxp.Orientribgn note 3sqpvar taxo 31-40qy 0.051 Patterns according to different methods

	Method	Sequence	Support
	MS MST MST R	xpvar taxo 31-40qy xp.ribgn 11-15sqpvar taxo 31-40qy xp.ribgn note 3sqpvar taxo 31-40qy xp.Orientribgn 11-15sqpvar taxo 31-40qy xp.	0.404 0.089 0.056 0.089

Acknowledgments

Thomas Lampert is gratefully acknowledged for helpful comments on the manuscript. This work was partly funded by french contract ANR11 MONU14.


Simon Struck
email: simon.struck@ovgu.de

Michael Lipaczewski
email: michael.lipaczewski@ovgu.de

Frank Ortmeier
email: frank.ortmeier@ovgu.de

Matthias Güdemann
email: matthias.gudemann@inria.fr

Multi-Objective Optimization of Formal Specifications

Even in the domain of safety critical systems, safety and reliability are not the only goals and a developing engineer is faced with the problem to find good compromises wrt. other antagonistic objectives, in particular economic aspects of a system. Thus there does not exist a single optimal design variant of a system but only compromises each "best" in its own rights . With the rising complexity, especially of cyber-physical systems, the process of manually finding best compromises becomes even more difficult.

To cope with this problem, we propose a model-based optimization approach which uses quantitative model-based safety analysis. While the general approach is tool-independent, we implement it technically by introducing well defined variation points to a formal system model. These allow enough variability to cover whole families of systems while still being rigorous enough for formal analysis. From the specification of this family of system variants and a set of objective functions, we compute Pareto optimal sets, which represent best compromises.

In this paper we present a framework which allows for optimization of arbitrary quantitative goal functions, in particular probabilistic temporal logic properties used for model-based safety analysis. Nevertheless, the approach itself is well applicable to other domains.

I. INTRODUCTION

The term "cyber-physical system" describes a system with a tight coupling between the computational units of a system and its physical components. Computational units are all kinds of programmable devices such as micro controllers, micro computers, embedded systems and computer systems. Physical components comprise all kind of sensors and actors, but also every mechanical interconnection between them. The overall system behavior cannot be determined by analyzing the program code or the physical components in isolation; only a combined analysis of both explains the system behavior.

To cope with the safety requirements of complex reactive systems, model-based safety analysis methods were developed [START_REF] Abdulla | Designing safe, reliable systems using SCADE[END_REF], [START_REF] Bozzano | Improving system reliability via model checking: the FSAP/NuSMV-SA safety analysis platform[END_REF]. These methods rely on a formal model which consists of the software control, i.e., the behavior of the computational units, and its environment, i.e., all the physical components interconnected with the computational units. This model is then analyzed for safety aspects using deduction, most often based on automatic model checking techniques.

However, for most cyber-physical systems, safety is not the only important aspect, even in safety critical domains. Other aspects are for example system costs or performance. A system with very high dependability may simply be too expensive to build, or its performance might be below expectations or even requirements. On the other hand, it will often be a goal to build a system with the highest possible performance and lowest cost, but at the same time with an appropriate safety to allow for necessary certifications. Because of the often antagonistic nature of these different aspects, there is typically no single optimal possibility to develop a cyber-physical system. In general there exist different variants, where each is a best compromise in the sense that only worsening one aspect can augment another. The identification of such compromises is no trivial task. It is traditionally performed manually by engineers during system development.

In this paper we propose a systematic model-based optimization approach to find the set of best compromises for a cyber-physical system. We implemented the approach using a formal system modeling language (SAML) and multi-objective optimization algorithms to identify best system specifications.

The antagonistic nature of reliability, availability, maintainability and safety (RAMS) is also context of standards like IEC62278/EN50126. In general, this can be seen as a multiobjective optimization problem, i.e. the identification of best compromises among the four different objectives. Thus our approach is capable of targeting requirements according to RAMS.

The remainder of the paper provides an introduction to multi-objective optimization in Section II. Section III provides a brief introduction of SAML and introduces our formal analysis techniques. We describe our optimization approach of formal specifications in Section IV and provide an evaluation of our approach on a case study in Section V. Finally, Section VI provides some related approaches and Section VII gives an outlook on future work.

II. MULTI-OBJECTIVE OPTIMIZATION

To optimize a system design wrt. different goals and to identify the best compromises, we use a very abstract view on a system. We consider the possible variations of a system as free parameters of an optimization problem and the different characteristics of the system as its objective functions. More precisely, we consider it to be a multi-objective optimization problem (MOP). Every possible configuration of the parameters describes a specific system variant for which the objective functions are evaluated. One of these is the probabilistic safety analysis, which consists of the computation of the hazard occurrence probability.

In general a MOP is described as follows: All k parameters of a system are denoted as a vector of decision variables x = {x 1 , . . . , x k } ∈ X. The l objectives of the optimization are functions i = 0 . . . l, f i : X → R. All l objectives together forming the vector-valued function f = (f 1 (x), . . . , f l (x))

T , the MOP is defined as the minimization of the vector-valued objective function:

min (f (x)) = min    f 1 (x) . . . f l (x)    (1)
For such vector functions, it is in general impossible to define a total order, in particular if the objective functions f i are antagonistic to each other. Consider for example a system that moves goods from one point to another. While designing such a system, the engineer needs to determine the maximum allowed speed. In this example, there might be two objectives functions that need to be considered: The traveling time and the systems hazard probability. Obviously with a higher speed, the traveling time decreases and the hazard probability increases. It is clear that there is no single best solution considering the minimization of the hazard probability and the traveling time at the same time.

To cope with antagonistic objective functions we rely on the Pareto dominance criterion. It is a partial order which allows for optimality comparison of different solutions. Further details are for example in [START_REF] Miettinen | Some methods for nonlinear multi-objective optimization[END_REF].

Definition 1 (Pareto Dominance):

A vector u = (u 1 , . . . , u n ) ∈ R n dominates another vector v = (v 1 , . . . , v n ) ∈ R n if and only if: ∀i ∈ {1, . . . , n} : u i ≤ v i ∧ ∃j ∈ {1, . . . , n} : u j < v i . (2)
If u dominates v this is also denoted as u ≺ v.

Given two vectors u and v, where u ≺ v means that all elements in u are less than or equal to the corresponding elements in v and at least one element in u is strictly smaller than the corresponding element in v. In other words: u is for at least one goal function strictly better than v and not worse in all others.

Applied to the example of the transportation system, the Pareto dominance considers a solution with a short traveling time and a high hazard probability to be equally good as a system with a long traveling time and a low hazard probability. Yet a solution with a high traveling time and a high hazard probability is not optimal. This fits the idea of the best compromises. All parameter settings x that lead to Pareto optimal solutions form the Pareto set: Definition 2 (Pareto Set): For a given multi-objective optimization problem, f (x) with x ∈ X, the Pareto Set, P * , is defined as:

P * := {x ∈ X | ¬∃x ′ ∈ X : f (x ′ ) ≺ f (x)}. (3)
The set of images of the elements of the Pareto set defines the Pareto front, as the set of solutions that do not dominate each other, and are not dominated by any other available solution:

Definition 3 (Pareto Front): For a given multi-objective optimization problem, f (x) with x ∈ X, and its Pareto Set P * the Pareto front PF * is defined as:

PF * := {u = f (x) | x ∈ P * } (4)
Finding the Pareto front is a non-trivial problem even for analytic functions. In the past we did some experiments to approximate the Pareto set using evolutionary algorithms [START_REF] Güdemann | Model-Based Multi-Objective Safety Optimization[END_REF]. They are well applicable and are widely used for MOP. Besides sophisticated search algorithms, more pragmatic approaches are also feasible. In Section V we use brute force and nondominated sorting to solve a MOP with a sufficiently small search space.

All elements in the Pareto optimal set are considered as equally good solutions. To select a single one for the realization of one system variant, further criteria are required. These could either be an expert's choice, or further heuristics. However, these are very domain specific and thus not discussed in this paper.

III. MODEL-BASED SAFETY ANALYSIS

Model-based safety analysis forms the basis for the proposed optimization approach. It starts from a formal description of a model, its physical environment and the possible failure modes and uses formal deduction techniques to analyze the system.

A. Specification

We express models in the safety analysis modeling language (SAML). Here, we only cover its syntax and semantics as far as it is required to understand the optimization approach. A complete definition is presented in [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF], [START_REF] Güdemann | Towards Model-driven Safety Analysis[END_REF]. It is derived from the PRISM language [START_REF] Norman | The PRISM language -semantics[END_REF], allows for modeling of finite state automata and uses a discrete time model. The state automata are executed in a synchronous parallel way, the transitions allow for the expression of combined probabilistic and nondeterministic behavior. Every automaton description is included in a module that consists of one or more state variables and a set of update rules. Figure 1 shows a brief example. The module myMod contains one state variable named state, defined for the interval of naturals from one to three and the initial value one.

Transitions are specified in update rules. Every update rule has an activation condition and one or more non-deterministic choices. The activation condition is a propositional logic formula. It denotes in which system states the update rule is active. In the example there is one update rule for each of the three possible states. If the activation condition evaluates to true, one of the following non-deterministic choices (denoted by the keyword choice) is selected. In the example, the first update rule consists of one choice and a trivial probability distribution. The third update rule consists of two choices. Within every non-deterministic choice a discrete probability distribution is given. The second update rule in the example shows a non-trivial probability distribution. The constant (prob) denotes the probability of reaching one of the two possible successor states. For a variable v, v ′ marks the value in the next time step, i.e., the successor state. If the module uses more than one state variable, the assignment is done in parallel.

Constants are used as a generic concept to use named aliases for certain numbers. Constants may be of floating point and integer types. Additional to the usage in the example as probabilities (prob) they may also be used in the propositional expressions of the activation conditions and in next state assignments.

Besides constants and modules, SAML supports a third construct called formulas which are propositional logic expressions over the state variables. In Summary, a SAML model consists of a set of constants, a set of formulas and a set of modules:

Definition 4 (SAML Model): A SAML model is syntactically defined as a tuple

S = (C, F , M) , (5)
where C is a set of constants, F is a set of formulas and M is a set of modules.

The semantics of a SAML model with one single module is a Markov decision process (MDP). An MDP consists of a finite set of states, one initial state and a non-empty set of probability distributions over the successor states. A labeling function assigns a set of atomic Boolean propositions to each state. For a detailed definition see, e.g. [START_REF] De Alfaro | Model checking discounted temporal properties[END_REF].

If a model contains more than one module and thus more than one state automaton, these are combined into one single module by parallel composition of all modules, creating the product automaton of all synchronous parallel automata.

Informally, a SAML model with a single module and a single state variable maps to an MDP in the following way: All possible values of the state variable form the set of states. The initial value (stated with the init keyword) denotes the initial state. The update rules of the module contain the probability distributions over the successor states. Multiple distributions can be assigned to every state. The exact semantics of SAML models can be found in [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF].

The MDP for the SAML model listed in Figure 1 is depicted in Figure 2. Here every value of the single state variable represents one state. The transitions are labeled with k : p where k is the index of the non-deterministic choice that is available in the originating state and p is the probability of this transition in the selected choice. A valid SAML model must fulfill the following requirements: (1) All the probabilities in a distribution must sum to one. This ensures that the distributions are valid probability distributions. [START_REF] Bozzano | Improving system reliability via model checking: the FSAP/NuSMV-SA safety analysis platform[END_REF] The activation conditions in one module must be pairwise exclusive. This assures that normalization of the probabilities is not necessary, their intended value does not change and also that there is no unwanted non-determinism in the model. [START_REF] Miettinen | Some methods for nonlinear multi-objective optimization[END_REF] The transition relation must be total. This assures that every modeled deadlock state is also deliberately integrated.

B. Analysis

For an actual analysis, SAML models are transformed into the input languages of state of the art model checkers using provably sound model transformations [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF]. Currently implemented are NuSMV1 and PRISM 2 . In this section we provide a brief introduction to model based safety analysis of SAML models. A complete explanation of model-based safety analysis is out of scope of this paper.

NuSMV is a symbolic model checker for Kripke structures. It performs qualitative analysis, for which all the probabilistic information in the SAML model is stripped and replaced with non-determinism. The qualitative safety analysis points out worst case scenarios. The properties are specified in computational tree logic (CTL) [START_REF] Clarke | Model Checking[END_REF]. If the model checker finds a violation of a property in the given model, it presents a sequence of states (starting from the initial state) that demonstrates the violation.

The results of the qualitative analysis are not directly relevant for the optimization procedure, but can identify infeasible system variants. Still, they can provide very valuable information about the behavior of the model to the system engineer. Counter examples can be produced which show a possible sequence of states, for further analysis by a human analyst or engineer.

Qualitative analysis is very useful for safety analysis. The deductive cause consequence analysis (DCCA) is a qualitative model-based safety analysis approach [START_REF] Güdemann | Computing ordered minimal critical sets[END_REF]. For DCCA, the analyst must identify unique failure modes of system components. The model is then extended with a set of automata which represents the possible occurrence for every failure mode. The functional part of the model is extended in such a way that every component (module) has an appropriate failure behavior if the corresponding failure automaton is in the failure state. It is important to notice that the failure extended model is path inclusion equivalent to the model without failure extension [START_REF] Ortmeier | Formal failure models[END_REF]. Based on the failure extended model and a specification of an hazard (a set of states which are considered as hazardous), the DCCA automatically calculates all possible failure combinations that can lead to a hazardous system state. In general, qualitative analysis is used in our optimization approach for identifying feasible solutions which satisfy a set of qualitative properties.

The quantitative safety analysis is performed using a probabilistic model checker. It is of fundamental importance for the optimization. Quantitative properties exploit the probabilistic information in the SAML models. Their result is not only a simple yes or no but a probability stating how likely it is that a trace of the model is chosen on which a certain temporal logic property holds . The specification of properties is expressed in PCTL probabilistic temporal logic [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF]. For safety analysis we use the probabilistic deductive cause consequence analysis (pDCCA) [START_REF] Güdemann | Probabilistic model-based safety analysis[END_REF], [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF]. It computes the occurrence probability of the system hazard, if the occurrence probabilities of all possible failures are specified. The approach allows for a sensible combination of per-time failure modes specified by failure rates and per-demand failure modes specified by failure probabilities.

IV. OPTIMIZATION OF A SPECIFICATION

A. Syntax and Semantics

Qualitative and quantitative analysis can be performed on SAML models. Using this, a system can be optimized for multiple goals by extending the specification with certain variation points. We now propose a modifier keyword (param) to SAML, used in the context of constants and modules. We named the new language optimizable-SAML (OSAML).

• Param Module: describes a set of interchangeable modules. Each of the modules must have the same signature in terms of the state space (i.e. Alternative modules must have the same state variables). However, the state transitions may differ. Thus variants of similar behavior can be expressed. • Param Constant: describes constants in a model that may be changed. Throughout one analysis this value is constant, but another analysis may be performed with different values for the constant. The idea behind param module is to model parts with a similar purpose but different architecture or technique. Consider the configuration of a car where several different tachometers are available; all do measure the speed, but one does this by measuring the rotational speed of the wheels and another is based on satellite navigation. Both will differ in accuracy, reliability and costs. In a formal model such a component might be modeled with a state variable representing the current reading. The rotational measurement can immediately present coarse results, while the satellite-based measurement presents very accurate results but with a slower update interval. In both cases the state variables are the same, but the update rules differ.

A very minimalistic example with one param constant and one param module is depicted in Figure 3. The first variant of the param module (variantA) is the same as the module depicted in Figure 1. The second variant (variantB) possesses the same definition for the state variable (state : [1..3]), but has other update rules. The constant used in the example from Figure 1 is now also turned into a param constant. This means, that it specifies no specific value but an interval of allowed/valid values. An OSAML model is a tuple of sets of all the elements. In addition to the SAML model there are two more sets for the param modules and the param constants:

Definition 5: An OSAML model is syntactically defined as a tuple

OS = C O , PC O , F O , M O , PM O , (6)
where C O is a set of constants, F O is a set of formulae, M O is a set of modules, PC O is a set of sets of param constants and PM O is a set of param modules.

Here PC O is a finite set of param constants. Every param constant may cover an infinite set of values. Similarly PM O is a finite set of param modules where every param module contains a finite number of (normal) modules.

On the semantic layer, an OSAML model is a family of SAML models (which matches a family of MDPs). Based on the assumption that a family can contain good and bad candidates there is no point in analyzing a whole family. For a meaningful analysis one specific variant for each param module and a certain value for every param constant must be chosen. Thus only members of the model family are analyzed. We call the process of selecting one specific SAML model from a family the instantiation of an OSAML model.

Definition 6 (Instance of an OSAML model): An SAML model S is an instance of a OSAML model OS if and only if all of the following rules apply:

1) For every param constant in the OSAML model OS a constant (with the same name) exists in the SAML model S.

2) The types of the param constants in the OSAML model OS must match with the corresponding constants in the SAML model S.

3) The values of the constants corresponding to the param constants must fit in the domain of the param constants. 4) For every param module in the OSAML model OS one module exists in the SAML model S that equals one module in the corresponding param module group. 5) Every constant, formula and module in the SAML model S is also in the OSAML model OS. 6) The SAML model S contains no additional constants, formulas and modules except the above mentioned. We omit a complete formal definition because it is rather complicated to express, while the meaning is pretty straight forward. Instantiating an OSAML model is done in the following way: In a first step, all modules, constants and formulas are transferred from the OSAML model to an empty SAML model that is going to become the instance. In a second step, every param constant is turned into an appropriate constant of the same type and with a value from the domain of the param constant. The specific values of the generated constants depends on the desired instance. The newly generated constants are added to the list of constants in the SAML model. In the third step, one module is picked out of every param module and added to the list of modules in the SAML model. Obviously, now every instance of an OSAML model can be analyzed as every other SAML model.

We extended the existing SAML grammar [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF] with rules for the param keyword and used the ANTLR [START_REF] Parr | The Definitive ANTLR Reference: Building Domain-Specific Languages, ser. Pragmatic Programmers[END_REF] parser generator for the creation of a OSAML parser. After parsing all constants, formulas, modules, param constants and param modules they are stored in an abstract syntax tree consisting of instances as Java objects. Given the abstract syntax tree the instantiation of a OSAML model is then a trivial procedure.

B. Optimization Framework

We extended the existing implementation of the SAML framework towards optimization. The software concept is depicted in Figure 4. An OSAML model is covered by a SamlSystemCluster instance. The generateInstance method performs the instantiation. The resulting SAML model is represented by a SamlSystem object.

All objectives for the optimization must implement the IObjective interface. In Figure 4 this is exemplary shown for the PRISM model checker based safety analysis. This class (PrismObjective) contains the model transformation and also handles the invocation of the model checker. ´«use´»

´«use´»

´«create´»

´«use´» Fig. 4. Simplified software architecture of the SAML optimizer All optimization algorithms must implement the IOptimizer interface. Currently there are implementations of a brute force enumeration algorithm and the fast non-dominated sorting genetic algorithm (NSGA2) [START_REF] Deb | A fast and elitist multi-objective genetic algorithm: NSGA-II[END_REF] available. The optimizer implementations use the SamlSystemCluster to create instances of the OSAML model and the IObjective interface to evaluate the optimization objectives. In general, the architecture of the framework provides two ways to specify objectives for the optimization.

• Use one or more CTL/PCTL formulas with the (already existing) PRISM or NuSMV objectives. • Generic objectives can be expressed in Java code that implements the IObjective interface. These objectives can access both the SAML instance and the parameter configuration used to generate the instance. Both methods can be used in the same optimization problem. It is therefore possible to optimize two ore more PCTL properties and several Java based objectives at the same time.

V. EVALUATION

We evaluated the approach on a case study of a redundant data processing system. The goal was to identify optimal design solutions by minimizing the hazard probability and the system costs at the same time.

A. Case Study

The following case study from model-based safety analysis is taken from literature. It was first presented in [START_REF] Walker | Compositional Temporal Fault Tree Analysis[END_REF].

The case study consists of the measuring of an input signal, its computational processing, and the generation of an output signal. For higher reliability certain components are redundant. A block diagram is depicted in Figure 5. The two sensors S1 and S2 measure an input signal. The sensor values are then processed by the arithmetical units A1 and A2. The first arithmetical unit uses both sensor signals and the second one only uses one sensor signal. The second arithmetical unit is disabled by default. A monitor M observes A1 and activates the second algorithmic unit if the first one fails. The output unit O selects the proper signal from the two redundant arithmetical units. A quantitative and qualitative analysis of the model may be found in [START_REF] Güdemann | Computing ordered minimal critical sets[END_REF] and [START_REF] Güdemann | Probabilistic model-based safety analysis[END_REF]. For optimization, we introduced three parameters in addition to the original case study. This turns the specification into a family of specifications. All parameters are expressed in the OSAML model with param constants and param modules. The parameters are as follows:

• Instead of using a separate sensor S1 and arithmetical unit A1, one unit that combines the functionality of the sensor and the arithmetical unit can be used. In Figure 5 this is indicated by the dashed box around S1 and A1. • Analogously to S1 and A1 also S2 and A2 may be combined into one unit. In the figure this is indicated with the second dashed box. • For the monitor there are 5 different realizations available, which differ in the probability:. 1 • 10 -4 , 5 • 10 -4 , 1 • 10 -5 , 5 • 10 -5 , 1 • 10 -6 and 5 • 10 -6 .

The failure of all other components is modeled with a probability of 10 -5 . We did not assign a distinct time model to our formal state-based model. Thus all stated probabilities are per-step probabilities. If the analysis results shall be evaluated for a real system, a sampling interval ∆t must be specified. The relation between failure rates and the per-step probability is defined in [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF].

The first two parameters are realized with the param module feature in OSAML. The third parameter can be expressed as param constant or also as param module. A param constant would allow the setting of arbitrary values between the lower and upper bounds. Thus we used a param module to express the logarithmic scaling of the failure probability parameter.

B. Pareto Optimization of the Case-Study

For multi-objective optimization of the case study, we introduced two different objective functions. The first one is the occurrence probability of the hazard within k time-steps, where k is also referred as mission time. This hazard is defined as the state when the system is not able to generate a proper output signal. With H describing the hazardous state (i.e., when the system does not report a proper output signal) the failure probability objective using pDCCA is defined as PCTL formula 3 :

P max k (H) := P max =? [trueU ≤k H] (7)
This objective is evaluated with the PRISM model checker. For the optimization we assumed a mission time k = 360 steps.

The second objective calculates the costs of the system. It is based on the assumption that a system with less components is cheaper and that components with lower failure probability are more expensive. The combination of a sensor and an arithmetical unit into one component is assumed to cost 1, whereas the separation in two components costs 2. For the Monitor module the costs are 1, 2, 4, 8 and 16 in sequence of decreasing failure probability. The overall system costs are then defined as the sum of the costs of all chosen system components. This objective function is implemented directly as Java method. This leads to two binary parameters and one parameter with up to five different values. All in all this gives a total of 20 different parameter configurations. The interesting question is, which of these configurations lead to Pareto optimal systems.

Due to the fact that the search space of the problem rather small we decided to solve the problem with brute force via the enumeration of the design space. However, the current implementation also contains an evolutionary algorithm for dealing with larger or infinite families.

The Pareto front for the example is depicted in Figure 6, the failure probability on the x axis and the system costs on the y axis. All ten exact parameter configurations of the Pareto optimal solutions are listed in Table I. For the example, only half of all available parameter configurations lead to optimal systems. This means that for every non optimal parameter configuration a different configuration exists, that is both cheaper and has a lower hazard probability at the same time. By further analyzing the results in Figure 6 it seems like there are dominated points in the front. However, Table I clarifies that they are all non-dominated. It turns out that the affected pairs of design variants ({4, 5}, {7, 8} and {9, 10}) only differ in the realization of S2 and A2. This leads only to a very small decrease in the hazard probability while causing a strong impact on the costs. All other parameter configurations differ rather strongly in both, the hazard probability and the system costs. This means that in theory, all these solutions must be considered equally good. In reality of course, one of the variants must be chosen for the implementation.

If there is only a minor impact on one objective while there is a strong impact on the other objective it is straight forward to pick the cheaper one. To choose a variant when a stronger compromise is required, further strategies are needed. One way is to prioritize the objective functions, e.g., setting a threshold for the hazard probability and then taking the system with the lowest cost which is below that threshold. Another approach is to identify regions where the increase of one objective function leads to a rather large deterioration in at least one other [START_REF] Branke | Finding knees in multiobjective optimization[END_REF]. According to this rule, solution number 5 seems to be a good choice. There are of course other methods and heuristics to chose a single design variant from a Pareto set, often dependent on the nature of the problem domain.

VI. RELATED WORK

A first safety optimization approach which uses hazard probabilities as objectives was presented in [START_REF] Ortmeier | Safety optimization: A combination of fault tree analysis and optimization techniques[END_REF]. This aposteriori approach used an analytic mathematical model of all objectives for the optimization. Even though this method is computationally efficient, it relies heavily on stochastic independence which is not a very realistic assumption in more complex models.

The underlying idea of the SAML based multi-objective optimization was presented in [START_REF] Güdemann | Model-Based Multi-Objective Safety Optimization[END_REF]. In this paper the authors applied an evolutionary algorithm on a parametric SAML specification. To speed up the optimization process, artificial neural networks where proposed. However, it was on a very prototypical level without a grammar nor specific semantics of the parameters used in the specification. OSAML continues this initial idea, and introduces proper syntax and semantic as well as a proper implementation and a framework for the objective functions.

Another optimization approach of formal specifications is presented in [START_REF] Meedeniya | Architecturedriven reliability and energy optimization for complex embedded systems[END_REF]. In this case Markov reward models of embedded systems are optimized along reliability and energy consumption objectives. However, this approach is limited to quantitative analysis. By using SAML as specification language we can perform qualitative and quantitative analysis on the same model.

A completely different approach of analyzing parametric Markov models is presented in [START_REF] Hahn | Param: A model checker for parametric markov models[END_REF]. The parametric Markov model and a PCTL property is transformed into a polynomial formula, that states the dependency between the parameter values and the analysis result. However, translating a parametric model into a polynomial formula is a complex problem. Even if it is solvable for complex case studies the complexity still remains an issue. Also the approach presented by Hahn et al. does not fully support Markov decision processes (as SAML does).

Design space exploration (DSE) is another important topic for identifying ideal system configurations. DSE based approaches are presented for example in [START_REF] Hegedus | A model-driven framework for guided design space exploration[END_REF] and [START_REF] Knorreck | Formal system-level design space exploration[END_REF]. These approaches use either UML or proprietary modeling languages for system description and a explanatory algorithm to seek optimal solutions. In contrast to our Approach, DSE uses a set of variation rules to alter the model. The proposed modeling formalisms support the expression of complex and large structural models, but miss the support for the verification of functional aspects and non-deterministic behavior. From a theoretical point of view our approach should be able to cover the objectives aimed by the DSE approaches.

VII. CONCLUSION AND OUTLOOK

We proposed an approach to optimize formal specifications wrt. objective function expressed either as PCTL properties or directly via Java methods. As the objectives can be antagonistic, there exists no single optimal solution but only best compromises. We defined optimality in terms of minimization of a multi-objective optimization problem and used the Pareto dominance criterion to compare different design variants.

On the specification layer we used the SAML specification language. To express whole families of specifications in one model we extended the language with the new keyword param. The new version of the language is named optimizable safety analysis modeling language (OSAML). Parameters of a specification family can be the variation of the values of constants or the specification of alternative modules.

On the analysis layer we rely on sound model transformations from SAML into the specification languages of NuSMV and PRISM [START_REF] Güdemann | A framework for qualitative and quantitative model-based safety analysis[END_REF]. This allows for using PCTL properties as objectives for the optimization process. Our software architecture allows the implementation of further Objectives in the Java programming language.

We evaluated our approach on a case study from model based safety analysis. Our optimization along costs and hazard probability identified 10 Pareto optimal candidates. Despite the rather small case study in this paper we also implemented a evolutionary algorithm to efficiently cover a large search space in more complex models.

Particularly motivated by the fact that probabilistic model checking is time consuming, we are going to further investigate optimization algorithms. Brute force enumeration, as used for the case study in this paper, is only feasible for very small search spaces. Nevertheless we believe that often a reduction of the search space is possible by exploiting meta information from the application domain of the optimization problem. Very often the parameters represent physical quantities of components in the system. Due to tolerances these cannot be manufactured with an arbitrary precision. Thus it is not necessary to operate with a finely grained search space, if this does not map to the parameters of the real system. A possible approach could be the specification of logarithmically scaled parameters (similar to the failure probability parameter of the monitor module in the case study).

A second approach for the reduction of model checking time is to integrate stochastic estimation algorithms into the optimization procedure. We are currently busy with the evaluation of artificial neural network and simulation based estimation of PCTL formulas. For artificial neural networks the basic idea (as laid out in [START_REF] Güdemann | Model-Based Multi-Objective Safety Optimization[END_REF]) is the usage of previous function evaluations to train the estimation network. It can then be used to estimate further function evaluations. Instead creating an independent estimation model we are evaluating the usage of the verification models in terms of Monte Carlo simulation as estimator. We expect a relevant speed-up due to the usage of estimation techniques.

Fig. 1 .

 1 Fig. 1. A minimalistic SAML model with only one module

Fig. 2 .

 2 Fig.2. The MDP for the SAML model in Figure1

Fig. 3 .

 3 Fig. 3. Minimalistic OSAML model

Fig. 5 .

 5 Fig. 5. A System with Redundant Components

Fig. 6 .

 6 Fig. 6. The Pareto Front

http://nusmv.fbk.eu/

http://www.prismmodelchecker.org/

The =? in the formula belongs to PCTL syntax and applies to the quantitative aspect of PCTL. It denotes that the formula evaluates to the occurrence probability of the subsequent CTL formula.

ACKNOWLEDGMENTS

Michael Lipaczewski is sponsored by the Deutschen Ministerium für Bildung und Forschung in the ViERforES project (BMBF, project-Nr.: 01IM08003C).

Simon Struck is sponsored by the German Research Foundation (DFG) within the ProMoSA project.