Nikolaos Kekatos

Marcelo Forets

Goran Frehse

Constructing Verification Models of Nonlinear Simulink Systems via Syntactic Hybridization*

In this paper, we present a methodology that facilitates the integration of formal verification techniques into model-based design. The focus is on set-based reachability analysis and on control systems that are described by hybrid dynamics and nonlinear components. Starting with a standard simulation model, e.g., in MATLAB/Simulink, we transform it into an equivalent verification model, formally a network of hybrid automata, in the SX format used by several reachability tools. A major obstacle here is that highly scalable reachability algorithms and tools exist for piecewise affine (PWA) dynamical models, but not for nonlinear ones. To obtain PWA overapproximations of nonlinear dynamics, we use an abstraction method known as hybridization. It consists in partitioning the state-space into a set of domains, and for each domain, approximating the nonlinear dynamics by simpler ones plus nondeterministic inputs to account for the abstraction error. Existing hybridization procedures operate on the composed (flattened) system, so the number of partitions is exponential in the number of variables. This quickly leads to intractably large models, even for small systems. To mitigate this problem, we decompose the original dynamics and carry out the statespace partitioning and PWA approximation on the components. The number of partitions in each PWA component is at most quadratic in the abstraction error so that an explosion in the number of partitions is largely avoided. Since the SX format can handle templates, several components may share the same abstraction. The result is a highly compact model that retains the modular structure of the original simulation model. If only a small subset of the partitions is reachable, the bottleneck of having excessively large PWA models can be avoided by composing the model on-the-fly during the reachability analysis.

I. INTRODUCTION

In model-based design (MBD), the plant and its controller are designed based on a model, typically within a simulation environment like MATLAB/Simulink. Any kind of nondeterminism in the system, like disturbances, measurement noise, parameter uncertainties, user input, or operating conditions, may have adverse effects on the performance. These effects can be difficult to predict during the design step. Therefore, the system is typically tested by simulating a large number of trajectories, each with a different choice for the nondeterministic quantities, and checking whether they satisfy the requirements. This process is generally incomplete since the number of different choices is prohibitively large or even infinite. Therefore, it can be hard to say with high *This work was partially supported by the European Commission under grant 643921 (UnCoVerCPS) and by the Metro Grenoble through the project NANO2017.

1 Nikolaos Kekatos, Marcelo Forets and Goran Frehse are with VERIMAG Laboratory, University Grenoble Alpes, Bâtiment IMAG, 700 Avenue Centrale, 38400 Saint-Martin-d'Hères, France. {Nikolaos.Kekatos, Marcelo.Forets-Irurtia, Goran.Frehse}@univ-grenoble-alpes.fr confidence whether a requirement is truly satisfied under all circumstances. Formal verification attempts to guarantee that requirements are satisfied through a rigorous mathematical analysis of the system. A widely used verification technique is set-based reachability analysis, which exhaustively simulates families of trajectories using geometric operations on sets.

There are two main obstacles to applying reachability analysis in MBD. First, the simulation model needs to be converted to a suitable formal model, such as a hybrid automaton. Second, the model must be amenable to existing reachability algorithms, in particular in terms of scale. Highly scalable algorithms are known for piecewise-affine (PWA) dynamical systems, but not for more complex nonlinearities. While a large class of nonlinearities can be approximated arbitrarily well by a PWA system, the resulting models can be very large, again running into scalability problems.

In this paper, we propose an approach to transform a simulation model into a compact, i.e., relatively small, verification model with PWA dynamics. To achieve this, we decompose the nonlinear system and perform the transformation component-wise. The resulting model can be fed to the verification tool SpaceEx [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF], or translated into formats for other verification tools using the HyST tool [START_REF] Bak | HYST: a source transformation and translation tool for hybrid automaton models[END_REF]. Since SpaceEx composes the model on-the-fly during the analysis, only the reachable partitions of the PWA approximations are instantiated. Figure 1 illustrates the difference between the traditional hybridization methods and the proposed one, for the case where the original nonlinear model is described in Simulink. Classical hybridization techniques that rely on state-space partitioning [START_REF] Asarin | Hybridization methods for the analysis of nonlinear systems[END_REF] create a PWA model with O(1/ n ) locations, where is the mesh size, and n is the dimension of state-space. During the reachability analysis, O(T /δ) locations are visited, where δ is the minimum dwell time and T the global time horizon. On the contrary, with syntactic hybridization, we get m PWA components, where m is the number of nonlinearities, and the total number of locations is O(m/ 2 ). The non-reachable locations need not be instantiated.

The rest of this paper is organized as follows. In Section II, we present the related literature. In Section III, we describe the steps of our proposed methodology and apply them to a simple Simulink model example. In Section IV, we introduce the compositional syntactic hybridization and apply it in an industrial case study. Finally, we draw conclusions and perspectives in Section V.

II. RELATED WORK

The translation of Simulink models into modeling formalisms, for which formal verification tools can be applied, has attracted considerable interest; a comprehensive survey can be found at [START_REF] Zhan | Formal Verification of Simulink/Stateflow Diagrams: A Deductive Approach[END_REF]. There exist translators from Simulink to Lustre [START_REF] Tripakis | Translating discretetime Simulink to Lustre[END_REF], NuSMV model checker [START_REF] Meenakshi | Tool for translating Simulink models into input language of a model checker[END_REF] and BIP [START_REF] Sfyrla | Compositional translation of Simulink models into synchronous BIP[END_REF]. However, all these tools apply discrete verification and do not consider continuous-time models. Filipovikj et al. [START_REF] Filipovikj | Simulink to UPPAAL statistical model checker: Analyzing automotive industrial systems[END_REF] transformed Simulink models into the input language of UP-PAAL Statistical Model checker. UPPAAL supports hybrid automata. However, it either restricts their continuous parts to simple dynamics or applies the Euler integration method. Its inaccurate integration results are therefore not conservative. Zuiliani et al. [START_REF] Zuliani | Bayesian statistical model checking with application to Stateflow/Simulink verification[END_REF] presented a statistical model checking approach that is applicable on Simulink/Stateflow models. Stanley Bak et al. introduced a translation process from Simulink/Stateflow to hybrid automata in [START_REF] Manamcheri | A step towards verification and synthesis from Simulink/Stateflow models[END_REF]. Both these works focus entirely on Stateflow diagrams and require the transformation of the Simulink model into a Stateflow one. However, this is not feasible for most large-scale systems designed with Simulink.

The translation of a Simulink/Stateflow model to a hybrid automaton is supported by the tools HyLink [START_REF] Manamcheri | A step towards verification and synthesis from Simulink/Stateflow models[END_REF] and GreAT [START_REF] Agrawal | Semantic translation of Simulink/Stateflow models to hybrid automata using graph transformations[END_REF]. However, these tools do not allow hierarchical modeling and can be applied to a small subset of Simulink blocks. Recently, Minopoli and Frehse presented SL2SX, a semiautomated tool for translation of Simulink models into hybrid automata [START_REF] Minopoli | SL2SX translator: From Simulink to SpaceEx models[END_REF]. The translator supports a large number of Simulink blocks, but is restricted by SpaceEx limitation to handle piecewise constant and affine dynamics. As a result, the user should analyze the missing blocks (unsupported or nonlinear) and decide how to replace or approximate them.

It is well-known that non-compositional methods for PWA approximations are not computationally efficient for complex systems, since an acceptable accuracy requires a very large number of pieces in the piecewise-affine approximation. Very recently, [START_REF] Deshmukh | Piecewise-affine approximations for a powertrain control verification benchmark[END_REF] presented an experimental comparison of a compositional approach, similar to that presented in this paper (called nested approximation there), against a simplex-partitioning PWA hybridization, showing that the former scales much better than the latter for increasing demands on precision. The compositional PWA-approximation is presented informally, and the paper does not discuss the implications in terms of the model size, nor preserving this compositionality in the generated model. The complexity of the approximation can further be reduced by focusing on a set of reference trajectories, as done in [START_REF] Bak | Scalable static hybridization methods for analysis of nonlinear systems[END_REF].

Much work has been done towards the verification of Simulink models [START_REF] Alur | Symbolic analysis for improving simulation coverage of Simulink/Stateflow models[END_REF], [START_REF] Zhan | Formal Verification of Simulink/Stateflow Diagrams: A Deductive Approach[END_REF]. A promising group of approaches can be categorized as verification by simulation [START_REF] Donzé | Breach, a toolbox for verification and parameter synthesis of hybrid systems[END_REF]. Donzé presented a MATLAB/Simulink based tool, Breach, which performs simulation-based verification (approximate reachability analysis) and conducts efficient signal monitoring of properties and requirements. Breach facilitates the computation and property investigation of large sets of trajectories, but it still cannot provide absolute confidence in the simulation results. Another MATLAB toolbox that is designed to be seamlessly integrated into the model based design process of MATLAB/Simulink is S-Taliro [START_REF] Annpureddy | Staliro: A tool for temporal logic falsification for hybrid systems[END_REF]. S-Taliro conducts fast and efficient simulations, but intrinsically relies on gridding, restricting the formal focus on falsification. The MATLAB/Simulink-based tool C2E2 [START_REF] Fan | Automatic reachability analysis for nonlinear hybrid models with C2E2[END_REF] generalizes simulation trajectories to families of trajectories by deriving a neighborhood around the simulated trajectory in which all trajectories have equivalent behavior.

All the above verification by simulation approaches have in common that the set of initial states must be sampled. Since the number of required samples can increase exponentially with the number of state variables, this can limit the approach to systems with low-dimensional initial states. Also, simulation-based verification techniques can be used to verify properties for a bounded-time horizon, but cannot be applied for unbounded time. The tool HySon [START_REF] Bouissou | HySon: Set-based simulation of hybrid systems[END_REF] performs set-based simulation directly on a Simulink model, and allows to compute a good approximation of the set of all possible executions. However, the technical details suggest that it may have its drawbacks when analyzing hybrid systems for an unbounded switching horizon. In addition, HySon is not publicly available.

The main contribution of this work is to introduce a compositional syntactic hybridization method. This method is suitable for Simulink models, facilitates the construction of verification models and takes advantage of the on-thefly composition of hybrid systems that is supported by the SpaceEx platform.

III. CONSTRUCTING VERIFICATION MODELS

In this section, we present the steps of our approach (portrayed in Fig. 2). To clearly illustrate our methodology, the proposed steps are applied on a rotational pendulum model.

A. Model-based design with Simulink

The modeling and the control design is undertaken with MATLAB/Simulink, through the interconnection of blocks, signals, and systems. Simulink [START_REF]MATLAB 9.0 and Simulink 8.7[END_REF] is a graphical programming environment for modeling, simulating and analyzing dynamical systems. It enables hierarchical modeling, keeping functionally related models together, and simplifying the overall design process by means of abstraction.

Example III.1 The Simulink model for the rotational pendulum is shown in Fig 3 . The system produces simulation traces of the pendulum angle over time, when it is released from rest.

Fig. 3: Simulink model for the rotational pendulum

Simulink relies on must semantics, also known as urgent or as-soon-as-possible (ASAP) semantics. That means that discrete events/transitions occur as soon as a given condition (guard) is satisfied. On the other hand, most formal tools for reachability analysis use may semantics, demonstrating a broader set of behaviors. For the purposes of this paper, we consider that SL2SX takes care of these semantic differences [START_REF] Minopoli | From simulation models to hybrid automata using urgency and relaxation[END_REF].

B. Estimation of the signal range

In this step, the goal is to get bounds on the behavior (min, max) of the input signals of the Simulink blocks that cannot be described by linear or hybrid dynamics. The smaller the ranges of the signals, the smaller the number of locations that is required by the PWA abstraction, given a desired error bound. There are different ways to estimate them, such as simulations, interval analysis, or Monte Carlo methods. In this paper, we use the Breach [START_REF] Donzé | Breach, a toolbox for verification and parameter synthesis of hybrid systems[END_REF] toolbox for a (not necessarily conservative) estimation of the signal range.

Note that the signal range serves only as an indication for the hybridization step that is presented in the next section. The approximation is equipped with out-of-range scopes. So, in case the range is shown to be insufficient during reachability computations, it is revised (enlarged).

Example III.2 For the rotational pendulum, we estimate the range of the signal that acts as an input in the nonlinear block. The range is then enlarged by a percentage. A set of simulations for uncertain initial conditions and a Sobol distribution (quasi-random number sequence) are shown in Fig. 4, where we plot the angle θ as a function of time.

D. Hybridization

This step aims to generate PWA approximations for the Simulink blocks that are not handled automatically by the translator, either because no exact translation is available (e.g. nonlinearities) or because translation cannot be applied (e.g. Embedded MATLAB Function). The PWA representation is a special class of hybrid automata admitting both discrete events, i.e. jumps, switches, and continuous dynamics in the form of piecewise affine functions.

After computing PWA approximations for the nonlinear Simulink blocks, we integrate them with the original XML file (constructed in the previous step). In this way, we get the complete SpaceEx model in SX format, combining the exactly translated blocks from SL2SX and the overapproximated blocks from syntactic hybridization. The resulting model can then be fed into the SpaceEx verification platform.

Example III. [START_REF] Zhan | Formal Verification of Simulink/Stateflow Diagrams: A Deductive Approach[END_REF] As for the rotational pendulum, the nonlinear function (sin) of the rotational pendulum is overapproximated by a PWA function. We consider the linearization domains to be boxes. The approximation errors are computed for each domain (location in the resulting hybrid automaton) and are added in the form of non-deterministic input in the invariants. The constructed SpaceEx base component is shown in Fig. 6.

E. Reachability analysis with SpaceEx

The reachability analysis is undertaken by one of the SpaceEx Analysis algorithms (e.g. STC, LGG or PHAVer). SpaceEx composes the individual components on the fly, instantiating only the part of the model that is relevant. SpaceEx supports safety verification problems, but it is also possible to check more complicated properties or control specifications. Different ways to transform control objectives to reachability problems are presented in [START_REF] Loukkas | An interval technique to check the performance of control laws applied to wind turbines[END_REF].

Example III.5 Computing the reachable sets of the rotational pendulum (for the STC SpaceEx scenario, a flowpipe tolerance of 0.01, a global time horizon of 1s), we get the phase portrait shown in Fig. 7. The pendulum is released from the most upward position with an uncertain but bounded initial speed.

IV. COMPOSITIONAL SYNTACTIC HYBRIDIZATION

The objective of the compositional syntactic hybridization is to approximate in a compositional manner the original model with a hybrid automaton with PWA dynamics. Three main steps are involved: syntactic decomposition, replacing the original system by an equivalent one with extra variables; hybridization, constructing a PWA approximation for each domain and providing a sound over-approximation of the original system by adding an error term; and finally HA composition, where the PWA model is transformed into a hybrid automaton in SX format. The two first steps are introduced in Section IV-A, where the mathematical details are presented. The third step is presented in Section IV-B and it is implementation dependent.

A. Syntactic PWA Aprroximation

In the following paragraphs, we explain the technical details involved, considering a nonlinear differential equation

dx dt = f (x), x ∈ R n . ( 1a
)
This ODE is assumed to be regular (f is Lipschitz of constant L > 0 over the state space X ⊂ R n ). The method can be extended to semi-explicit differential-algebraic equations (DAE's).

1) Syntactic Decomposition: The decomposition consists in constructing a new system where nonlinear terms are replaced by auxiliary variables,

dx dt = g(x, y), y ∈ R m , (2a)
y = h(x, y). ( 2b
)
Here y is a vector of auxiliary variables, g(x, y) ∈ R n is linear in both x and y, and h(x, y) ∈ R m includes all the nonlinear terms, m, of the original system, as explained in detail below. Notice that we have replaced the original system by a linear ODE in a higher-dimensional space, R n+m , coupled with a set of nonlinear algebraic constraints. Moreover, this step is exact. Furthermore, let V i ⊆ {x 1 , . . . , x n } for i ∈ {1, . . . , m} be the variables involved in the i-th nonlinearity, and let p i = |V i | denote the number of variables in such expression. It should be noted that with a sufficient number of auxiliary variables, we can assume that h i (x) satisfies 1 ≤ p i ≤ 2 for all i.

2) PWA Approximation: We consider a set of nonoverlapping domains, R ij , which cover the operational range of the variables in V i , where j is a label for each individual domain. For each R ij , we perform a PWA linearization of h i . Hence, (2a)-( 2b) is replaced by

dx dt = g(x, y), y ∈ R m , (3a)
y = ĥ(x, y), ( 3b
)
where ĥ is a vector of PWA functions.

Let op denote the operating point in the domain R ij . Using Taylor's formula with Lagrange remainder, for each

1 ≤ i ≤ m, ĥi (x, y) = h i,op + ∂h i ∂x op (x -x op ) + ∂h i ∂y op (y -y op ), (4)
and

h i (x, y) -ĥi (x, y) = 1 2 (x -x op ) T ∂ 2 h i ∂x 2 ξ (x -x op )+ (5) 1 2 (y -y op ) T ∂ 2 h i ∂y 2 ξ (y -y op ) + (x -x op ) T ∂ 2 h i ∂x∂y ξ (y -y op ),
where ξ = (ξ x , ξ y ) ∈ R n+m is an intermediate point in the interval ξ x ∈ {x op + a(x -x op ), a ∈ [0, 1]}, and similarly for ξ y . The right-hand side of Eq. ( 5) is the Lagrange remainder, whose resulting values over the domain R i are used to estimate the approximation error [START_REF] Berz | Computation and application of Taylor polynomials with interval remainder bounds[END_REF]. The linearization errors h are computed by evaluation of the Lagrange remainder, and satisfy y = h(x, y) ∈ ĥ(x, y) ⊕ h . In this paper, we assume that R i are boxes. In the case of a box, the point which minimizes the absolute value of the Lagrange remainder is its center [START_REF] Althoff | Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization[END_REF]. Several interesting alternatives exist, notably simplices [START_REF] Asarin | Hybridization methods for the analysis of nonlinear systems[END_REF].

Example IV.1 For n = 4, consider the polynomial vector field

f = [x 1 -x 2 x 3 x 4 , x 1 x 2 -x 4 , -x 3 x 4 , x 2 -x 3 ]
. Introducing the auxiliary variables y 1 = x 3 x 4 , y 2 = x 1 x 2 and y 3 = x 2 y 1 , then g(x, y) = [x 1 -y 3 , y 2 -x 4 , -y 1 , x 2 -x 3 ] is a linear ODE and h(x, y) = [x 3 x 4 , x 1 x 2 , x 2 y 1 ] is a nonlinear algebraic equation of degree two with m = 3 elements. Consider a PWA approximation, f , based on a rectangular partitioning of the state space, with elements of size in each dimension. Then, f leads to O(1/ 4 ) elements, while the PWA approximation of h only to O(m/ 2 ) elements.

Example IV.2 Now, we consider a 9-dimensional genetic model adapted from the one presented in [START_REF] Klipp | Systems biology in practice: concepts, implementation and application[END_REF], [START_REF] Chen | Decomposed reachability analysis for nonlinear systems[END_REF]. The model is described by polynomial dynamics of the form,

f (x) = [3x 3 -x 1 x 6 , x 4 -x 2 x 6 , x 1 x 6 -3x 3 , x 2 x 6 -x 4 , 3x 3 + 5x 1 -x 5 , 5x 5 +3x 3 +x 4 -x 6 (x 1 +x 2 +2x 8 +1), 5x 4 +x 2 - 0.5x 7 , 5x 7 -2x 6 x 8 + x 9 -0.2x 8 , 2x 6 x 8 -x 9 ]
. Introducing the auxiliary variables

y 1 = x 1 x 6 , y 2 = x 2 x 6 , y 3 = x 6 x 8 , then g(x, y) = [3x 3 -y 1 , x 4 -y 2 , y 1 -3x 3 , y 2 - x 4 , 3x 3 + 5x 1 -x 5 , 5x 5 + 3x 3 + x 4 -y 1 -y 2 -2y 3 - x 6 , 5x 4 + x 2 -0.5x 7 , 5x 7 -2y 3 + x 9 -0.2x 8 , 2y 3 -x 9
] is a linear ODE and h(x, y) = [x 1 x 6 , x 2 x 6 , x 6 x 8 ] is a nonlinear algebraic equation of degree two with m = 3 elements. Consider a PWA approximation, f , based on a rectangular partitioning of the state space, with elements of size in each dimension. Then, f leads to O(1/ 9 ) elements, while the PWA approximation of h only to O(m/ 2 ) elements. Instead of gridding a 9-dimensional state-space, we only have to approximate with PWA functions three second order statespaces. This example highlights the usability of the syntactic approach for the cases that repeated nonlinearities appear.

B. Compositional Hybridization

Through the syntactic PWA approximation, we have produced a Linear ODE with PWA algebraic constraints. In order to feed this model to one of the available reachability tools, we must describe it as a network of hybrid automata. Each hybrid automaton corresponds to the PWA approximation of one nonlinearity. Each piece of the PWA approximation corresponds to one location in the corresponding automaton.

In the SX file format, used by SpaceEx and other tools, a model consists of components, which are either hybrid automata or networks of hybrid automata. A component can be instantiated inside a network, possibly remapping variables to other variables or replacing them with constant values. Note that an ODE or an algebraic constraint can be trivially embedded in a hybrid automaton with a single location. The ODE becomes the flow-constraint of the location and the algebraic constraint its invariant.

Expressing the PWA approximation to this setting, the linear ODE is modeled by a single (trivial) automaton. Each PWA constraint y i = ĥi (x, y) corresponds to a hybrid automaton with one location per piece. The locations of adjacent pieces are connected through transitions. The approximation error is expressed by extra variables with range h , and the error threshold µ > 0 is an upper bound on the maximum value it can take (in some chosen norm || • ||).

Example IV.3 Getting back to Example IV.2, we model the genetic system with Simulink (see Fig. 8). The nonlinearities correspond to products and are highlighted with different colors. Blocks with the same color represent the same nonlinear operator.

After translating the Simulink model and estimating the signal ranges, we are ready to perform syntactic hybridization on the components. As already mentioned in Section IV-A, we don't need to partition a 9-dimensional space, but only to construct 3 two-dimensional PWA approximations for the products. Note that the approximations correspond to base components in SpaceEx and can be seen as templates (building blocks). In this respect, it is possible to be reused. As a result, only one SpaceEx base component would be required to describe the three red blocks, one for the blue blocks, and another one for the yellow blocks.

Standard algorithms for reachability analysis, such as those inside SpaceEx, take as input a single hybrid automaton with ODE dynamics. To go from the multi-component input model to this form, the reachability tool performs two operations: First, it combines the components through a process called parallel composition. Second, it eliminates the algebraic constraints to obtain an ODE. In principle, parallel composition means building the product automaton, whose locations consist of the cross-product of the locations of the components. A model with m components of k locations each thus has a product automaton with m k locations. However, tools such as SpaceEx construct the product automaton on the fly, instantiating only the reachable locations. Similarly, the conversion from linear DAE to ODE is only carried out on the instantiated locations. The conversion can be carried out efficiently by Gauss-Jordan elimination [START_REF] Donzé | Modular, hierarchical models of control systems in SpaceEx[END_REF]. The underlying theory is explained in [START_REF] Gerdin | Identification and estimation for models described by Differential-Algebraic Equations[END_REF].

Example IV.4 (continuation of example IV.3.) Conducting reachability analysis with SpaceEx for a global time horizon of 1s, with a flowpipe tolerance of 0.01 and the STC scenario (octagonal directions), we compute the reachable sets, as illustrated in Fig. 9. In principle, the outlined procedure enables us to approximate the reachable set of the original dynamics with arbitrary precision. Let Φ(t, x) denote the trajectory starting from x evaluated at time t. The reachable set of the system from a set of initial points X 0 ⊆ X during the interval [0, t] is defined as

Reach(t, X 0 ) = {y = Φ(τ, x) : τ ∈ [0, t], x ∈ X 0 }. (6)
The approximate system converges to the original system, as expressed by the following theorem.

Theorem 1 (see [START_REF] Asarin | Reachability analysis of nonlinear systems using conservative approximation[END_REF]) The Hausdorff distance between the reachable set of (2a)-(2b) and the reachable set computed through hybridization, (3a)-(3b), from time 0 to a final time T > 0 satisfies

d H Reach f (T, X 0 ), Reach f (T, X 0 ) ≤ 2µ L e LT -1 , ( 7
)
where µ is the error threshold, L is the Lipschitz constant of the original, nonlinear function and f is the PWA approximation.

C. Case Study: Wind Turbine

The wind turbine benchmark from the ARCH workshop poses a challenging and relevant industrial model [START_REF] Schuler | Hybrid modelling of a wind turbine (benchmark proposal)[END_REF]. It is designed with MATLAB/Simulink and it is a largescale model with many nonlinearities. We used syntactic hybridization to transform the model into an approximative PWA one in SX format and we conducted the approximation component-wise. Ten nonlinear blocks were approximated and Table I presents the corresponding results. The resulting model only has 72 locations in all components combined. The standard hybridization would lead to an O(1/ n ) number of locations. Taking into account that the dimension of the state-space is 7 and considering 5 locations per state variable, the standard method would yield 5 7 = 78125 locations. That indicates a considerable difference between the hybridization methods in terms of the model size.

V. CONCLUSIONS

Model transformation plays an important role in bridging the gap between industrially relevant models and verification tools [START_REF] Alur | Symbolic analysis for improving simulation coverage of Simulink/Stateflow models[END_REF]. This work aims to assist the application of hybrid system reachability tools to models designed with MATLAB/Simulink. We propose a methodology to construct verification models out of Simulink systems. We make use of the SL2SX translator to handle the mechanical aspects of the translation to hybrid automata and the Breach toolbox to get bounds on the signal ranges. For the blocks that are not exactly translated, e.g. nonlinearities, we propose a method, which we call compositional syntactic hybridization.

Unlike standard state-space hybridization methods, we do not operate over the fully composed (flattened) model, but perform the PWA approximations component-wise. In this way, we can obtain a significant reduction in terms of model size. The constructed verification model consists of a network of hybrid automata and is described in the SX format. It can then be fed into the SpaceEx platform or other verification tools through the HYST translator. Using SpaceEx for the reachability computations, we can take advantage of the onthe-fly composition and instantiate only the reachable parts of the approximation. Note that our compositional hybridization can be applied not only to the dynamics, but also to algebraic and initial constraints. On an industrial benchmark, the wind turbine, our approach leads to a very compact model that is orders of magnitude smaller than a standard hybridization model. The next step is to improve the reachability tools so that they can take advantage of these compositional models. The primary objective is to instantiate as few locations and transitions as possible during the analysis. There are three issues to address. The first is the on-the-fly composition and instantiation of the models, which can reduce the number of instantiated locations of the product automaton. The second is a compositional pre-processing of the components, where we utilize the pre-image of the target invariant when checking which transitions are enabled. The third direction would be to perform compositional mapping of the initial states. In the case of SpaceEx platform, the identification of the initial conditions is done through enumeration. However, enumeration of the locations of the product automaton is an operation that does not scale. Applying compositional reasoning would allow us to identify the initial locations and instantiate as few locations as possible.

Fig. 1 :

 1 Fig.1: Constructing verification models from a Simulink system, with m = 3 nonlinear blocks. The syntactic approach can lead to approximations of smaller size and less instantiated locations.

Fig. 2 :

 2 Fig. 2: Methodology for constructing verification models; input: Simulink, output: PWA decomposed model (left). The syntactic hybridization steps are presented (right). Each step (in Latin) is analyzed in the following sections.

Fig. 4 :

 4 Fig. 4: Range estimation of the input signal of the nonlinear block (sin) of the rotational pendulum θ (rad) vs. time (s); conducted with Breach simulations.

Fig. 5 :

 5 Fig. 5: SpaceEx model (SX) of the rotational pendulum constructed by the translator. This is the network component, and we highlight, in red, the block which corresponds to a trigonometric function.

Fig. 6 :

 6 Fig. 6: PWA of the nonlinear component (sin) of the rotational pendulum. Only 6 locations are shown (out of 40). Here the non-deterministic input w 1 represents the approximation error.

Fig. 7 :

 7 Fig. 7: Reachability results of the rotational pendulum for a global time horizon of 1s. Phase portrait of the angular position (rad) and angular speed (rad/s).

Fig. 8 :

 8 Fig. 8: Genetic model in Simulink. Colored blocks correspond to nonlinearities. Red: x 1 x 6 , cyan: x 2 x 6 and yellow: x 6 x 8 .

Fig. 9 :

 9 Fig. 9: Reachable sets of the genetic model with SpaceEx. The variables x 1 and x 2 are displayed.

TABLE I :

 I Wind Turbine Benchmark -breakdown into blocks. Blocks, 1-9, are approximated syntactically; the 10 th block is exactly translated; 11 th is replaced by a nondeterministic input and the rest, 12-17, are not necessary due to semantic differences.

ACKNOWLEDGEMENT

The authors would like to thank Willem Hagemann and Stefano Minopoli for valuable discussions.


Elodie Kobeissi
email: elodie.kobeissi@ens-cachan.fr

Pierre-Yves Piriou
email: pierre-yves.piriou@edf.fr

Jean-Marc Faure
email: jean-marc.faure@ens-cachan.fr

Formal Verification of Safety Analysis Models of Repairable and Reconfigurable Systems

Keywords: Model Based Safety Analysis/Assessment, Formal verification, Dynamic properties, Generalized BDMP, Stability, Computation Tree Logic. 

This paper proposes a method to formally check whether formal properties hold on a dynamic model which has been designed by experts for Model Based Safety Analysis/Assessment. As repairable and reconfigurable systems are considered, this model is assumed to be described in the Generalized Boolean-logic Driven Markov Processes (GBDMP) modelling framework. Translation rules are given to obtain a formal model that describes correctly the evolution of the initial model with the semantics of the verification tool. The approach is exemplified on a simple case of standby redundancy.

INTRODUCTION

Model Based Safety Analysis/Assessment (MBSA) techniques are aiming to compute dependability attributes (instantaneous and asymptotic reliability and availability, minimal cut sequences, for instance) by analysing models which have been previously designed by safety/dependability experts [START_REF] Bouissou | A new formalism that combines advantages of fault trees and Markov models: Boolean logic Driven Markov Processes[END_REF][START_REF] Lipaczewski | Comparison of modeling formalisms for safety analyses: SAML and AltaRica[END_REF]. Accurate and trustworthy computation results require that the initial models are flawless; unfortunately, as any model based on expertise, these models may contain flaws coming from ambiguous or misunderstood specifications, human errors, etc. The overall aim of this work is to contribute to tackle out (or at least limit) this issue by proposing an approach of formal verification of models that will be used later for MBSA.

The modelling framework which was selected for this work is GBDMP (Generalized Boolean-Logic Driven Markov Processes). This framework has been developed for MBSA of dynamic, repairable, and reconfigurable systems [START_REF] Piriou | Generalized Boolean logic Driven Markov Processes: A powerful modelling framework for Model-Based Safety Analysis of dynamic repairable and reconfigurable systems[END_REF]; it particularly enables fine modelling of different reconfiguration strategies and of the failure of these strategies. Construction rules of well-formed GBDMP models have also been defined in the above reference to ensure syntactic consistency. However, these rules do not guarantee that a well-formed model satisfies dynamic properties, i.e. properties that consider evolutions of the model during time. Verification of such properties require to translate the model into the language of a verification tool then to check whether the properties hold or not on the model (Figure 1). Last, only logical time is considered in this work because functional correctness must be checked before time correctness; hence, the properties to check will be stated in Computation Tree Logic (CTL). Formal verification will be performed by using the widespread NuSMV model-checker; therefore, the formal model of the GBDMP will be a transition system in the language of this tool.

The outline of the paper is the following. A state of the art on formal verification of models for MBSA is given in the next section. The GBDMP modelling framework is shortly presented in section 3 while section 4 focuses on translation of a GBDMP model into the NuSMV language. Examples of positive and negative verification of dynamic properties are detailed in section 5. Finally, concluding remarks and perspectives for future works are drawn up in the last section. The objective of [START_REF] Thums | Model-checking FTA[END_REF] is similar, but the semantics of fault trees is given in CTL and the components behaviour is modelled with timed automata. Moreover, two kinds of gates are defined: decomposition gates with the Boolean semantics and cause consequence gates, where time may elapse between the causes and the consequence. A last contribution in this domain that deserves to be mentioned is the work reported in [START_REF] Ortmeier | Formal fault tree analysis -Practical experiences[END_REF] where fault trees are described in Interval Temporal Logic (ITL) and state charts model the behaviour of the system components.

To allow modelling of failures that depend on the order of fault events, a modelling framework named State/Event Fault Trees (SEFTs) is proposed in [START_REF] Kaiser | State/event fault treesa safety analysis model for software controlled systems[END_REF]. This framework subsumes fault tree, deterministic state machines and Markov chains. Analysis of models developed in this frame is performed by translating the component models into Deterministic and Stochastic Petri Nets (DSPNs). This framework is implemented in the ESSaRel (Embedded Systems Safety and Reliability Analyser) tool. Last, [START_REF] Bozzano | Safety assessment of AltaRica models via symbolic model-checking[END_REF] proposes recently to check properties of AltaRica OCAS models by using the NuSMV verification tool while [START_REF] Sharvia | Integrating model checking with HiP-HOPS in model-based safety analysis[END_REF] focuses on verification of HiP-HOPS models with the same tool. Despite their interest, none of these works has addressed models for dynamic, repairable, and reconfigurable systems. This study is aiming at filling this gap.

BRIEF REMINDER ON GBDMP

Overall description

A GBDMP model is a 3-tuple composed of:

• an extended fault tree that describes the physical and functional structure of the considered critical process; this fault tree comprises classical gates and leaves as well as switches to introduce explicitly reconfigurations of the process.

• a set of Switched Markov Processes (SMP) which represent the functional and dysfunctional behaviours of the components of the critical process; a SMP is associated to every leaf l  L, where L is the set of leaves of the fault tree.

• a set of Moore machines which model reconfiguration strategies; a Moore machine is associated to every switch s  S, where S is the set of switches of the fault tree.

The three components of this tuple are interconnected. Hence, syntactic rules to build well-formed GBDMP models have been stated; the input/output alphabets of a Moore machine must be consistent with the inputs/outputs of the associated switch, for instance. It will be assumed in what follows that the model to formally verify is well-formed. For space reasons, it is not possible to give a complete syntax of GBDMP; detailed presentations can be found in [START_REF] Piriou | Generalized Boolean logic Driven Markov Processes: A powerful modelling framework for Model-Based Safety Analysis of dynamic repairable and reconfigurable systems[END_REF] and [START_REF] Piriou | From safety analysis of reconfigurable systems to design of faulttolerant control strategies[END_REF].

Example

The example of Fig. 2 is composed of two groups of redundant components {C1a, C1b, C1c} and {C2a, C2b} in a series. The nature of these components does not matter; the only assumptions are that each component can be activated and deactivated and may fail and be repaired when active and inactive. The configuration of each group is managed by a controller (D1 for the first group and D2 for the second one) that performs a reconfiguration strategy:

• for the first group, C1a and C1b are by default active and C1c inactive. Then, if one of the active components fails while the other two ones are not faulty, the failed component is deactivated, and the inactive component activated. • for the second group, if C2a is not faulty, C2a and C2b must be respectively active and inactive, whatever the dysfunctional state (faultless/faulty) of C2b; if C2a is faulty, C2a and C2b must be respectively inactive and active.

Fig. 2. Standby redundancy example

The GBDMP model of this example, with the above reconfiguration strategies, is shown at Fig. 3. In the fault tree, the leaves C1a, C1b and C1c are linked by dashed arrows to the switch S1; this means that the SMP associated to these leaves send/receive variables to/from the Moore machine. This latter is associated to this switch to model the first strategy where the activation status (active/inactive) of the components C1a, C1b and C1c depends on the failure statuses of these components. A similar explanation can be given for the leaves C2a and C2b and the switch S2. It must be noted nevertheless that no variable is sent to S2 from the leaf C2b because the activation status of C2a does not depend on the failure status of C2b, for the second strategy.

Moreover, the success of a reconfiguration strategy is not guaranteed because the component where it is implemented, a Programmable Logic Controller for instance, may fail. Two failure modes of this controller will be considered hereafter:

• Frozen: no more switching is possible, i.e. the configuration of the process remains the same;

• Bad contact: there is no connection between D1 (D2) and the group {C1a, C1b, C1c} ({C2a, C2b}).

Evolutions of a GBDMP model

A GBDMP model describes a stochastic system whose components may fail and be repaired. However, when focus is put on formal verification of functional (and not quantitative) dynamic properties of this model, the random transitions of its SMP must be replaced by deterministic transitions labelled by the appropriate event. Hence, the transitions that model failures (repairs) will be labelled with failure (repair) events, and not failure (repair) rates, in what follows. With this modelling, the global state of a GBDMP model is completely defined by two sets of variables:

• The set of state variables of the Moore machines.

The state variable of the Moore machine associated to the switch s ∈ S will be noted Us.

• The set of state variables of the Switched Markov Processes. The state variable of the SMP associated to the leaf l ∈ L will be noted Xl.

Moreover, for each node (leaf or gate) n of the fault tree, three other variables can be computed from these two sets of state variables, as detailed in [START_REF] Piriou | Generalized Boolean logic Driven Markov Processes: A powerful modelling framework for Model-Based Safety Analysis of dynamic repairable and reconfigurable systems[END_REF]):

• Fn: Boolean variable that represents the failure status of the node n (Fn=True/False means that n is faulty/faultless).

• Rn: Boolean variable that represents the requirement status of the node n (Rn=True/False means that n is required/not required to perform the process function).

• Mn: Integer that represents the activity status of the node n (Mn=k means that n is in the k th activity mode).

All these variables (state variables and variables associated to a node) are interdependent and must be computed from the dependency graph that can be built from the extended fault tree. Therefore, the state space of a GBDMP model includes two kinds of states: stable states and unstable states. A state is stable when the state variable of every SMP is consistent with the failure and activity statuses of the corresponding leaf l (if Fl is True and Ml equal to k for instance, the SMP associated to l must be in a failure state of the k th activity mode), unstable otherwise.

A GBDMP model evolves from state to state on occurrence of two types of events: spontaneous events and provoked events. A spontaneous event occurs only when the current state is stable; it corresponds to a failure or repair event or an operator request, like a phase change. The spontaneous events are associated to the transitions in solid lines in the SMP. A provoked event occurs when the state is unstable; it is the consequence of a spontaneous event, like the transition of the activity mode of a spare component from inactive to active. The provoked events are associated to the transitions in dashed lines in the SMP. The spontaneous events are asynchronous (their occurrence dates are different) whereas the provoked events which are consequences of a given spontaneous event occur synchronously (at the date of the spontaneous event which caused them). The global state of the GBDMP model is stable when the Boolean variable Stable which is the conjunction of every Stable_Xl and Stable_Um variable is True.

The main issue in this translation is modelling of the spontaneous events which are asynchronous in the GBDMP framework while NUSMV considers by default that two events may occur simultaneously. To avoid incorrect evolutions, an event selector must be introduced (Figure 5). The role of this automaton is to prevent simultaneous failure/repair events, from the SMP associated to the leaves of the fault tree, to occur. From the initial state, Choice, when the global state of the GBDMP model is stable, only one state of the event selector can be activated and then only one spontaneous event can occur. When a non-initial state of this automaton is active, a variable l_evol where l is the name of the corresponding leaf, is then set; this variable means that the SMP associated to this leaf can evolve and must be introduced as a guard in every transition of the SMP. The initial state of the event selector becomes active again once the SMP has evolved. To sum up, introducing this automaton and these guards in the NUSMV code permits to describe correctly the evolutions of the GBDMP model with the semantics of the verification tool.

Notations

Seven properties that must be satisfied by the GBDMP model will be given hereafter. They concern the behaviour of the components, redundant components of the critical process and controllers where the reconfiguration strategies are implemented, which are modelled by SMP, the reconfiguration strategies themselves, modelled with Moore machines, and the top event of the extended fault tree. These properties will be exemplified on the example of figure 3 by using CTL expressions. Hence, the notations of the state and path quantifiers in this logic must be reminded:

• X represents the next state of the current state,

• F means "for at least one state (there exists in the Future) of a path", • G means "for every state (Globally) of a path", • A means "along All paths" from the current state, • E means "along at least one (there Exists) path" from the current state. With these notations, EF , where is a logic expression, means "there exists a path where there exists at least one state where  is verified" and AG "for every path and for every state of the path,  is verified", for instance.

Properties presentation

Properties on the behaviour of the components

The first property concerns the behaviour of a redundant component (C1a, C1b, C1c, C2a, C2b in the example); the three following properties the behaviour of a controller where a reconfiguration strategy is implemented.

Property 1: A redundant component can be activated (1a) and deactivated (1b); it can also fail (1c) and be repaired (1d). This property can be formalized by the four CTL formulas below, where EF ( & AX!) means that there exists in the future a state where  is true and for all immediate next possible states,  is false:

EF (Ml=0 & AX Ml>0) (1a) EF (Ml >0 & AX Ml =0) (1b) EF (!Fl & AX Fl) (1c) EF (Fl & AX!Fl) (1d)
where Ml represents the activity status of the leaf l associated to the component and Fl the failure status of this component; it is assumed that Ml=0 means that the component is inactive.

Property 2: A controller may fail in two ways (two failure modes) (2a) and be repaired when failed (2b).

For the component associated to the leaf D1 for instance, this property is formalized by the following two expressions:

EF (!FD1 & AX ((FD1 & (XD1=2 | (XD1=1))); (2a) EF ((FD1 & (XD1=2 | (XD1=1)) & AX (!FD1)); (2b)
where XD1=1(2) means that the failure mode is bad contact (frozen).

Property 3: When the failure mode is bad contact, every component of the group of redundant components controlled by the controller is inactive.

For the example of figure 2, this property is formalized by two CTL expressions (one for each group of components).

AG ((XD1=1 & Stable

) → (MC1a=0 & MC1b=0 & MC1c=0)); (3-1) AG ((XD2=1 & Stable) → (MC2a=0 & MC2b=0));
(3-2) It must be noted that the variable Stable is introduced in these statements because the property holds only in the stable states of the model. This remark can be made for the other following properties.

Property 4: When the failure mode is frozen, every component of the group of redundant components which is controlled by the controller remains in its current state.

AG (XD1=2 → (Stable_XC1a & Stable_ XC1b & Stable_ XC1c)); (4-1) AG (XD2=2 → (Stable_ XC2a & Stable_ XC2b));
(4-2)

Properties on the reconfiguration strategies

To check whether the two reconfiguration strategies which have been defined at 3.2 have been correctly modelled in the GBDMP, the properties 5 and 6 are to be verified. Property 5 will be expressed in an event-oriented form whereas property 6 will be given in a state-oriented form to show that informal specifications can be formalised in these two forms. In the expression of property 5, the construction "A [! Stable U p]" is used to specify the assertion "p will be True in the next stable state".

Property 5: If the controller D1 is faultless, by default (no component is faulty), C1a and C1b are arbitrarily activated (5a). Then, if one of the active components fails while the other two are faultless, disable the failed component and activate the inactive component ((5b-1) to (5b-3)).

AG ((!FD1 & Stable

& !FC1a & !FC1b & !FC1c) → (MC1a=1 & MC1b=1 & MC1c=0)); (5a) AG ((!FC1a & AX FC1a) & (MC1a=1 & !FC1b & !FC1c & !FD1)→ AX (A [!Stable U (MC1a=0 & MC1b=1 & MC1c=1)])); (5b-1) AG ((!FC1b & AX FC1b) & (MC1b=1 & !FC1a & !FC1c & !FD1)→ AX (A [!Stable U (MC1a=1 & MC1b=0 & MC1c=1)])); (5b-2) AG ((!FC1c & AX FC1c) & (MC1c=1 & !FC1a & !FC1c & !FD1)→ AX (A [!Stable U (MC1a=1 & MC1b=1 & MC1c=0)])); (5b-3)
Property 6: If the controller D2 is faultless and the model stable, whatever the dysfunctional state of C2b, if C2a is not faulty, then C2a and C2b are respectively active and inactive (6a); if C2a is faulty, C2a and C2b are respectively inactive and active (6b).

AG ((

!FD2 & Stable & !FC2a) → (MC2a=1 & MC2b=0)); (6a) AG ((!FD2 & Stable & FC2a) → (MC2a=0 & MC2b=1)); (6b)

Property on the top event

The top event of the extended fault tree is False when the critical system of figure 2 performs correctly its function, True otherwise. The following property can then be stated.

((!FD1 & !FD2 & Stable ) → (count (!FC1a, !FC1b, !FC1c) ≥ 2 & count (!FC2a, !FC2b) ≥ 1) ↔ !TE ); ( 7
)
where count is a function which counts the number of variables that are True in a set of Boolean variables.

Detection of errors in the GBDMP model

As the model of figure 3 is correct, every property which has been defined at the previous section holds on this model. To check whether errors during modelling can be detected by the verification tool, two modifications have been brought separately:

• The top gate (OR in the correct model) has been replaced by a AND gate, which is obviously a basic error. • The label of the transition from q3 to q1 in the Moore machine M1 has been changed. The first error is detected because Property 7 does not hold anymore; a counterexample that shows that it is possible in this case to have a False TE while two components in the group {C1a, C1b, C1c} are faulty (Figure 6) is provided by the model checker. Such a counterexample surely eases analysis the GBDMP model and detection of the modelling error. When only the second error is introduced, both properties 5a and 7 are not satisfied. A counterexample shows that the three components of the group remain inactive after the controller D1 has been repaired. Analysis of the Moore machine that describes the strategy implemented in this controller permits to detect the erroneous label. In figure 7, the states in solid lines correspond to stable states of the system, while the dotted ones correspond to unstable states. Fig. 7. Counter example for the second modelling error.

CONCLUSIONS

This paper has shown that modelling errors in a model developed for MBSA of repairable and reconfigurable systems can be detected by using formal verification techniques. Solutions to translate this model in the language of the selected verification tool while keeping the initial semantics have been proposed and some examples of errors that can be detected have been given. The experiments show that the approach scales well. Future work is aiming at facilitating verification by non-experts in temporal logic; development of libraries of CTL formulas which are built automatically from the knowledge of the structure of the critical process and the reconfiguration strategies is under investigation.

Fig. 1 .

 1 Fig. 1. Principle of the work

Fig. 3 .

 3 Fig. 3. GBDMP model of the system of Figure 2, a) Extended Fault Tree; b) Moore Machine associated to switch S1; c) Moore Machine associated to switch S2; d) SMP associated to leaf De; e) SMP associated to leaf SF This possible issue is modelled by the link between the leaf D1 (D2) and the switch S1 (S2). The activation status of a component of the groups {C1a, C1b, C1c} and {C2a, C2b} depends on the failure status of the reconfiguration controller.

Fig. 4 .

 4 Fig. 4. Part of the NuSMV code for the model of figure 3 4. TRANSLATING GBDMP MODELS IN NUSMV LANGUAGE NuSMV [Cimatti 2002] is a symbolic verification tool that checks properties on a model in the form of a transition system in a specific syntax. The aim of this section is to give the main translation rules to obtain this model from a wellformed GBDMP model. A part of the result obtained for the GBDMP of figure 3 is given at figure 4. Only the variables that define the global state (state variables of the Moore machines and SMP) are first defined in the VAR section. Translation of the evolutions of these variables is straightforward by using the operator 'next' that defines the next state. The status variables of the nodes are then computed and the stability condition is obtained, for every leaf of the fault tree, from the values of the state variables and the status variables of this leaf, as follows: Stable_Xl := ((Ml =0 & Xl in {set of states in mode 0}) | … | (Ml =k & Xl in {set of states in mode k})) & ((Fl & Xl in {set of failure states}) | (!Fl & Xl in {set of non-failure states}))

Fig. 5 .

 5 Fig. 5. Event selector for the model of figure 3 5. FORMAL VERIFICATION OF DYNAMIC PROPERTIES

Fig. 6 .

 6 Fig.6. Counter example for the first modelling error.

Property 7 :

 7 If every controller is faultless, the global function is performed if at least two components among {C1a, C1b, C1c} are faultless and at least one component among {C2a, C2b} is also faultless. This property can be formalized as follows by introducing a Boolean variable TE (Top Event); this variable is merely the complement of the failure status of the top gate of the fault tree.

	AG


C Le Bas

S Sahuguede

A Julien-Vergonjanne

A Behlouli

P Combeau

L Aveneau

Human Body Impact on Mobile Visible Light Communication Link

Keywords: Visible Light Communication, Mobile indoor optics, Channel characterizations, Body model I

Shadowing impact for indoor mobile visible light communications (VLC), is linked to the person wearing the optical receiver, depending on the body features and varying according the data rate. Here, we investigate VLC performance considering two body models with different reflectivity values. In addition, we consider a mobile receiver with random orientations to account body movements. Results show that for low data rates, a 2D body model with perpendicular receiver orientation allows suitable performance estimation. For high data rates, 2D body model can be used but link efficiency is deteriorated depending on receiver orientation variations whatever body reflectivity values are.

INTRODUCTION

As energy saving is a topic of interest, much research on Visible Light Communications (VLC) based on Light Emitting Diodes (LEDs) has been developed in recent years [START_REF] Chvojka | Channel Characteristics of Visible Light Communications Within Dynamic Indoor Environment[END_REF][START_REF] Shlomi | Visible Light Communication[END_REF][START_REF] Elgala | Indoor Optical Wireless Communication: Potential and State-of-the-Art[END_REF]. LEDs have several attractive characteristics, so they are used in many contexts in order to replace classic lighting methods like incandescent light. Actually, they are energy efficient, and they can be easily combined to perform dimming for brightness, warm management or cool lighting. In addition, LEDs structure allows fast switching abilities, so they emerge to be the best candidates for VLC networking, which main goal is to provide both illumination and communication.

VLC and more generally Optical Wireless Communications (OWC) have a number of advantages over existing wireless communication ones and constitute a solution to overcome spectrum congestion of radio frequency (RF) communications.

First, OWC spectral domain offers a high and not regulated bandwidth. In addition, optical components are cheap so it is possible to deploy a low cost network. Furthermore, this technology does not generate electromagnetic interferences, so optical network will not interfere with RF existing ones, making its application particularly promising in sensitive environment like hospitals or airplanes [START_REF] Murai | A novel visible light communication system for enhanced control of autonomous delivery robots in a hospital[END_REF][START_REF] Chan | Wireless optical links for airplane applications[END_REF].

Several outdoor applications of OWC have been considered such as intra and extra vehicular communications or to be a complement to RF in the 5G technology [START_REF] Pergoloni | Coverage optimization of 5G atto-cells for visible light communications access[END_REF]. Nevertheless, it is for indoor applications that much research was conducted and a first standard, IEEE 802.15.7, has been recently proposed, allowing regulating the use of VLC [START_REF] Boucouvalas | Standards for indoor Optical Wireless Communications[END_REF].

For indoor environments, the high power provided by the room lighting allows obtaining high Signal to Noise Ratios (SNR), increasing the interest of this kind of scenario. In addition, due to the light confinement, transmitted data are restricted to room area, contributing to information security enhancement.

Most current indoor VLC studies consider a fixed receiver with emitters placed on the ceiling. These configurations offer a high quality of service (QoS) but it has been shown that shadowing effect due to the presence of objects and/or moving person can degrade VLC performance [START_REF] Chvojka | Channel Characteristics of Visible Light Communications Within Dynamic Indoor Environment[END_REF][START_REF] Komine | A study of shadowing on indoor visible light wireless communication utilizing plural white LED lightings[END_REF][START_REF] Kizilirmak | Impact of repeaters on the performance of indoor visible light communications[END_REF].

In addition, when the receiver is mobile, for example worn by a person, the body movements induce receiver orientation variations which can also contribute to performance degradation because of link blocking cases. Methods like angle diversity can overcome this problem [START_REF] Hussein | 10 Gbps Mobile Visible Light Communication System Employing Angle Diversity, Imaging Receivers, and Relay Nodes[END_REF][START_REF] Nuwanpriya | Indoor MIMO Visible Light Communications: Novel Angle Diversity Receivers for Mobile Users[END_REF] but it is necessary to have a more complex receiver. This is difficult in the case of a system worn by a person.

In this paper, we study a classical VLC configuration constituted of several emitters fixed on the room celling but assuming a mobile and simple receiver worn by a moving person. Our contribution is to determine body impact on the performance including orientation and position variations of the receiver. The analysis is provided thanks to ray-tracing channel model with different body profiles and considering body/room reflectivity features. Performance is obtained as a function of the data rate.

The paper is organized as follows: in part II we describe the studied environment and the different body models used in our simulation. Part III presents the results obtained taking into account body impact. Finally, conclusions are given in part IV.

II. DESCRIPTION OF THE STUDIED SYSTEM

A. System Overview

The configuration studied in this paper has been considered in many VLC studies [START_REF] Ghassemlooy | Opitcal Wireless Communications System and Channel Modelling with MATLAB[END_REF][START_REF] Grubor | Bandwidth-efficient indoor optical wireless communications with white light-emitting diodes[END_REF][START_REF] Bas | Impact of receiver orientation and position on Visible Light Communication link performance[END_REF][START_REF] Tan | Simulation of MIMO channel characteristics for indoor visible light communication with LEDs[END_REF][START_REF] Barry | Simulation of multipath impulse response for indoor wireless optical channels[END_REF]. As seen in Fig. 1, it consists in an empty room (5mx5mx3m), where four VLC transmitters, each constituted of 3600 LEDs, are placed on the celling (coordinates are indicated in TABLE I). The total emitted power is 144W and ensures a minimum of 400 lux in the room [START_REF] Tan | Simulation of MIMO channel characteristics for indoor visible light communication with LEDs[END_REF].The receiver can be at any position on a 2D plan at 1.2m from the floor. The walls and the ceiling are considered as plaster surfaces and the floor as plastic one. Moreover, we assume a random uniform distribution of receiver positions within the 2D plan in order to model the receiver mobility. As mobility also implies the variations of the receiver orientations (for example due to the body movements), we define azimuthal and longitudinal angles respectively φ and θ, as illustrated in Fig. 1, whose values are uniformly and randomly chosen within φ ∈ [0, 2π] and θ ∈ [0,

B. Channel model

Usually, studied VLC configurations are static considering a receiver oriented towards the ceiling (θ=0 in Fig. 1). In this case, the VLC channel is mainly characterised by Line Of Sight (LOS) optical paths.

However, in indoor environment, it is also important to take into account optical reflections from surfaces which can impact the received power. As we consider receiver mobility and variations of receiver orientation, this is all the more significant. The global channel impulse response can be expressed as:

ℎ 𝑡𝑜𝑡𝑎𝑙 (𝑡) = ∑ ℎ 𝑖 (𝑡) = 𝑁 𝑖=1 ∑ ∑ ℎ 𝑖 (𝑘) (𝑡) ∞ 𝑘=0 𝑁 𝑖=1 (1)
Where N is the number of sources and ℎ 𝑖 (𝑘) (𝑡) is the k bounces contribution to the channel impulse response according to the i th source. It has been shown that a suitable model for indoor dynamic VLC channel can be obtained considering three reflections (k=3) of the optical signal [START_REF] Bas | Impact of receiver orientation and position on Visible Light Communication link performance[END_REF].

An accurate way to calculate ℎ 𝑡𝑜𝑡𝑎𝑙 (𝑡) uses a ray tracing algorithm relying on the Monte Carlo method. For this aim, we use a ray-based simulator developed at the XLIM Laboratory, which reduces computation time by implementing some techniques of variance reduction such as importance sampling or next event estimation methods [START_REF] Behlouli | Efficient simulation of optical wireless channel, Application to WBANs with MISO link[END_REF].

For our simulation set-up, each of the N impulse responses is determined over duration of 60ns with a step of 0.15ns and 2 million launched rays.

Various environmental elements such as walls and objects contribute to the optical reflections. The human body, which can wear the receiver device, can also have a great influence on the propagation loss.

In this paper, to study the impact of the human body on the VLC performance, we have to model its geometry on one side and its reflectance properties on the other side.

. A first question concerns the level of details to be used by simulation and its impact on impulse response. To answer this question, we consider two geometrical models, one in 2D corresponding to the profile of a generic human body (cf. Fig. 2 (a)), and another one being a highly detailed 3D mesh [18] (cf. Fig. 2 (b)). These two models represent a body of 1.72m height and 45cm width. The 3D one presents a thickness of 25cm and is composed of 1796 facets (triangles). corresponding to the reflectivity at normal incidence. The value of ρ varies from 0 to 1. A weak value represents very absorbing surfaces, while a high value is associated to very shiny materials. The surface of a human body exposed to light propagation is composed of two main parts, the skin and the clothes, and their ρ values vary according to the light wavelength. Fig. 3 resumes the corresponding range of some ρ values according to the light wavelength [START_REF] Glassner | Principles of Digital Image Synthesis[END_REF]. In conclusion, to take into account the whole range of ρ values connected to the different human skins and clothes, and so to obtain the most representative conclusions, all light simulations presented in this paper are obtained concerning two extreme mean values of , respectively = 0.1 and = 0.7. Finally, the reflectance properties of the rest of the environment (walls, ceiling and floor) are presented in the TABLE II. They have been fixed according to [START_REF] Lee | Indoor Channel Characteristics for Visible Light Communications[END_REF].

A. Impact on the optical gain

For the simulation we have considered 961 receiver positions uniformly chosen within the 2D plane. In addition the receiver is worn by a person whose body is modelled by the 2D and 3D previously described shapes and supposed to be moving with the receiver. So, we have assumed that the body is oriented towards the receiver, and located at 20cm of it. At last, for the orientation variations, 18 values for each receiver position have been simulated. They have been chosen uniformly within the interval defined in Part.II A.

The probability density function (PDF) of the optical gain H(0) considering five different scenarios is reported in Fig. 4 and5. The optical gain H(0) is defined as:

𝐻(0) = ∫ ℎ 𝑡𝑜𝑡𝑎𝑙 (𝑡)𝑑𝑡 +∞ -∞ (2)
The studied scenarios consider receiver mobility. This means that we simulate 961 receiver positions and that receiver orientation varies as defined in section II:

- First, it may be observed by comparing scenario (b) and (c) in Fig. 4 and scenario (d) and (e) in Fig. 5, that for a given value of reflectivity coefficient, the results obtained with 2D and 3D body models are quite similar. This means that it is not necessary to use a complex model which is a significant achievement given the diversity of human body. Then, if we compare scenario (a) and (b) in Fig. 4, it can be noted that as expected, a surface with a high value of ρ doesn't significantly modify the optical gain distribution. On the contrary, we can observe that results for a weak ρ value in Fig. 5 exhibit different conclusions. Indeed, the optical gain PDF is more spread of about 5dB in scenarios (d) and (e) compared to scenario (a) where the body is not taken into account. Those results show that the body has an impact on the optical gain, mostly depending on its ρ value and not on its shape. To discuss performance, we will use in the following a 2D model with the two extreme values respectively = 0.1 and = 0.7.

B. Impact on the channel performance

To establish the performance, we first determine the SNR for the scenarios (a), (b) and (d). In the case of On-Off Keying modulation (OOK) the SNR can be expressed as [START_REF] Ghassemlooy | Opitcal Wireless Communications System and Channel Modelling with MATLAB[END_REF]:

𝑆𝑁𝑅 𝑂𝑂𝐾 = 2(𝑅.𝑃 𝑟-𝑠𝑖𝑔 )² 𝜎 𝑠ℎ𝑜𝑡 2 + 𝜎 𝑡ℎ𝑒𝑟𝑚𝑎𝑙 2 (3)
Where R is the photodetector sensitivity, σ 𝑠ℎ𝑜𝑡 2 represents the shot noise power and σ 𝑡ℎ𝑒𝑟𝑚𝑎𝑙 2 the thermal one. In our case, we consider that σ 𝑡ℎ𝑒𝑟𝑚𝑎𝑙 2 is very small compared to σ 𝑠ℎ𝑜𝑡 2 . So, the noise power is obtained as:

𝜎 𝑠ℎ𝑜𝑡 2 = 2𝑞𝑅(𝑃 𝑟-𝑠𝑖𝑔 + 𝑃 𝑟-𝑖𝑠𝑖 )𝐵 + 2𝑞𝐼 𝐵 𝐼 2 𝐵 ( 4
)
Where q is the electron quantum charge, IB is the photocurrent due to background radiations, I2 is the bandwidth noise factor and B = 1/T is the bandwidth of the signal. We set I2 = 0.562 and IB = 5100µA [START_REF] Moreira | Optical interference produced by artificial light[END_REF]. Pr-sig and Pr-isi are respectively the received power within the symbol duration T and the power linked to Inter Symbol Interference (ISI). They are defined as:

𝑃 𝑟-𝑠𝑖𝑔 = ∫ ∑ ℎ 𝑖 (𝑡)⨂𝑥(𝑡) 𝑁 𝑖=1 𝑇 0 𝑑𝑡 ( 5
)
𝑃 𝑟-𝑖𝑠𝑖 = ∫ ∑ ℎ 𝑖 (𝑡)⨂𝑥(𝑡)𝑑𝑡 𝑁 𝑖=1 ∞ 𝑇 (6)
x(t) is the transmitted optical pulse. The symbol ⨂ denotes the convolution. Also, the average transmitted power Pt is given by:

𝑃 𝑡 = lim 𝑇 ⇾ ∞ 1 2𝑇 ∫ 𝑥(𝑡)𝑑𝑡 +𝑇 -𝑇 (7)
We use the outage probability criterion to discuss the performance considering mobility [18]. It corresponds to the probability that the communication link would not ensure an aimed Quality of Service (QoS) defined here by a targeted SNR value named 𝑆𝑁𝑅 0 related to the BER for OOK modulation by:

𝐵𝐸𝑅 = 1 2 𝑒𝑟𝑓𝑐(√ 𝑆𝑁𝑅 0 2 ) (8)
The outage probability 𝑃 𝑜𝑢𝑡 can be expressed as:

𝑃 𝑜𝑢𝑡 = 𝑃 𝑟 (𝑆𝑁𝑅 < 𝑆𝑁𝑅 0 ) (9)
We have reported in Fig. 6 and Fig. 7 the outage probability for two data rates: Rb1 = 1Mbps and Rb2 = 350Mbps respectively for scenarios (a);(b) and scenarios (a);(d).

First, we can see that to achieve a Pout of 10 -1 , the impact of the body is negligible whatever the data rate and the body reflectivity are.

Secondly, for lower values of Pout, the results show that body impact is more significant when the data rate increases. For example in Fig. 6, for a Pout of 10 -3 , we have determined the corresponding SNR0 value for each scenario.

The difference in terms of SNR is ΔSNR = 0.4dB between scenario (a) and (b) for Rb1 whereas ΔSNR = 21dB for Rb2. In addition, this effect is all the more significant as the reflectivity value is weak: in Fig. 7, ΔSNR = 8dB between scenario (a) and (d) for Rb1 and ΔSNR = 36dB for Rb2. Nevertheless, the body shape and ρ value are not the only parameters having penalizing effect. Indeed, we have already shown in [START_REF] Pergoloni | Coverage optimization of 5G atto-cells for visible light communications access[END_REF] for the same VLC configuration that receiver orientation variations have a strong impact on the performance. But, this had been performed without considering body.

So, in order to determine the impact of orientation variations with a body, we have considered another scenario (f) and compared the results with scenario (b):

-(f) mobile receiver with 2D body model with =0.7 and fixed receiver orientation (θ = 0)

We have reported in Fig. 8 the outage probability for Rb1 and Rb2 for the two scenarios (b) and (f). It can be observed that the achievable performance between both scenarios at Rb1 are similar, meaning that for low data rates, considering a 2D body with only perpendicular receiver orientation is an assumption allowing suitable performance estimation.

However, for high data rates, the impact of the receiver orientation variations cannot be neglected. Actually, performance with Rb2 is highly degraded in scenario (f) compared to (b). For example, for a Pout of 10 -3 , the ΔSNR value between scenario (b) and (f) is quite equal to 30dB. This high value of ΔSNR can be explained by the fact that considering a strong ISI because of high data rates implies cases in some orientations where the received signal is almost blocked by the body and so the optical gain is quite equal to zero. This blocking effect has been observed considering a high body ρ2 value (ρ2 = 0.7) and it is likely that the effect will be even more significant for smaller values.

To investigate body blocking impact for high data rates regarding body ρ value, we thus consider the two extreme values ρ1 and ρ2 and that receiver orientations randomly vary. We have reported in Fig. 9, the evolution of the minimal SNR0 value required for an outage probability of 10 -3 and a wide range of data rate from 1kbps to 500Mbps. The studied scenarios are (a), (b) and (d).

As expected, we can observe that for a high ρ value (scenario (b)) and for data rates lower than 250Mbps, the impact of the body on the SNR0 value is not significant. However, from 300Mbps, the SNR0 value rapidly decreases and so performance is degraded. As we have seen above, this is linked to ISI impact in some orientations. For a weak ρ value in case (d), the SNR0 value is always lower than in case (b) even for low data rates. For example, for Rb = 100Mbps, the SNR0 gap between the two cases is around ΔSNR = 15.3dB. However, until around 300Mbps the requested SNR0 values for the targeted Pout are still high enough to ensure low BER values with OOK modulation. Finally, from 300Mbps we can also remark as previously that the SNR0 values highly decrease. This means that to design an efficient high speed VLC link it is important to model the presence of the body along with receiver orientations to take into account blocking cases.

These results confirm that the choice of the body's ρ value is very important, because regarding its value the body can or cannot have an impact on the channel performance, even at low data rates.

IV. CONCLUSION

The main goal of this paper was to investigate the impact of the body on a VLC channel link. This channel was considered mobile, meaning that positions and orientations of the receiver are randomly varying. Moreover, we have considered that the receiver is worn by a person, whose body is also taken into account. Our study was based on the comparison of two body models: one in 2D and the other one in 3D. We have also considered two extreme ρ values for the body (ρ1 = 0.1 and ρ2 = 0.7).

First, by comparing the VLC gain distributions obtained with the two body models, we have shown that the body shape does not have a significant impact whatever the reflectivity is. Thus, it is not necessary to use a complex model to study the body impact on the VLC link performance.

Performance has been then analyzed in terms of outage probability considering different values of data rates corresponding to different ISI cases. We have shown that to achieve low QoS (outage probability of 10 -1 ) body impact is not significant whatever the data rate or the body reflectivity are.

However, for high QoS (outage probability of 10 -3 ), performance is degraded depending on the body ρ value along with the data rate one.

In addition, we have included the receiver orientation impact in our analysis. The results have first shown that for low data rates, considering only perpendicular receiver orientation is an assumption allowing suitable performance estimation. But, for high data rates, it is necessary to model receiver orientation variations along with the body features.

Finally, we have evaluated the minimal SNR value required for a given outage probability for a wide range of data rates (1kbps to 450Mbps).

We have seen that for a high body ρ value and data rates lower than 300Mbps, body impact is not significant, whereas from higher data rates performance can be highly and rapidly degraded due to body blocking impact. For weak body reflectivity values, the results are similar but more penalizing.

Those results have been obtained for a specific environment, but they highlight that the impact of the body on a mobile VLC channel depends on the data rate and on the ρ value of the body. This is a first study which has to be completed in particular by studying the impact on the performance of the distance between the receiver and the body.

Fig. 1 .

 1 Fig. 1. VLC scenario configuration

2 .

 2 Geometrical models of the body in (a) 2D, (b) 3DThe second question concerns the reflectance properties of the environment. Each facet composing the body (2D or 3D) is characterized by a Lambertian Bidirectional Reflectance Distribution Function (BRDF)[18] which is a purely diffuse model. This one is defined by a single parameter ρ,

Fig. 3 .

 3 Fig. 3. Evolution of some ρ values according to the light wavelength On one side, white human skin ρ values go from 0.16 to 0.86 when dark skin values go from 0.04 to 0.35. The corresponding mean values are 0.51 and 0.19 respectively. On the other side, dark blue jeans ρ values go from 0.5 to 0.75 which correspond to a mean value of 0.62, while faded jeans ρ values go from 0.55 to 0.72 with a mean value of 0.63.In conclusion, to take into account the whole range of ρ values connected to the different human skins and clothes, and so to obtain the most representative conclusions, all light simulations presented in this paper are obtained concerning two extreme mean values of , respectively = 0.1 and = 0.7. Finally, the reflectance properties of the rest of the environment (walls, ceiling and floor) are presented in the TABLE II. They have been fixed according to[START_REF] Lee | Indoor Channel Characteristics for Visible Light Communications[END_REF].

  a): mobile receiver without considering body -(b): mobile receiver and 2D body model with =0.7 -(c): mobile receiver and 3D body model with =0.7 -(d): mobile receiver and 2D body model with =0.1 -(e): mobile receiver and 3D body model with =0.1

Fig. 4 .

 4 Fig. 4. Probability density function of the optical gain for scenarios (a), (b) and (c)

Fig. 5 .

 5 Fig. 5. Probability density function of the optical gain for scenarios (a), (d) and (e)

Fig. 6 . 7 Fig. 7 .

 677 Fig. 6. Outage probability for Rb1 =1Mbps and Rb2 =350Mbps considering 2D body and ρ2 = 0.7

Fig. 8 .

 8 Fig. 8. Outage probability for scenarios (b) and (f) for Rb1 =1Mbps and Rb2 =350Mbps

Fig. 9 .

 9 Fig. 9. Minimal SNR value at Pout = 10 -3 for data rates from 1kbps to 450 Mbps for scenarios (a), (b) and (d)

TABLE II .

 II Reflectance properties of the room

	Geometrical elements	BRDF	ρ values
	Walls/Ceiling	Lambertian	0.73
	Floor	Lambertian	0.35
		III. HUMAN BODY IMPACT