Cette Thèse

Philippe Schnoebelen

Albert Cohen

Paul Feautrier

Laure Gonnord

Jacqueline Altimira

Mehdi Amini

Corinne Ancourt

Karim Barkati

Pierre Beauguitte

Danielle Bolan

Catherine Le Caër

Fabien Coelho

Laurent Daverio

Wajdi Farhani

Imré Frotier De La Messelière

Emilio Jesús

Gallego Arias

Florian Gouin

Pierre Guillou

Pierre Jou- Velot

Dounia Khaldi

Nelson Lossing

Robert Mahl

Claire Medrala

Amira Mensi

Benoît Pin

Antoniu Pop

José Afonso Sanchez

Ronan Saillard

Arnaud Spiwack

Claude Tadonki

Karel De Vogeleer

Haisheng Wang

Pierre Wargnier

interpretation is a theory to approximate computer program semantics, based on monotonic functions on ordered sets, in particular lattices. It can be defined as a partial program execution to obtain information on its semantic (for instance, control structure, data flow) without having to fully process it.

The main function of abstract interpretation is static analysis, automatic extraction of information about possible executions of a program. This chapter defines basic notions of abstract interpretation, that are used later in the document.

Introduction (English version follows.)

Un système critique pour la vie ou pour la sécurité est un système dont la défaillance ou le dysfonctionnement peut grièvement blesser ou tuer des gens, détruire ou faire subir des dégâts importants à des équipements, ou causer des dommages environnementaux. Les méthodes et les outils de l'ingénierie de la sécurité permettent d'appréhender les risques de cette nature. Un système critique est conçu pour perdre moins d'une vie par milliard (10 9 ) d'heures de fonctionnement [START_REF]System Design and Analysis[END_REF]. Sa conception fait général appel à l'évaluation probabiliste des risques, une méthode qui combine l'analyse des modes de défaillance, de leurs effets et de leur criticité (AMDEC) avec l'analyse par arbre de défaillance. Aujourd'hui, les systèmes de sécurité critiques sont de plus en plus souvent informatisés.

L'ingénierie logicielle pour les systèmes critiques est particulièrement difficile. Elle repose sur plusieurs éléments. Le premier est le génie et la gestion des procédés. Le second est le choix des outils et d'un environnement appropriés pour le système. Cela permet au développeur de tester efficacement le système par simulation et d'observer son efficacité. Le troisième consiste à se soumettre à toutes les exigences légales et réglementaires de sécurité, telles que les exigences de la FAA pour l'aviation (DO-178B [START_REF]Software Considerations in Airborne Systems and Equipment Certification[END_REF][START_REF] Johnson | DO-178B, Software Considerations in Airborne Systems and Equipment Certification[END_REF]). La définition d'une norme selon laquelle doit être développé un système contraint les concepteurs à respecter ces exigences. L'industrie aéronautique a réussi à produire des normes standard pour la production de logiciels d'aviation critiques. Des normes similaires existent dans l'industrie automobile (MISRA C [START_REF] Misra | Guidelines for the Use of the C Language in Vehicle Based Software[END_REF], ISO 26262 [START_REF] Kucharski | ISO 26262 Software Compliance: Achieving Functional Safety in the Automotive Industry[END_REF]), médicale (IEC 62304 [START_REF] Jordan | Standard IEC 62304 -Medical Device Software -Software Lifecycle Processes[END_REF]) et nucléaire (IEC 61513). L'approche standard consiste à soigneusement coder, inspecter, documenter, tester, vérifier et analyser le système.

Correctement appliquée, cette approche permet d'atteindre un très haut niveau de sécurité : des logiciels critiques ont été utilisés durant des décennies dans l'industrie aérospatiale sans causer de désastre majeur. Malheureusement, ces procédés de développements sont excessivement lourds et coûteux : la vérification et la validation des logiciels critiques représentent ainsi plus de 50% du coût de son développement [WDD + 12]. En raison de leur complexité et de leur coût, ces contraintes sont parfois mal respectées, ce qui conduit à des bogues avec des conséquences parfois dramatiques. Un exemple récent est le défaut du contrôleur électronique de la commande d'accélération des véhicules Toyota Camry, qui a causé plusieurs morts [START_REF] Yoshida | Toyota Case: Single Bit Flip That Killed[END_REF]. De plus, les tests seuls ne sont pas suffisants : étant donné le nombre considérable de scénarios de tests, il est impossible de couvrir tous les cas et il n'y a aucune garantie qu'un bogue n'ait pas pu échapper aux tests.

Ces éléments expliquent l'intérêt académique et industriel actuel pour les méthodes formelles, qui fournissent des garanties mathématiques de correction, avec différents niveaux d'automatisation. Une approche possible est de certifier un système de production et un compilateur, puis de générer le code du système à partir des spécifications. Une autre consiste à utiliser des méthodes formelles pour générer les preuves qui garantissent que le code répond à sa spécification.

Contenu et contributions

Ce manuscrit est organisé en deux parties, chacune d'elles correspondant à une contribution distincte. La première partie porte sur les systèmes de contrôle. Un système de contrôle est un dispositif ou un ensemble de dispositifs qui gère, commande, dirige ou régule le comportement d'autres dispositifs ou systèmes (chapitre 1). Des systèmes de contrôle industriels sont utilisés en production, pour contrôler des équipements ou des machines. La propriété formelle qui nous intéresse est la stabilité : en termes simples, c'est la propriété selon laquelle toutes les solutions d'un système dynamique qui commencent à proximité d'un point d'équilibre restent par la suite dans son voisinage. La théorie du contrôle énonce qu'un contrôleur linéaire est stable si et seulement si il admet une invariant quadratique -géométriquement, un ellipsoïde, appelé fonction quadratique de Lyapunov (chapitre 2), et propose des méthodes pour trouver un invariant de cette nature à partir de la description algébrique du système. Nous présentons un framework formel, générique, pour traduire des preuves de stabilité de Lyapunov codées en MATLAB avec des nombres réels vers du code C de bas niveau en arithmétique machine (chapitre 3) ; ainsi qu'un outil pour générer automatiquement les invariants de stabilité sur de l'arithmétique à virgule flottante (chapitre 4).

Il arrive aussi qu'il soit nécessaire de reconstruire automatiquement les invariants sur un programme donné. C'est pourquoi la seconde partie porte sur l'analyse des relations linéaires (chapitres 5, 6, 7), un framework d'interprétation abstraite couramment utilisé, consacré à la génération automatique d'invariants sous forme d'inégalités linéaires sur les variables numériques d'un programme. Nous présentons ALICe, un ensemble d'outils pour comparer les techniques de calcul automatique d'invariants affines sur les boucles (chapitre 8). Il s'accompagne d'un benchmark construit à partir de 102 cas de test provenant de la littérature sur les problèmes d'invariants et de terminaison de boucles, et s'interface avec trois outils d'analyse de programmes, utilisant des algorithmes différents : Aspic, iscc et PIPS. Les résultats expérimentaux montrent l'importance du codage du modèle et les performances en retrait de PIPS en présence de boucles concurrentes dans le programme analysé. Pour étudier ces problèmes, nous utilisons deux techniques de restructuration de programme (chapitre 9), prouvées correctes en Coq, et nous présentons plusieurs améliorations apportées aux algorithmes de PIPS pour le calcul d'invariant (chapitre 10).

* * *

A life-critical system or safety-critical system is a system whose failure or malfunction may result in death or serious injury to people, loss of or severe damage to equipment or environmental harm. Risks of this sort are usually managed with the methods and tools of safety engineering. A critical system is designed to lose less than one life per billion (10 9 ) hours of operation [START_REF]System Design and Analysis[END_REF]. Typical design methods include probabilistic risk assessment, a method that combines failure mode and effects analysis (FMEA) with fault-tree analysis. Safety-critical systems are increasingly computer-based.

Software engineering for critical systems is particularly difficult. There are three aspects which can be carried out to support this discipline. First is process engineering and management. Secondly, selecting the appropriate tools and environment for the system. This allows the system developer to effectively test the system by emulation and observe its effectiveness. Thirdly, address any legal and regulatory requirements, such as FAA requirements for aviation (DO-178B [START_REF]Software Considerations in Airborne Systems and Equipment Certification[END_REF][START_REF] Johnson | DO-178B, Software Considerations in Airborne Systems and Equipment Certification[END_REF]). By setting a standard for which a system is required to be developed under, it compels the designers to stick to the requirements. The avionics industry has succeeded in producing standard methods for producing critical avionics software. Similar standards exist for automotive (MISRA C [START_REF] Misra | Guidelines for the Use of the C Language in Vehicle Based Software[END_REF], ISO 26262 [START_REF] Kucharski | ISO 26262 Software Compliance: Achieving Functional Safety in the Automotive Industry[END_REF]), medical (IEC 62304 [START_REF] Jordan | Standard IEC 62304 -Medical Device Software -Software Lifecycle Processes[END_REF]) and nuclear (IEC 61513) industries. The standard approach is to carefully code, inspect, document, test, verify and analyze the system.

Correctly applied, this approach allows to achieve a very high level of safety: critical software has been used for decades in the aerospace industry without any major disaster. Unfortunately, such development processes are excessively heavy and expensive: verification and validation for critical software represent more than 50% of software development costs [WDD + 12]. Due to their complexity and cost, these constraints are sometimes poorly respected, which leads to bugs with dramatic consequences. A recent example is the electronic throttle control defect linked to unintended acceleration in Toyota Camry vehicles [START_REF] Yoshida | Toyota Case: Single Bit Flip That Killed[END_REF] that caused deaths. Furthermore, testing alone is not enough: considering the huge number of test scenarios, exhaustive testing is impossible and there is not any guarantee that the processed test cases did not fail to reveal a specific bug.

These elements explain the current academic and industrial interest in formal methods; they provide mathematical insurances of correctness, with various levels of automation. One possible approach is to certify a production system, a compiler, and then generate the system's code from specifications. Another one is to use formal methods to generate proofs that the code meets requirements.

This thesis focuses on model checking and abstract interpretation (with a secondary use of proof assistants), an efficient family of formal methods to exhaustively and automatically check whether a model meets a given specification, given in the form of numerical properties. More precisely, these techniques compute invariants on a program, i.e. conditions on its states that are true during its execution, or during some portion of it. They are logical assertions that hold true during a certain phase of execution, and they are used to check whether the specification holds. Therefore, generating accurate invariants is a key element to prove the correctness of a program with respect to its specification.

Contents and Contributions

This manuscript in organized in two parts, each of them corresponding to distinct contributions. The first part is about control systems. A control system is a device, or set of devices, that manages, commands, directs or regulates the behavior of other devices or systems (Chapter 1). Industrial control systems are used in industrial production for controlling an equipment or a machine. The formal property we are interested in is stability: in simple terms, that all solutions of a dynamical system that start out near an equilibrium point stay in its neighborhood forever. Control theory states that a linear controller is stable if and only if it admits a quadratic invariant -geometrically speaking, an ellipsoid, called quadratic Lyapunov functions (Chapter 2), and offers practical ways to find such an invariant from the algebraic description of the system. We present a formalized, generic framework to translate Lyapunov stability proofs on MATLAB real-valued code to low-level C code using machine arithmetic (Chapter 3), and a tool that automatically generates stability invariants in the case of floating-point arithmetic (Chapter 4).

It is also sometimes necessary to automatically reconstruct the invariants on given implementations. This is why the second part focuses on linear relation analysis (Chapters 5, 6, 7), a widespread abstract interpretation framework devoted to the automatic discovery of invariant linear inequalities on the numerical variables of a program. We present ALICe, a toolset to compare automatic computation techniques of affine loop scalar invariants (Chapter 8). It comes with a benchmark that we built using 102 test cases found in the loop invariant and termination bibliography, and interfaces with three analysis programs, that rely on different techniques: Aspic, iscc and PIPS. Experimental results show the importance of model coding and the poor performances of PIPS on concurrent loops. To tackle these issues, we use two restructuring techniques whose correctness is proved in Coq (Chapter 9), and present various improvements made to PIPS algorithms for invariant computation (Chapter 10).

Part I Linear Systems

Chapter 1

Control-Command Critical Systems

Un système dynymatique est un concept en mathématiques où une règle fixe décrit la dépendance temporelle d'un point dans un espace géométrique. Les modèles mathématiques qui décrivent le mouvement du balancier d'une horloge, l'écoulement de l'eau dans un tuyau ou le nombre de poisson chaque été dans un lac sont des exemples de systèmes dynamiques.

À tout instant, le système possède un état défini par un ensemble de nombres réels (un vecteur), qui peut être représenté par un point dans l'espace d'état approprié. Une légere altération de l'état corespond à un légere altération de ces valeurs. La règle d'évolution du système dynamique est une règle fixe qui décrit quel état futur peut découler de l'état actuel. Cette règle est déterministe : dans un intervalle de temps donné, l'état actuel ne peut mener qu'à un seul état futur.

La théorie du contrôle (ou automatique) est la branche interdisciplinaire de l'ingénierie et des mathématiques qui traite du comportement des systèmes dynamiques avec des entrées. L'entrée extérieure d'un système est appelée la référence (ou consigne). Quand une ou plusieurs variables de sortie d'un système doivent suivre une certaine consigne au cours du temps, le contrôleur manipule l'entrée du système de manière à obtenir l'état souhaité en sortie. Il s'agit d'un dispositif, historiquement mécanique, hydraulique, pneumatique et/ou électronique, plus récemment un microprocesseur ou un microcontrôleur, qui surveille et modifie les conditions de fonctionnement du système dynamique. La régulation de la température, de la pression, du débit ou de la vitesse sont des applications typiques des contrôleurs. Le régulateur à boules sur la figure 1.1 est l'un des premiers exemples de contrôleur mécanique, que l'on rencontre souvent sur les moteurs à vapeur : il contrôle la vitesse du moteur en régulant la quantité de vapeur admise, afin de maintenir une vitesse quasi-constante, indépendamment de la charge ou de la fourniture en combustible.

Les sujets abordés en théorie du contrôle incluent la stabilité, la contrôlabilité et l'observabilité, la spécification du contrôle, l'idenfication et la robustesse du modèle. Parmi ceux-ci, un objectif courant est de calculer les solutions pour les actions de correction du contrôleur qui garantissent la stabilité du système, c'est-à-dire que le système restera sur le point qui aura été défini au lieu d'osciller autour de ce point. Une définition précise de la stabilité est donnée dans la section 1.3. Bien que le champ d'application principal de l'automatique soit l'ingénierie des systèmes de contrôle, qui traite de la conception de systèmes de contrôle industriels, d'autres applications dépassent largement de ce cadre. L'automatique peut être utilisée pour étudier tout système avec rétroaction et trouve des applications dans des domaines aussi divers que la physiologie, l'életronique, la climatologie, les écosystèmes, la navigation, les réseaux de neurones ou encore la génétique.

* * *

A dynamical system is a concept in mathematics where a fixed rule describes the time dependency of a point in a geometrical space. Examples include the mathematical models that describe the swinging of a clock pendulum, the flow of water in a pipe, or the number of fish each springtime in a lake.

At any given time a dynamical system has a state defined by a set of real numbers (a vector) that can be represented by a point in an appropriate state space (a geometrical manifold). Small changes in the state of the system are represented by small changes in the numbers. The evolution rule of the dynamical system is a fixed rule that describes which future state can follow from the current state. The rule is deterministic; in other words, for a given time interval only one future state follows from the current state.

Control theory is the interdisciplinary branch of engineering and mathematics that deals with the behavior of dynamical systems with inputs. The external input of a system is called the reference. When one or more output variables of a system need to follow a certain reference over time, a controller manipulates the input to a system to obtain the desired output state of the system. It is a device, historically using mechanical, hydraulic, pneumatic or electronic techniques often in combination, but more recently in the form of a microprocessor or microcontroller, which monitors and physically alters the operating conditions of a given dynamical system. Typical applications of controllers are to hold settings for temperature, pressure, flow or speed. The centrifugal governor in Figure 1.1 is an early example of mechanical controller, commonly seen on steam engines, that controls the speed of an engine by regulating the amount of steam admitted, so as to maintain a near-constant speed, irrespective of the load or fuel-supply conditions. Topics in control theory include stability, controllability and observability, control specification, model identification and robustness. Among them, an usual objective is to calculate solutions for a proper corrective action from the controller that result in system stability, that is, the system will hold the set point and not oscillate around it. A precise definition of system stability is given in Section 1. 3.

Although a major application of control theory is control systems engineering, which deals with the design of process control systems for industry, other applications range far beyond this. As the general theory of feedback systems, control theory is useful wherever feedback occurs; a few examples are in physiology, electronics, climate modeling, machine design, ecosystems, navigation, neural networks, and gene expression.

Control-Command Systems

A control-command systems is an autonomous, self-controlled system that includes:

• The regulated parameter.

• The controller.

• Sensors to feed input data to the controller.

• Actuators through which the controller acts on the system to reach the desired state.

A complex example of control-command system is a modern airliner: parameters like the plane's altitude and speed are regulated by an on-board computer, using a set of sensors (speedometer, altimeter…), and adjusting the flight path by playing on actuators like engine thrust and flight control surfaces.

Plant Modeling

To design a control-command system, control engineers need to know precisely the effect actuators have on the system. For a physical system, it is required to design a model of the system behavior in its environment, using laws of mechanics and physics. This model is called a plant.

Linear systems are a particular class of control-command systems that can be described by linear equations (linear differential equations or difference equations). They are much easier to study than the general, nonlinear cases. Therefore, a nonlinear system is often linearized around an equilibrium point or trajectory (when possible), to obtain a linear system that accurately represents the nonlinear system in the neighborhood of this equilibrium point or trajectory.

From a mathematical perspective, the internal state x p ∈ R n of a linear plant (a vector containing the system characteristics in the environment: for a plane, it could be the speed, altitude, etc.) is ruled by a differential equation ẋp = A p x p + B p u p where A p ∈ R n×n and B p ∈ R n×p are constant, real-valued matrices, and u p ∈ R p is the vector of system inputs (speed, altitude…).

In the setting defined by the internal state x p and the input values u p , the output y p (interpreted as the system action) follows the equation

y p = C p x p + D p u p .
where C p , D p are constant, real-valued matrices.

Controller Design

Control theory offers various methods to develop a controller corresponding to a particular plant model. In the usual case of a linear controller, the controller internal state x c follows the equation

ẋc = A c x c + B c u c (1.1)
where u c is the controller input, that typically comes from input sensors y p (speed and altitude measurements) and commands y d (cockpit controls), and equals the output of the physical system, minus the command (y py d ). The controller computes the output

y c = C c x c + D c u c (1.2)
that is applied to actuators and is equal to the plant input u p (Figure 1.3).

Controller Implementation

Once the controller is designed, they are several ways to implement it:

• Hardware implementation: either with a purely mechanic implementation -as is the centrifugal governor in Figure 1.1 -or with analog electronics.

• Software implementation: running a program on a numerical controller, as on most modern devices.

In this document, we focus exclusively on software controllers. To implement a controller, it is necessary to discretize its internal state equation ( 1.1), in order to be able to recompute its internal state x c at a fixed frequency:

x c k+1 = A ′ c x c k + B ′ c u c k . (1.3)
Typical frequencies for embedded systems are about 100 Hz.

Eventually, the discrete controller equations (1. [START_REF] Feron | From Control Systems to Control Software[END_REF] This is a specific, numerical code, consisting of a main loop and nested computations; it looks very different from the source code of a compiler for example. Of course, the true controller of a large scale system is something much more complex: real-scale software may require hundreds of thousands of lines of C code (about 170,000 source lines of generated C code for an Airbus A340, and 760,000 lines for an Airbus A380).

Open-Loop and Closed-Loop Controllers

An open-loop controller (also called a non-feedback controller), is a type of controller that computes its input into a system using only the current state and its model of the system.

A characteristic of the open-loop controller is that it does not use feedback to determine if its output has achieved the desired goal of the input. This means that the system does not observe the output of the processes that it is controlling. Consequently, a true open-loop system may not compensate for disturbances in the system.

An open-loop controller is shown in Figure 1. 2.

Controller (state x c ) u c y c Figure 1.2: Open-loop controller
To overcome the limitations of the open-loop controller, control theory introduces feedback. A closed-loop controller uses feedback to control states or outputs of a dynamical system. Its name comes from the information path in the system: process inputs (e.g., fuel flow in the jet engines) have an effect on the process outputs (e.g., engine thrust), which is measured with sensors and processed by the controller; the result (the control signal) is fed back as input to the process, closing the loop.

Closed-loop controllers have the following advantages over open-loop controllers:

• Disturbance rejection.

• Guaranteed performance even with model uncertainties, when the model structure does not match perfectly the real process and the model parameters are not exact, unstable processes can still be stabilized.

• Reduced sensitivity to parameter variations.

• Improved reference tracking performance.

A closed-loop controller is shown in Figure 1. 3.

Plant (state x p ) Controller (state x c ) y d y c u p y p u c

Stability Proofs

Stability theory addresses the stability of solutions of differential equations and of trajectories of dynamical systems under small perturbations of initial conditions.

When focusing on stability, two types of properties are of major interest: closed-loop stability and open-loop stability.

• Closed-loop stability expresses that the controller manages to stabilize the plant, assuming that the physical system is correctly described by the model: if the command y d is bounded, then the internal states x p and x c and the outputs y p and y c of the plant and the controller remain bounded.

• Open-loop stability means that all variables in the controller are bounded, assuming its inputs are bounded: if the controller input u c is bounded, its internal state x c and its output y c remain bounded. In [START_REF] Feron | From Control Systems to Control Software[END_REF], Feron investigates how control-system domain knowledge and, in particular, Lyapunovtheoretic proofs of stability for the high-level theoretical modeling, can be migrated towards computer-readable and verifiable certificates of control software behavior by relying on Floyd's and Hoare's proof system [START_REF] Peled | Software Reliability Methods[END_REF], applied to MATLAB pseudo-code (see Sections 2.2 and 2.3).

Feron presents both open-loop and closed-loop implementations of his controller, i.e. depending whether the system action on the environment and its feedback are modeled or not. But errors resulting from the use of floating-point arithmetic are not considered.

In this part, we address this issue by presenting an automatic approach to translate Lyapunovtheoretic stability proof invariants on pseudo-code with real arithmetic, as provided in Feron's article, to similar invariants on machine code that take into account rounding errors introduced by floating-point arithmetic. We use them to verify whether the stability conditions still hold, in which case system stability with floating-point numbers is established. This approach is based on a generic theoretical framework, the soundness of which is proved in Coq, described further in this document, and is implemented as an automatic tool.

This part of the dissertation is organized as follows. In the remainder of this chapter, we describe the second-order dynamical system example used by Feron in [START_REF] Feron | From Control Systems to Control Software[END_REF], with the corresponding controller. Next, Feron's analysis of the open-loop controller with real numbers is presented.

In the next chapters, we introduce our generic translation scheme (Chapter 3, page 29) and its implementation as a Python library, and use it thereafter to rewrite Feron's analysis with floating-point arithmetic (Chapter 4, page 37). Finally, the case of the closed-loop system is handled. This part concludes with a discussion on the generality of this approach.

Lyapunov Stability Theory

Various types of stability may be discussed for the solutions of differential equations describing dynamical systems. The most important type is the one concerning the stability of solutions near to a point of equilibrium. This may be discussed by the theory of Lyapunov. In simple terms, if the solutions that start out near an equilibrium point x e stay near x e forever, then

x e is Lyapunov stable. More strongly, if x e is Lyapunov stable and all solutions that start out near x e converge to x e , then x e is asymptotically stable.

Lyapunov functions are scalar functions that may be used to prove the stability of an equilibrium. The existence of Lyapunov functions is a necessary and sufficient condition to establish the stability of a discrete, time-invariant closed system, meaning that its state variable will remain bounded during the execution [START_REF] Haddad | Nonlinear Dynamical Systems and Control: A Lyapunov-Based Approach[END_REF].

Informally, a Lyapunov function is a function that takes positive values everywhere except at the equilibrium in question, and decreases (or is non-increasing) along every trajectory of the control system's differential equation f . It is a function V : R n → R that satisfies the following properties:

V (0) = 0 and ∀x ∈ R n \ {0}, V (x) > 0 and lim x →∞ V (x) = ∞ (2.1) ∀x ∈ R n , V (f (x)) -V (x) ≤ 0. (2.2)
In the case of discretized linear systems, i.e. systems whose internal state is defined according to Equation (1.3):

x k+1 = Ax k + Bu k and assuming that u k ∞ ≤ 1, (2.3)
it can be shown that the Lyapunov function V can be taken as a quadratic form of the state variable of the system:

V (x) = x ⊺ P x with P ∈ R n×n a symmetric matrix such that P ≻ 0 (2.4)

A ⊺ P A -P ≺ 0 (2.5)
where the notation P ≻ 0 means that P is a positive-definite matrix, i.e., for all non-zero vector x, ensuring that x ⊺ P x > 0. Equation (2.3) expresses that values coming from input sensors are usually bound into a given range. The bound value 1 is chosen with no loss of generality since the matrix B can always be altered to account for different bounds.

Motivating Example by Feron

We consider the first system described in the article of Feron [START_REF] Feron | From Control Systems to Control Software[END_REF]. It is a dynamical system composed of a single mass and a single spring shown in Figure 2. 1.

The position input y of the mass is available for feedback control. The signal y d is the reference signal, that is, the desired position to be followed by the mass. We do not show (and Feron does not) that the distance |y d -y| is bounded. y d has to be understood as a command. A discrete-time MATLAB implementation of the controller, using real numbers, is provided in [START_REF] Feron | From Control Systems to Control Software[END_REF]. The source code is shown in Listing 2. 1. The process to produce this code from the system modeling is covered in detail in Feron's paper. In this code, the skip statement is a null operation: when it is executed, nothing happens.

It is useful as a placeholder when a statement is required syntactically to be surrounded by invariants, but no code needs to be executed.

Apart from the mechanical system state observation y and the desired system output y d , variables in this code are:

• x c = xc 1 xc 2 ∈ R 2 is the discrete-time controller state; • y c ∈ [-1, 1]
is the bounded output tracking error, i.e. the input (yy d ) is passed through a saturation function to avoid variable overflow in the controller;

• u ∈ R is the mechanical system input (and the controller output at the same time), i.e. the action to be performed according to the controller.

Constants A c , B c , C c and D c are the discrete-time controller state, input, output and feedthrough matrices. The commands send and receive are basically I/O: they respectively send and receive data given in the commands' first argument through a specific channel given by the commands' second argument.

Open-Loop Stability Proof with Reals

The stability proof of this system relies on Lyapunov theory. As seen in Section 1.3 and Section 2.1, a system is Lyapunov stable if all states x c reachable from an initial starting state belonging to a bounded neighborhood V of an equilibrium point x e remain in V . Lyapunov theory provides constraints that must be satisfied by such a V . On linear systems, they are equations that can be solved using linear matrix inequalities [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. Commonly, V is an ellipsoid.

In our case, to prove Lyapunov stability, we need to show that at any time, x c belongs to the set E P chosen by Feron: This set is the full ellipse shown in Figure 2.2, centered around 0 and slightly slanted:

E P = {x ∈ R 2 | x ⊺ • P • x ≤ 1}, P = 10
x c ∈ E P ⇐⇒ 0.6742x c 2 1 + 0.0856x c1 x c2 + 2.4651x c 2 2 ≤ 1000.
x c1 Most of the steps of the proof rely on algebraic arguments. For example, the invariant loosening on Line 9:

skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) ⊺ -1 , % u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -
% x c ∈ E P , y 2 c ≤ 1 9 skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 means x c ∈ E P ∧ y 2 c ≤ 1 =⇒ ( xc yc ) ∈ E Qµ ,
with Q µ = µP 0 2×1 0 1×2 1-µ and µ = 0.9991. The correctness of this assertion stems from the fact that, given any value of µ ∈ [0, 1], the domain E Qµ is a solid ellipsoid, centered around 0, and whose intersection with the planes

y c = 1 and y c = -1 is equal to E P . Consequently, the solid bounded cylinder C = {( xc yc ) | x c ∈ E P ∧ y 2 c ≤ 1} is included within E Qµ (Figure 2.3). E Qµ E P C x c2 y c 1 -1 x c1 Figure 2.3: Inclusion of E P within E Qµ
The next two invariants

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) ⊺ -1 , % u 2 ≤ (C c D c ) • Q -1 µ • (C c D c
) -1 also rely on similar algebraic arguments and theorems. Other invariants are trivial. Finally, only the very last loosening

% x c ∈ E R skip; % x c ∈ E P i.
e. E R ⊂ E P , that closes the loop, is not purely algebraic since its validity relies on the numerical parameters A c , B c , C c , D c . This assertion needs to be checked to ensure the correctness of the proof statements. This can be done either numerically, or algebraically for at most two-dimensional systems like this one.

We have checked this proof using Mathematica 8 [START_REF]Wolfram Mathematica[END_REF]. Our proof notebook is available online †1 . †1 Mathematica source file is available at: http://www.cri.mines-paristech.fr/people/maisonneuve/ lyafloat/resources/lyafloat_stability.nb, and the corresponding PDF file at: http://www.cri. mines-paristech.fr/people/maisonneuve/lyafloat/resources/lyafloat_stability.pdf.

Chapter 3

Proof Scheme with Limited-Precision Numbers

Nous voulons vérifier que la preuve de stabilité est toujours valide sur un dispositif de commande réaliste, en utilisant des nombres à précision limitée (par exemple, mais pas nécessairement, des nombres à virgule flottante). Avec une arithmétique en précision limitée, les valeurs des constantes sont légèrement modifiées et les calculs sont susceptibles de produire des erreurs d'arrondi. D'un point de vue théorique, il est impossible de passer de nombres réels à des nombres de précision limitée sans affecter le comportement du contrôleur. Ainsi, la preuve de la stabilité ne peut pas être conservée dans le cas général. Par ailleurs, le programme peut encore être stable si les erreurs d'arrondi sont suffisamment petites et que l'inclusion finale E R ⊂ E P leur laisse assez de place. Nous étudions comment les invariants de preuve peuvent être modifiés pour être appliqués à une sémantique en précision limitée. L'objectif est d'obtenir, à partir du schéma de preuve en nombres réels, un schéma de preuve qui convienne à une arithmétique en précision limitée et dont la correction, non garantie, pourra être vérifiée aussi aisément que dans la preuve d'origine.

* * *

We would like to check that the stability proof still holds on a realistic controller device, using limited-precision numbers (for instance, but not necessarily, floating-point numbers). When using limited-precision arithmetic, the values of constants are slightly altered and calculations are likely to produce rounding errors.

In absolute terms, it is impossible to switch from real to limited-precision numbers without affecting the behavior of the controller. Thus, the stability proof cannot be preserved in the general case. On the other hand, the program can still be stable if rounding errors are small enough and the final inclusion E R ⊂ E P leaves enough room. We study how proof invariants can be tweaked so that they apply to a limited-precision semantic. Our goal is to derive from the proof scheme for real numbers a proof scheme suited for limited-precision arithmetic, whose correctness, although not guaranteed, can be checked as easily as in the original proof.

Formalism

In this section, an abstract formalism is defined to handle a program and its proof, with both real and limited-precision arithmetic. This formalism is implemented in Coq.

Let X be the (finite) set of variables in the program.

% d i % d ′ = p(d, i)

Values

Let F ⊂ R be the target set of limited-precision representations of real numbers (e.g., floatingpoint numbers). We assume that each real number c ∈ R has a finite-precision representation c ∈ F and that Equation (3.1) is satisfied:

c ∈ F =⇒ c = c. (3.1)
Throughout the rest of the document, we will use the uniform notation K to refer to either the set R or F.

Functions Symbols

We consider a set of function symbols F R on real numbers. Each real function symbol f ∈ F R is associated to a finite-precision counterpart f ∈ F F . Functions symbols in F R and F F will be simply referred as functions when there will be no ambiguity.

Each function symbol f in F R or F F has an arity n ∈ N and can be evaluated into a corresponding mathematical function:

f ∈ F K =⇒ eval(f ) ∈ K n → K.
Functions f and f have the same arity: arity f = arity f .

As for now, no other assumption is made as to the behavior of f compared to f . In particular, there is no guarantee that, applied to the same inputs, they will return the same value.

Notice that function symbols are not identified to their evaluations. The reason is that functions with similar evaluations on real arithmetic may behave differently on a limited-precision paradigm, e.g. function f : x → 2 100 + x -2 100 compared to the identity function Id : x → x: eval(f ) = eval(Id) but eval( f ) = eval(Id).

Valuations

A valuation v on K is a function that maps variables in X to values in K. The set of valuations on

K is noted Val K . Val K ∋ v : X → K.
Notice that Val F is a subset of Val R .

Each valuation v in Val R is associated to a valuation v in Val F , defined as follows:

v : x ∈ X -→ v(x).
Using this definition, it can easily be shown using Equation (3.1) that for any valuation

v ∈ Val R , v ∈ Val F =⇒ v = v.
Let v 1 and v 2 be two valuations in Val R . The distance dist between v 1 and v 2 is the valuation that maps a variable in X to the distance in absolute value between its associated values in v 1 and v 2 :

Val R ∋ dist(v 1 , v 2 ) = x ∈ X -→ |v 1 (x) -v 2 (x)|.
We say that v 1 is lower than v 2 , and note

v 1 ≤ v 2 , if ∀x ∈ X, v 1 (x) ≤ v 2 (x).
Finally, the sum of the two valuations, noted v 1 + v 2 , is the valuation defined as follows:

v 1 + v 2 : x ∈ X -→ v 1 (x) + v 2 (x).

Domains (Invariants)

d = {v | v ∈ d}.
As for valuations in F, it can be shown that

d ∈ D F =⇒ d = d.
We also define an operation to overapproximate a domain with respect to a valuation (intuitively, an error valuation ). The extension of a domain d ∈ D R using a positive valuation v ∈ Val R , noted d ⊕ v, is defined as follows:

d ⊕ v = {v 2 ∈ Val R | v 1 ∈ d ∧ dist(v 1 , v 2 ) ≤ v}.

Expressions

We consider a very simple language with variables, constants and function calls. This language is expressive enough to represent most of control-command softwares, which corresponds to our needs. The expressions E K of this language are constructed as follows:

E K ∋ e ::= x ∈ X variable | c ∈ K constant | f (e 1 , . . . , e n ) function call with f ∈ F K , arity f = n and ∀i ∈ [1, n], e i ∈ E K
An expression e ∈ K can be evaluated in the environment described by the valuation v ∈ Val K .

The result is a value in K.

eval(e, v)

=          v(x) if e = x ∈ X c if e = c ∈ K eval(f )(c 1 , . . . , c n ) if e = f (e 1 , . . . , e n ) with ∀i ∈ [1, n], c i = eval(e i , v)
where the notation eval is overloaded to handle expressions.

Instructions

An instruction is the affectation of an expression value to a variable. The set of instructions is noted I K :

I K = X × E K We note x := e the instruction (x, e) ∈ I K .
As for expressions, an instruction (x := e) ∈ I K can be evaluated in an environment v ∈ Val K . The result is a valuation in Val K .

eval(x

:= e, v) = y ∈ X -→ eval(e, v) if x = y v(y) if x = y
This definition can be extended to evaluate the image of a domain d ∈ D K , the result being also a domain of D K . eval(x := e, d) = {eval(x := e, v) | v ∈ d}

Invariant Propagation

An invariant propagator, or simply propagator, is a function p that takes as input a domain in D K , known as precondition, an instruction in I K , and returns a domain in D K called postcondition.

p : D K × I K -→ D K A propagator is expected to return a valid postcondition for a given precondition and instruction. Formally, this can be described by the correctness condition

∀(d, i) ∈ D K × I K , ∀v ∈ d, eval(i, v) ⊂ p(d, i)
that can also be written:

∀(d, i) ∈ D K × I K , eval(i, d) ⊂ p(d, i). (3.2)
The set of (correct) propagators is noted P K .

Propagators are an abstract representation of the argument theorems used to propagate Floyd's and Hoare's invariants en route to prove a program. Unlike these arguments, propagators have to be defined on the whole domain D K × I K , that is, for any precondition and instruction. It is well known that most mathematical theorems make assumptions on their inputs. To make propagators total functions, the universal postcondition can always be returned if the domain is outside the scope of the corresponding argument, for instance when the hypothesis of some Lyapunov theorem is not verified.

Translation Scheme

The general translation scheme that is implemented is the following (see Figure 3.2). Each instruction i in the abstract controller code with real numbers is turned into a corresponding instruction ī on the concrete controller, where real constants and operators are replaced by their limited-precision counterparts. A post-condition d′ on the limited-precision side is computed using the same proof argument p as in the original code, enlarged to include the error term bound v err that may occur in limited-precision arithmetics in ī. The same propagator p is used because Floyd's and Hoare's invariants are still based on real numbers, even if the code is implemented with a limited-precision arithmetic. v err is computed on an intermediate instruction form ĩ, the meaning of which will be explained later. The precondition d might have been previously replaced by d as the postcondition of a preceding instruction. The error term bound v err depends on the operators of the instruction ī and on the information on the program variables given by d. The limited-precision arithmetic error is not necessarily bounded even when the variables are, e.g. consider division. In this case, we are not able to compute a suitable d′ and the proof translation attempt fails.

% d i % d ′ = p(d, i) % d ī % d′ ⊃ p( d, ĩ) ⊕ v err
Also, d′ will be used as a precondition in the next instruction. An effort has to be made so that d′ has a similar shape to d ′ , in order to fit with the proof argument p ′ used in this instruction, if possible. In our case, we will strive to maintain the invariant under the form of an ellipsoid. Thus, d′ is an overapproximation of p( d, ī) ⊕ v err in the domain of p ′ . In consequence, the automated analysis falls short if such a d′ cannot be found.

Translation Steps

In this section, we consider an instruction i = (x := e) in I R , along with a precondition d and a postcondition d ′ in D R . d, i and d ′ are linked through a propagator p in P R :

d ′ = p(d, i).
As outlined in Section 3.2, we study how the instruction i can be translated into an equivalent instruction ī ∈ I F on limited-precision arithmetic, while using information on the precondition d and propagator p of i to compute an interesting postcondition d′ for ī (given a modification d of d). This is a two-step process:

1. First, the program constants c ∈ R are converted into their limited-precision representations c ∈ F while keeping absolute-precision functions, giving an instruction ĩ ∈ I R . This induces a modified application of the propagator p.

2.

Second, the functions f ∈ F R that appear in the code are changed into their limitedprecision counterparts f ∈ F F . This induces the introduction of rounding errors that are added to the invariant obtained at step 1.

The resulting code relies on limited-precision arithmetic only, which is what is being sought. The translation scheme is described in Figure 3. 3.

% d i % d ′ = p(d, i) % d ĩ % d′ = p( d, ĩ) % d ī % d′ ⊃ p( d, ĩ) ⊕ v err

Converting Constants

We construct an expression ẽ ∈ E R , obtained from e by converting its constants from R to F. Formally, ẽ is defined recursively as follows:

ẽ =      x if e = x ∈ X c if e = c ∈ R f (ẽ 1 , . . . , ẽn ) if e = f (e 1 , . . . , e n ) (3.3)
We let ĩ = (x := ẽ).

Converting Functions

The next step is to convert real functions f that appear in the expression into their limitedprecision counterparts f , as defined in Section 3.1.2. The resulting expression ē ∈ E F is obtained with:

ē =      x if e = x ∈ X c if e = c ∈ R f (ē 1 , . . . , ēn ) if e = f (e 1 , . . . , e n ) (3.4)
Compared to the definition of ẽ in Equation (3.3), we just have turned functions symbols f into f . As previously, we also let ī = (x := ē).

Unlike ĩ, instruction ī relies on different functions than i (taken from F F instead of F R ): the underlying invariant propagation argument in propagator p is very likely not to be applicable to it, as the mathematical theorems used in p typically rely on properties of real arithmetic operators. In this case, p( d, ī) would fall back on the universal domain, which of course is valid but prevents us from doing any interesting proof further in the program. We want to obtain a more precise and still valid postcondition invariant d′ .

Functional Arithmetic Errors

To circumvent this issue, we propose when it is possible to enlarge the invariant d′ = p( d, ĩ) found in Section 3.3.1 to take into account arithmetic errors introduced when converting functions from F R to F F . Such errors are called functional arithmetic errors, but we will refer to them as arithmetic errors, or simply errors, later in this document.

Equation (3.5) defines the arithmetic error on an expression

ẽ ∈ E R modulo a valuation v ∈ Val F .
The result is a value in R.

err(ẽ, v) = | eval(ē, v) -eval(ẽ, v)| (3.5)
This definition is extended to define the arithmetic error induced by an instruction ĩ = (x := ẽ) in I R on a valuation v ∈ Val F . The result is a valuation in Val R defined as follows:

err(x := ẽ, v) = y ∈ X -→ err(ẽ, v) if x = y 0 if x = y .
Finally, a valuation v err ∈ Val F is an upper bound of the arithmetic error induced by instruction

ĩ on a domain d ∈ D R if: ∀v ∈ d, err(ĩ, v) ≤ v err .
∀ d′ ∈ D R , d′ ⊕ v err ⊂ d′ =⇒ eval(ī, d) ⊂ d′ holds.
In other words, if d′ is a postcondition for an instruction ĩ on real numbers, on a given precondition domain d, and d′ is a superset of d′ loose enough to support the arithmetic error v err , then d′ is a valid postcondition for ī, the limited-precision counterpart of instruction ĩ.

A proof in Coq, limited to binary functions (which include all the functions that we need in this document), has been implemented and is available at http://www.

* * *

We would like to check whether the stability proof in Feron's program still holds on a controller implemented with floating-point numbers using the approach of the previous chapter. In this chapter, we use the IEEE Standard for Floating-Point Arithmetic [IEE08] encoded on 64 bits †2 , as most of today's floating-point units do. In this standard, both addition and multiplication are correctly rounded depending on the active rounding mode, which allows to bound the rounding error depending on the values of operands. The case of alternative numeric representations is discussed in Section 4.5.

Automatic Translation

We have developed a program to automatically perform these computations, following the translation scheme framework described in the previous chapter. It is a module called LyaFloat, written in Python and built upon Python libraries SymPy (to handle symbolic mathematics) and Mpmath (for arbitrary-precision floating-point arithmetic).

Here is the listing of a script that automatically computes the floating-point output invariant ellipsoid E R (in variable ERbar) corresponding to E R in the original proof, and tests whether E R ⊂ E P . Two cases are possible:

• either E R ⊂ E P , then the program is Lyapunov-stable on a floating-point architecture;

• or E R ⊂ E P : as E R was obtained through overapproximations, we cannot conclude about the program behavior. row_insert(2, zeros(1, 3)) In this open-loop case, our program LyaFloat is able to check the inclusion. Thus, the stability of the open-loop system with a 64-bit IEEE 754 compliant implementation is formally proven to hold using our proof translation scheme.

Qmu[2,2] = 1 -mu EQmu = Ellipsoid(Qmu) # Symbols xc1, xc2, yc = symbols("xc1 xc2 yc") Xc = Matrix([[xc1], [xc2]]) Yc = Matrix([[yc]]) Zc = Matrix([[xc1], [xc2], [yc]]) # Constant matrices Ac = Matrix(
Notice that the precision of the floating-point arithmetic can be set in setfloatify. This allows us to check that the controller is still stable on a 32-bit only architecture, or to compute that the minimum required precision is 17 bits. This can be useful if, instead of adapting the controller code and/or its proof to the hardware, we have to select a controller device that suits the proof.

In the remainder of this section, we describe the operations performed by our tool to compute E R . It follows the two-step scheme described in Section 3. 3. Alternatively to our program, an implementation of this proof using Mathematica is available online: see Footnote †1, page 27.

Handling Constants

Following Section 3. Apart from constants, the proof scheme for the first part of the program is unchanged:

5 xc = zeros(2,1); % x c ∈ E P 6 receive(y, 2); receive(yd, 3); % x c ∈ E P 7 while (1) % x c ∈ E P 8 yc = max(min(y -yd, 1), -1); % x c ∈ E P , y 2 c ≤ 1 9 skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.

Handling Functions

The second step of our proof translation scheme is to convert function symbols from real to floating-point. We show how it is done instruction by instruction.

Invariant on u

The next instruction after Section 4.2 in the original proof scheme is: We can notice that entering this instruction, the values of matrices C c , D c and E Qµ are known, and the values of x c and y c are bounded by the precondition

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 Listing 4.
( xc yc ) ∈ E Qµ , that is 0.000673593x c 2 1 + 0.000085523x c1 x c2 + 0.00246288x c 2 2 + 0.9991y 2 c ≤ 1. (4.1)
From Equation (4.1), we deduce:

         |x c1 | ≤ 3 • 10
< 20.1614 |y c | ≤ 100 √ 9991 < 1.00046 (4.2)
Here we are able to find algebraic solutions, but this may be impossible with ellipsoids of higher dimension. Still, we would be able to find bounds using numerical methods.

In floating-point arithmetic, rounding errors created by addition and multiplication operators can be bounded when the operands are known or bounded by Equation (4.2), provided that overflow, underflow, and denormalized numbers do not occur [START_REF] Goldberg | What Every Computer Scientist Should Know About Floating-Point Arithmetic[END_REF][START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF].

Here, we need to compute

C c x c + D c y c = C c(0,0) x c1 + D c y c
where all values in the right-hand term are known or bounded. Thus, a constant ε can be computed that bounds the absolute rounding error created when computing u. We obtain:

ε = 5.90 • 10 -12
This way, starting from the algebraic result obtained on real numbers

|u| ≤ (C c D c ) • Q -1 µ • (C c D c ) ⊺
we can ensure that with floating-point numbers, the following inequality holds:

|u| ≤ Ū = (C c D c ) • Q -1 µ • (C c D c ) ⊺ + ε
which leads to the invariants:

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ Ū 2 Listing 4.6: Computation of u, modified

Invariant on x c

The next instruction, considering changes of constants, is:

11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) ⊺ -1 , Listing 4.7: Computation of x c , original code
where R is defined the same way R is, using floating-point terms A c , B c instead of the realvalued counterparts A c , B c , and E R is the ellipsoid built upon R. Again, this invariant holds independently of matrix values.

Here, we compute the values assigned to x c = xc 1 xc 2 :

x c 1 := A c(0,0) x c1 + A c(0,1) x c2 + y c x c 2 := A c(1,0) x c1 + x c2
Using the same method as above, absolute rounding errors introduced by floating-point operations can be bounded on dimensions x c1 and x c2 by constants

ε 1 = 7.42 • 10 -15 and ε 2 = 3.62 • 10 -15 .
These constants must be taken into account in the postcondition. Then the postcondition can be replaced by

% x c ∈ E R
where E R is an ellipse that includes E R plus the rounding error terms (see Figure 4.1). As mentioned in Section 3.2, replacing the ellipse in the postcondition by another ellipse has the advantage of introducing little change in the stability proof sketch (instead of using a different domain, which would involve different theorems), which can greatly facilitate tweaking the rest of the proof in longer codes. Formally, E R must satisfy:

∀x c ∈ E R , ∀x c ′ ∈ R 2 , |x c ′ 1 -x c1 | ≤ ε 1 ∧ |x c ′ 2 -x c2 | ≤ ε 2 =⇒ x c ′ ∈ E R (4.3) E R E R x c1 x c2 ε 1 ε 2 Figure 4.1: Relation between E R and E R
At the end of the proof scheme, the system is stable with floating-point numbers if and only if the inclusion E R ⊂ E P holds. To succeed, E R should be as narrow as possible with respect to Equation (4.3). This is not a clear criterion, as several shapes are possible for E R with no clear winner. Indeed, different criteria can be taken into consideration: area, bounding box, shape, or a combination of criteria… We propose to define E R as the smallest homothety of E R centered around 0 that satisfies Equation (4.3). It can be computed rather easily, for any number of dimension; we give details for two dimensions.

Let a, b, c be the coefficients of E R :

E R = {(x c1 , x c2 ) | ax c 2 1 + bx c 2 2 + cx c1 x c2 ≤ 1}.
a, b and c are known, positive values. Then there exists k ≥ 0 such that

E R = {(x c1 , x c2 ) | ax c 2 1 + bx c 2 2 + cx c1 x c2 ≤ k}. As E R is wider than E R , k ≥ 1.
We need a condition on k that guarantees Equation (4.3).

We consider a point (x c1 , x c2 ) located on the border of E R :

ax c 2 1 + bx c 2 2 + cx c1 x c2 = 1. (4.4)
By construction, for any values e 1 , e 2 ∈ R such that

|e 1 | ≤ ε 1 ∧ |e 2 | ≤ ε 2 the relation (x c1 + e 1 , x c2 + e 2 ) ∈ E R must hold, that is to say: a(x c1 + e 1 ) 2 + b(x c2 + e 2 ) 2 + c(x c1 + e 1 )(x c2 + e 2 ) ≤ k. It develops into (ax c 2 1 + bx c 2 2 + cx c1 x c2 ) + (2ae 1 + ce 2 )x c1 + (2be 2 + ce 1 )x c2 + (ae 2 1 + be 2 2 + ce 1 e 2 ) ≤ k, that is 1 + (2ae 1 + ce 2 )x c1 + (2be 2 + ce 1 )x c2 + (ae 2 1 + be 2 2 + ce 1 e 2 ) ≤ k due to Equation (4.4).
Greatest values for the left-hand term are reached with

|e 1 | = ε 1 ∧ |e 2 | = ε 2 ,
depending on the signs of x c1 and x c2 . As the ellipse E R is symmetric around the origin point (0, 0), we can set e 1 = ε 1 , which leaves only two cases to study. Finally, we numerically verify that greatest values of the term are reached when e 2 = ε 2 . This is the only case we detail here.

We can write:

1 + αx c1 + βx c2 + γ ≤ k with values α = (2aε 1 + cε 2 ), β = (2bε 2 + cε 1 ) and γ = (aε 2 1 + bε 2 2 + cε 1 ε 2 ).
We know that x c1 and x c2 are bounded, thus so is the term αx c1 + βx c2 : we can compute a minimum bound δ such that αx c1 + βx c2 ≤ δ. So it is sufficient that k satisfies:

k ≥ 1 + γ + δ.
Consequently, the smallest homothety of E R that satisfies Equation (4.3) is obtained with k = 1 + γ + δ; we take it as our definition of E R . The instruction becomes:

11 xc = Ac*xc + Bc*yc; % x c ∈ E R Listing 4.8: Computation of x c , modified
In our case, starting from the ellipse that gives E R .

E R = {(x c1 , x c2 ) | 0.00269007x c 2 1 + 0.000341414x c1 x c2 + 0.

End of Proof Scheme

Then, what remains of the stability proof scheme is only to handle inputs and outputs. No change is required and this becomes:

% x c ∈ E R 12 send(u, 1); % x c ∈ E R 13 receive(y, 2); % x c ∈ E R 14 receive(yd, 3); % x c ∈ E R 15 skip; % x c ∈ E P 16 end
Listing 4.9: End of proof scheme

The only step to check carefully is at line 15. At this stage, the final assertion E R ⊂ E P holds and is checked successfully.

Closed-Loop Stability Proof Scheme with Floating-Point Numbers

We now show how the proof of state boundedness of the closed-loop system specifications can be migrated to the level of the controller code and executable model of the system. To be more precise, we exploit the invariance of the ellipsoid E P to develop a proof of proper behavior, that is, stability and variable boundedness, for the computer program that implements the controller as it interacts with the physical system. Unlike the developments related to the open-loop controller, this proof necessarily involves the presence of a model of the physical system. In [START_REF] Feron | From Control Systems to Control Software[END_REF], Feron chooses to represent the physical system and the computer program by two concurrent programs. The code for controller dynamics is unchanged, same as in Section 2. 2.

The pseudocode to represent the physical system dynamics is shown below. In this scheme, the computer program representation of the physical system is to remain unchanged, since it only exists for modeling purposes and does not correspond to any actual program, whereas the controller code is allowed to evolve to reflect the various stages of its implementation.

Establishing proofs of stability of the closed-loop system at the code level is necessarily tied to understanding the joint behavior of the controller and the plant. The entire state space therefore consists of the direct sum of the state spaces of the controller and the physical system. The approach described in the previous sections is used to document the corresponding system of two processes. One interesting aspect of these processes is their concurrency, which can complicate the structure of the state transitions. However, a close inspection of the programs reveals that the transition structure of the processes does not need to rely on the extensions of Hoare's logic to concurrent programs: one program at a time is running, through the blocking nature of the receive primitive.

Feron's stability proof with real numbers is much longer than for the open-loop system. We do not detail it, the interested reader is referred to [START_REF] Feron | From Control Systems to Control Software[END_REF] for full information. To be noticed, the resulting invariants are not much more complex than those available from the study of the controller alone in the open-loop case. On the good side, as already mentioned, the Hoare formalism is not significantly affected by the concurrent structure of the closed-loop system.

A floating-point representation of the closed-loop system consists of keeping the listing corresponding to the physical system in its real-number settings, while replacing the controller part with the corresponding floating-point implementation, as we did in the first part of this chapter.

Using similar techniques to the study of the controller alone, proof invariants can be tweaked to take into account constant changes and rounding errors resulting from the use of floating-point arithmetic in the controller and real arithmetic in the plant.

Unfortunately, these invariants are not sufficient to show that the stability condition holds at the end of the loop body in the case of 64-bit floating-point numbers. In this case, our tool cannot prove the system stability on a double-precision floating-point architecture: either the system is not stable with the floating-point based controller, and in this case the proof parameters (E P , µ, etc.) must be chosen more carefully by the controller designer, or stability holds but we lost it by overapproximating the errors. Tuning the precision of floating-point arithmetic, as described in Section 4.1, teaches us that the controller stability holds on a quadruple-precision platform (i.e., with 128 bits). Thus, an alternative solution might be to use a microcontroller device with superior precision.

Alternative Limited-Precision Arithmetics

Our general idea is to replace some of the invariants in the original proof scheme by wider ones that include rounding errors, with the hope that the stability condition is strong enough and still holds. This approach is made possible by the fact that the rounding errors introduced by the operations used in the code are bounded on bounded inputs and bounded state variables.

In previous sections, we mostly dealt with the case of a floating-point representation of real numbers. They are not available on all architectures, especially on microcontrollers that are commonly used to implement control systems. In this section, we quickly discuss alternative real-number representations.

• We can deal with fixed-point arithmetic the same way we do with floating-point, as long as we stand far enough from extremal values that can lead to overflows or extremely large error terms;

• Another way to represent real numbers is to use two integers, a numerator and a denominator. Considering that the input values are exact, the elementary operations do not introduce rounding errors but can easily lead to overflows, e.g. when computing

p 1 q 1 + p 2 q 2 = p 1 q 2 + p 2 q 1 q 1 q 2 .
A strategy must be used to prevent overflows by introducing approximations: in this case, the question is to quantify the errors introduced by these approximations.

In our example, we exclusively used additions and multiplications: divisions are not involved in linear control. Still, programs with divisions can also be analyzed, if the numerator can be shown to be far enough from zero: it is a supplementary constraint, but it is reasonable to assume that it should be respected on a realistic control system that uses divisions. Differentiable, periodic functions such as (sin) can be computed with a table of values (abacus) and an interpolation function, thus with bounded error. In the same way, functions not periodic, but restricted to finite domains, can also be approximated. Other functions, such as tangent or square root [START_REF] Neron | A Quest for Exactness: Program Transformation for Reliable Real Numbers[END_REF], could raise more issues.

Part II

Linear Relation Analysis

* * *

In this chapter, we give the main definitions and notations related to programs, their semantics, and the properties we want to check about them. First, we introduce transition systems (Section 5.1) that will constitute the standard semantic model of programs, then we define the properties of theses systems, and the general problematic of verification (Section 5.2). In Section 5.3, we introduce interpreted automata, a very general program model.

Transition Systems

Transition systems are the usual semantic models of programs. We give here definitions and notations related to transition systems.

Definition 5.1 (Transition System

). A transition system is a triple (Q, →, Q init ) where Q is a set of states, Q init ⊆ Q is a subset of states called initial states, and → ⊆ Q × Q is a binary relation over Q called transition relation.
Considering two states q, q ′ ∈ Q, the relation (q, q ′ ) ∈ → is usually noted: q → q ′ .

Definition 5.2 (Successor State, Predecessor State).

A state q ′ is a successor of a state q if q → q ′ . Likewise, q is a predecessor of q ′ .

Definition 5.3 (Path, Trace).

A path from a state q to a state q ′ is a sequence (q 0 , . . . , q n ) such as q 0 = q, q n = q ′ and: ∀i < n, q i → q i+1 . A trace is a path whose origin is a state of Q init .

Definition 5.4 (Accessible State, Coaccessible State).

A state q ′ is accessible from a state q (likewise, q is coaccessible from q ′ ) if there exists a path between q and q ′ , that is if (q, q ′ ) belong to the transitive closure of →.

Verification

We consider the framework of property verification (as opposed, for example, to program comparison [START_REF] Boudol | Process Calculi, from Theory to Practice: Verification Tools[END_REF][START_REF] Fernandez | An Implementation of an Efficient Algorithm for Bisimulation Equivalence[END_REF]). More specifically, we consider a linear semantics (as opposed, for example, to arborescent semantics [START_REF] Lamport | Sometime is Sometimes Not Never : On the Temporal Logic of Programs[END_REF]): a behavior of a transition system is a trace, a property is a set of traces. We are essentially interested in safety properties.

Properties

As said, we consider that a property of a transition system is a set of traces. In this paragraph, we assume that these traces are infinite. Let Q ω be the set of infinite sequences of states in Q.

Two important classes of properties are usually distinguished [Lam77, AS87], mainly because different techniques are used to verify them. Informally:

• A safety property expresses that something bad will never happen during the execution.

• A liveness property expresses that something good will necessarily happen at some point as the system is executing.

The violation of a safety property can be detected by a finite execution. Indeed, if a trace σ = (q 0 , q 1 , . . . ) violates the property, then there exists a state q i in σ where something bad occurs and this is definitive since the property supposes that it never happens. So, any trace having (q 0 , q 1 , .

σ: σ ∈ P ⇐⇒ ∀i ≥ 0, ∃σ ′ ∈ Q ω such as σ[. . . i] • σ ′ ∈ P (5.1)
where σ[. . . i] is the finite sequence of the i first terms of σ, and • is the sequence concatenation operator.

Usual examples of safety properties include the sequential behavior of transition systems, mutual exclusion in products of transition systems, the absence of division by zero, the validity of array accesses or pointer dereference, or the absence of deadlocks.

If P is a liveness property, no finite execution is definitively bad, since something good can always happen in the future. Every finite execution can be prolonged in a trace that satisfies P : Definition 5.6 (Liveness Property). P ⊆ Q ω is a liveness property if and only if

∀σ ∈ Q * , ∃σ ′ ∈ Q ω such as σ • σ ′ ∈ P. (5.2)
A usual example of liveness property is program termination.

Example. We consider a program with two threads. Each thread can be either waiting to be executed (0), working (1) or terminated (2). The state space of the program is

Q = {0, 1, 2} × {0, 1, 2} and the unique initial state is Q init = {(0, 0)}.
The mutual exclusion property P mutex = {σ ∈ Q ω | ∀q ∈ σ, q = (1, 1)} ensures that the two threads cannot be working at the same time. It is a safety property, as defined by Equation (5.1): indeed, given a trace σ ∈ P , each finite prefix σ[. . . i] can be extended, for instance with the infinite prefix ((0, 0), (0, 0), . . . ), into a trace that verify P mutex .

The termination property In this dissertation, we are only interested in verifying safety properties.

P join = {σ ∈ Q ω | ∃q ∈ σ, q = (2,

Accessibility and Verification

As mentioned above, a safety property expresses that something bad never happens . Accessibility properties constitute a subset of safety properties, expressing that a certain set of states -representing something bad , and noted Q err -cannot be reached. Their verification involves the computation of accessible states (from the initial states Q init ), or of coaccessible states (from the error states Q err ).

We introduce a few usual definitions relative to accessibility. Let X be a subset of Q.

Definition 5.7 (Postconditions, Preconditions). We note:

Post(X) = {q ′ ∈ Q | ∃q ∈ X, q → q ′ } and Pre(X) = {q ∈ Q | ∃q ′ ∈ X, q → q ′ }.
Post(X) is the set of successor states of all states in X. Similarly, Pre(X) is the set of all predecessors of states in X.

Definition 5.8 (Accessible States, Coaccessible States).

We note acc(X) (respectively coacc(X)) the set of all accessible states (respectively coaccessible states) from some state in X. As acc(X) (respectively coacc(X)) is the smallest solution of the fixpoint equation

Y = X ∪ Post(Y ) (respectively Y = X ∪ Pre(Y )), we have: acc(X) = n≥0 Post n (X) and coacc(X) = n≥0
Pre n (X).

It appears that there are two ways to verify an accessibility property:

• Compute acc(Q init ) and verify that acc(Q init ) ∩ Q err = ∅.

• Compute coacc(Q err ) and verify that

Q init ∩ coacc(Q err ) = ∅.
Therefore, the crux of verifying accessibility properties is fixpoint computation. In some cases, this compution can be done exactly: in these cases, we speak of model checking [QS82, CES86] (Section 5.2.3). But in the general case, only approximations of fixpoints can be computed automatically: then, we speak of conservative verification (Section 5.2.4).

Model Checking

Model checking mostly concerns finite-state systems. In this case, iterative fixpoint computations always converge and all properties are, theoretically, decidable. Efforts on model checking are about, on the one hand, the verification of very general properties -both of safety and of liveness, and with linear as well as arborescent semantics (temporal logic) -and, on the other hand, the optimization techniques allowing to work on huge sets of states. As the verification is exact, its conclusion consists either in a satisfaction verdict, or in an error verdict that can come with a counterexample.

Moreover, model checking is also used in the case of infinite state systems, when they can be abstracted (ignoring some aspects) into finite systems [START_REF] Graf | A Tool for Symbolic Program Verification and Abstraction[END_REF][START_REF] Clarke | Model Checking and Abstraction[END_REF]. In this case, we talk of abstract model checking. It should be noted that in case of failure, this kind of verification cannot conclude that the property is not satisfiable, as the failure can be due to the approximations that were made during abstractions.

Conservative Verification

If every type of property requires fixpoint computations, approximate computations can be enough in the case of accessibility properties: if we compute a superset Q ′ of acc(Q init ) (respectively, of coacc(Q init )) and if it does not intersect Q err (respectively, Q init ), then the property is certainly satisfied (Figure 5.1). However, if the intersection is not empty, we do not know whether the property is violated or if the overapproximation used is too inaccurate. Conversely, if we compute a subset of acc(Q init ) (respectively, of coacc(Q err )) that intersects Q err (respectively, Q init ), the property is certainly violated. In both cases, we speak of conservative verification: in case of success, the verification provides a guaranteed result, but in case of failure, it is inconclusive.

acc(Q init ) Q ′ Q err Figure 5.1: Conservative Verification
The approximation of fixpoints is the central theme of the theory of abstract interpretation that is presented in Chapter 6.

Program Models

Transition systems constitute a very low-level model. We need a higher-level model of programs, but, as we do not want to rely on a particular programming language, we choose a relatively abstract and general model: interpreted automata. We give in Section 5. 3.2 a semantics of interpreted automata in terms of transition systems. Also, the translation of a safety property into an accessibility property is defined at the interpreted automaton level (Section 5.3.5).

Interpreted Automata

An interpreted automaton [bor95] is composed of a finite set of control points K, including one initial point. Each control point is mapped to a set of guarded commands, that act on variables and lead to other control points.

Let X be a finite set of typed variables (for example, integer and boolean variables). A valuation is a function mapping each variable to one value of its type. Let Val be the set of valuations on X.

Definition 5.9 (Guard, Action, Guarded Command).

• A guard g is a function g : Val → B = {⊤, ⊥}, giving the value of a logic formula on a valuation V . In practice, the guard is assimilated to its truth set (the set it is a characteristic function of): {V | g(V ) = ⊤}. The set of guards is noted Guards.

• An action a is a binary relation on valuations (a ⊆ Val × Val). The set of actions is noted Act.

• A command is a quadruplet c = (k, g, a, k ′ ) ∈ K × Guards × Act ×K. Intuitively, when the automaton is at control point k, and if the current valuation V satisfies the guard g (i.e. g(V ) = ⊤), the command may be executed, and when executed it moves to the control point k ′ , with a valuation

V ′ such that (V, V ′ ) ∈ a.
The actions are by default nondeterministic. In particular, a read action (noted x := ? ) maps a valuation V to any valuation V ′ that differs from V only by the value of x (in other words,

x can take any value). Thus, the value of some variables (input variables) can be completely unspecified at some control points. At transition level, the choice of the next state can also be nondeterministic.

Definition 5.10 (Interpeted Automaton

). An interpreted automaton is a triple α = (K, Com, k init ) where K is a finite set of control points, Com is a finite set of commands, and k init ∈ K is the initial control point.

Interpreted automata are a refinement of finite state machines: they have state and transitions, but also variables that are involved in transitions. Also, the terminology is different: transitions in finite state machines are called commands and states are called control points.

An example of interpreted automaton is given in Figure 5. 2. It corresponds to the interpreted automaton

(K, Com, k init ) with K = {k init , k 1 }, Com = {c 1 , c 2 , c 3 } and c 1 = k init , (λx, y : cos y ≤ x), {(x, y : x + 1, y) | x, y ∈ Val }, k 1 c 2 = k 1 , (λx, y : ⊤), {(x, y, x, y + 2) | x, y ∈ Val }, k 1 c 3 = k 1 , (λx, y : x = y), {(x, y, x, y x ) | x, y ∈ Val }, k init k init k 1 c 1 : cos y ≤ x? → x := x + 1

Semantics in Terms of Transition Systems

We give the semantics of interpreted automata in terms of transition systems. Let (K, Com, k init ) be an interpreted automaton. It is associated to a transition system (Q, →, Q init ) defined as follows:

• The set of states is the set of pairs (k, V ) where k is a control point and V is a valuation:

Q = K × Val.
• Initial states have k init as control point, and any valuation:

Q init = {k init } × Val.
• The transition relation is derived from the command semantics:

(k, V ) → (k ′ , V ′ ) ⇐⇒ ∃(k, g, a, k ′ ) ∈ Com, g(V ) = ⊤ ∧ (V, V ′ ) ∈ a.

Example

We illustrate the modeling of a program by an interpreted automaton on a very small example, from [START_REF] Halbwachs | Delay Analysis in Synchronous Programs[END_REF]. We want to model a car 1. Whose speed is always lower or equal to 2 m/s. 2. That stops after 4 seconds. 3. And should not collide into a wall after 10 meters when leaving the track.

Our purpose is to show that, under assumptions (1) and ( 2), the car stops before hitting the wall.

The simulation program has two boolean inputs m and s (respectively true when a meter is covered and a second is elapsed, the events are supposed to be not simultaneous), read at every transition. The program counts the distance traveled (d, a counter of meters), the elapsed time (t, a counter of seconds) and the instant speed (v, a counter of meters by second, reset at each second), stops the car if it goes too fast, or hits the wall after 10 m are covered.

The corresponding automaton is represented in Figure 5. 3. As in [START_REF] Merchat | Réduction du nombre de variables en analyse de relations linéaires[END_REF], the automaton can be slightly modified to only take into account numerical aspects (Figure 5.4). We get rid of the inputs, and we can see that state k fast is not accessible. It turns out that, as the system is finite-state, the corresponding transition system is also finite. On this transition system, it is easy to verify that the state k wall is unreachable. Remarks.

k 0 k 1 k fast k wall k stop d := 0; t := 0; v := 0 t ≥ 4 d ≥ 10 v ≥ 3 s? → v := 0; t := t + 1 m? → v := v + 1; d := d + 1
k 0 k 1 k fast k wall k stop d := 0; t := 0; v := 0 t ≥ 3 d ≥ 10 v ≥ 3 t ≤ 2? → v := 0; t := t + 1 v ≤ 1, d ≤ 8? → v := v + 1; d := d + 1
• As this example is finite-state, model checking can be used.

• If the constants appearing in the hypotheses (2 m/s, 4 s, 10 m) are increased, the number of states in the transition system tends to explode, making model checking difficult. If these constants are parameters, given symbolically, then model checking cannot be used.

Collecting Semantics

To verify an accessibility property in an interpreted automaton, it is necessary to compute, or approximate, either the set acc({k init } × Val) of accessible states in the associated transition system, either the set coacc({k err } × Val) of coaccessible states, assuming that error states are characterized by a control point k err (in the car example, the control point k wall ).

To this end, instead of directly applying the semantics defined in Section 5. 3.2, we can try to take advantage of the control structure of the automaton to break up the problem.

For every k ∈ K, the set of accessible valuations (respectively, coaccessible valuations) at point

k is noted A k (respectively, C k ): A k = {V | (k, V ) ∈ acc({k init } × Val)} and C k = {V | (k, V ) ∈ coacc({k err } × Val)}.
Of course, we have:

acc({k init } × Val) = k∈K A k and coacc({k err } × Val) = k∈K C k .
Then, the sets (A k ) k∈K (respectively, (C k ) k∈K ) can be caracterized as the smallest solution of a forward (respectively, backward) system of semantics equations, defined as follows:

• Forward system:

A k init = Val and ∀k = k init , A k = (k ′ ,g,a,k)∈Com a(A k ′ ∩ g) (5.3)
where a(X) = {V ′ | ∃V ∈ X, (V, V ′ ) ∈ a} (image of X under the relation a).

• Backward system:

C kerr = Val and ∀k = k err , C k = (k,g,a,k ′ )∈Com a -1 (C k ′ ) ∩ g (5.4)
where a -1 is the inverse relation of a, and

a -1 = {(V ′ , V ) | (V, V ′ ) ∈ a}.
For instance, in Figure 5.5, the equation associated to point k is the following: Intuitively, in any non-initial point k, a valuation V accessible in k is obtained through the command (k ′ , g, a, k) on a valuation V ′ accessible in k ′ , satisfying g, and connected to V by a ((V ′ , V ) ∈ a). Similarly, in any point k different from the target k err , a valuation V coaccessible in k must lead, through the command (k, g, a, k ′ ), to a valuation V ′ coaccessible in k ′ , that is to say that V must satisfy g, and that (V, V ′ ) must belong to a.

A k = a 1 (A k 1 ∩ g 1 ) ∪ a 2 (A k 2 ∩ g 2 ) ∪ a 3 (A k ∩ g 3 ). k 1 k k 2 g 1 ? → a 1 g 2 ? → a 2 g 3 ? → a 3
The advantage of this semantics (often called collecting semantics ) is that it gives an equation system of fixed points, instead of a single equation. We will show in Chapter 6 that this partitioning allows to greatly increase the accuracy of computations in the case of an approximated computation of the solution.

Composition of Interpreted Automata and Observers

The purpose of this section is to show how safety properties can be turned into accessibility properties, thanks to the notion of synchronous observer [START_REF] Halbwachs | Synchronous Observers and the Verification of Reactive Systems[END_REF]. The idea is to associate to a safety property an interpreted automaton called observer, that works as a recognizer for traces satisfying the property. Appropriatly composed with the automaton to verify (that is to say, in synchronous parallelism), the observer detects every violation of the property by joining an error control point. This way, the verification consists into making sure that the composition of the initial automaton with its observer cannot reach a state where the observer is in its error point, which is an accessibility property.

First, we define the parallel synchronous composition on interpreted automata: intuitively, two automata composed with parallel synchronisation will execute commands together, provided that the values of variables coincide on both sides -the idea being that a variable is either local to an automaton (and so, left unspecified in the other), or computed in one automaton and read (and so, left unspecified) in the other.

Definition 5.11 (Synchronous product). The synchronous product

A × A ′ of two inter- preted automata A = (K, Com, k init ) and A ′ = (K ′ , Com ′ , k ′ init ) with a common set of variables X is the automaton (K × K ′ , Com × , (k init , k ′ init ))
where

Com × = {((k 1 , k ′ 1 ), g ∧ g ′ , a ∩ a ′ , (k 2 , k ′ 2 )) | (k 1 , g, a, k 2 ) ∈ Com, (k ′ 1 , g ′ , a ′ , k ′ 2 ) ∈ Com ′ }.
In other words, a command of the product from a composed state (k 1 , k ′ 1 ) consists in the (synchronous) composition of a command for each component: the synchronous composition of two commands (g, a) and (g ′ , a ′ ) on a valuation V is only possible if both guards are satisfied ((g ∧ g ′ )(V )), and its result V ′ is a common result of both actions ((V, V ′ ) ∈ (a ∩ a ′ )).

Let A be an interpreted automaton, working on a set of variables X. Let P be a safety property. Definition 5.12 (Observer). An observer of P is an interpreted automaton O working on a partitioned set of variables X ⊎ X ′ and such that • Variables in X are left unspecified in O.

• O has a distinctive control point, noted k err , that the automaton reaches (and never leaves) if and only if the evolution of variables in X violates P . Formally, consider

((k 0 , V 0 ⊎ V ′ 0 ), (k 1 , V 1 ⊎ V ′ 1 ), . . . , (k n , V n ⊎ V ′ n )) a trace in O, then k n = k err ⇐⇒ ((k 0 , V 0 ), (k 1 , V 1 ), . . . , (k n , V n )) ∈ P.
K err = {(k, k err ) | k ∈ K}.
Thus, every safety property for which an observer can be provided, can be turned into an accessibility property.

Chapter 6

Abstract Interpretation

L'interprétation abstraite est une théorie d'approximation de la sémantique des programmes informatiques fondée sur les fonctions monotones pour ensembles ordonnés, en particulier les treillis. Elle peut être définie comme une exécution partielle d'un programme pour obtenir des informations sur sa sémantique (par exemple, sa structure de contrôle, son flot de données) sans avoir à en faire le traitement complet.

La principale fonction de l'interprétation abstraite est l'analyse statique, l'extraction automatique d'informations sur les exécutions possibles d'un programme. Ce chapitre définit les notions de base de l'interprétation abstraite, qui sont utilisées dans la suite du document.

Introduction

In the previous chapter, we have seen the kind of systems and properties we want to verify. In particular, we have seen

• How a safety property can be translated into an accessibility property (Section 5.3.5).

• How verifying an accessibility property can boil down to computing a set of accessible states acc(Q init ) or coaccessible states coacc(Q err ) (Section 5.2.2).

• Finally, than these sets of states were solutions of fixed-point equations (Section 5.2.2)

acc(Q init ) = Q init ∪ Post(acc(Q init )) coacc(Q err ) = Q err ∪ Pre(coacc(Q err ))
or of systems of fixed-point equations (Section 5. 3.4)

A k init = Val and ∀k = k init , A k = (k ′ ,g,a,k)∈Com a(A k ′ ∩ g) C kerr = Val and ∀k = k err , C k = (k,g,a,k ′ )∈Com a -1 (C k ′ ) ∩ g.
Usual theorems ensure that the considered fixed-point equations have least solutions, which are well-characterized: according to Kleene's fixed-point theorem, an equation X = F (X) -F being a continuous function on a complete lattice -has a least solution

lfp(F ) = n≥0 F n (⊥) (6.1)
where ⊔ is the least upper bound operator, and ⊥ is the least element in the complete lattice.

In the present case, the complete lattice is the set of parts of Q ordered by inclusion, the upper bound operator is the union of sets of states, and the least element is the empty set.

In general, solving Equation (6.1) is impossible. When the set of system states Q is infinite, computing the ascending Kleene chain F (⊥), F (F (⊥)), F (F (F (⊥))), . . . raises two kinds of issues:

• It requires to be able to represent arbitrary sets of states and to compute on them.

• The iterative computation may not terminate.

However, it was pointed out in Section 5.2.4 that accessibility properties could be proved in a conservative way, by computing an overapproximation of the solution. The approximate resolution of fixed-point equations is the main objective of abstract interpretation, a domain that was formalized by Patrick and Radhia Cousot in the late 1970s [START_REF] Cousot | Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approximation of Fixpoints[END_REF], whose main principles will be recalled in the first part of this chapter (Section 6.2). The second part provides an overview of abstract interpretations dedicated to numerical properties (Sections 6.3, 6.4 and 6.5).

Principles of Abstract Interpretation

We have identified two obstacles to the exact resolution of fixed-point equations, that are:

• The impossibility to represent complex values -in general, infinite sets of statesand to compute on those values.

• The impossibility, in the general case, to compute the limit of infinite iterations.

Abstract interpretation proposes approximate solutions to these two problems. On the one hand, a notion of abstract domain is proposed, that allows to compute on abstract values simpler than sets of states. On the other hand, depending on the choice of these abstract values, the finiteness of iterations can be guaranteed; if not, techniques are provided to extrapolate and overapproximate the limit.

Abstract Domains

An abstract interpretation is first characterized by a domain of abstract values. These values are used to approximate, or abstract, concrete values (sets of states) that are too complex. The notion of approximation is formalized by an ordering relation: the abstract domain must be a complete lattice, as the concrete domain. The initial fixed-point equation on concrete values is then associated to an abstract fixed-point equation, to be solved in the abstract domain.

Abstract values must be chosen so that their computer representation is possible, as well as the operations necessary to the iterative computation of the abstract equation.

We have seen that the concrete lattice is (P(Q), ⊆, ∪, ∩, Q, ∅), that is the power set of Q, ordered by inclusion, with ∪ and ∩ as upper and lower bound operations, and Q and ∅ as greatest and least elements. As the theorical framework of abstract interpretation is more general, we note (L, ⊑, ⊔, ⊓, ⊤, ⊥) this concrete lattice. The lattice of abstract values is noted

(L ♯ , ⊑ ♯ , ⊔ ♯ , ⊓ ♯ , ⊤ ♯ , ⊥ ♯ ).
The relation between concrete and abstract values may be handled by the notion of Galois connection: a Galois connection between L and L ♯ is a pair of functions α : L → L ♯ (abstraction) and γ : L ♯ → L (concretization) such that (Figure 6.1): As there exists a Galois connection between concrete and abstract lattices, the following theorem can be applied to the computation of fixed-point overapproximations:

∀x ∈ L, ∀y ∈ L ♯ , α(x) ⊑ ♯ y ⇐⇒ x ⊑ γ(y).
Theorem 6.1 (Fixed-point overapproximation). If L and L ♯ are connected by a Galois connection (α, γ), if F is a continuous function in L [GHK + 03] and G is a continuous function in L ♯ , if ∀x ∈ L, α(F (x)) ⊑ ♯ G(α(x)) (6.2)
then lfp(F ) ⊑ γ(lfp(G)).

This theorem means that, to compute an overapproximation of the least fixed-point of a concrete function F , we can choose or build an abstract function G satisfying Equation (6.2), compute the least fixed-point of G in the abstract lattice, and concretize the result.

Abstract Operations

It can be shown that the best choice for the abstract function G is α • F • γ. However, this function is not computable in the general case, since it requires computations in the concrete domain. On the other hand, we aim to have an automatic method to compute G from F . We take advantage of the fact that the function F , in the case of program analysis, is constructed by composition of elementary functions on sets of states (the preconditions or postconditions of elementary instructions in the programming language), of unions and of intersections. Similarly, an abstract function G can be mechanically constructed from elementary abstract functions, associated to language instructions, and of upper and lower bounds operations (⊔ ♯ and ⊓ ♯ ). Therefore, once the abstract lattice and the associated Galois connection are chosen, designing an abstract interpretation requires to define these abstract functions, and to provide effective algorithms.

Thus, in the case of interpreted automata, and for forward analysis (Equation (5.3), page 55), we must give

• An interpretation of guards in the abstract lattice, that is to say an effective way to associate to each guard g an abstract value g ♯ such that α(g) ⊑ ♯ g ♯ (or: g ⊑ γ(g ♯ )).

• An interpretation of actions, that is to say an effective way to associate to each action a an abstract function a ♯ such that ∀y ∈ L, α(a(γ(y))) ⊑ ♯ a ♯ (y) (or: ∀x ∈ L, a(x) ⊑ γ(a ♯ (x))).

So, the concrete equation ( 5. 3)

A k = (k ′ ,g,a,k)∈Com a(A k ′ ∩ g)
is associated to the abstract equation

A ♯ k = ♯ (k ′ ,g,a,k)∈Com a ♯ (A ♯ k ′ ⊓ ♯ g ♯ ).

Widening and Decreasing Sequence

In many cases, the iterative fixed-point computation in the abstract lattice converges with a finite number of iterations: it is obviously the case when the abstract lattice is finite, but also when simply its depth is finite (i.e., every strictly increasing chain is finite). In these cases, the problem of infinite iterations is solved. However, it may be wise for accuracy purposes [START_REF] Cousot | Abstract Interpretation and Application to Logic Programs[END_REF] to use more expressive lattices that do not satisfy these properties. In this case, the idea is to guess the limit of an infinite sequence using its first terms. To do this, the abstract interpretation comes with a widening operator, used to extrapolate the limit of the iterative sequence.

Definition 6.2 (Widening).

A widening operator is a function ∇ : L ♯ ×L ♯ → L ♯ that satisfies the following properties:

1. ∀y 1 , y 2 ∈ L ♯ , (y 1 ⊔ ♯ y 2 ) ⊑ ♯ (y 1 ∇y 2 ).
2. For any increasing sequence (x n ) n≥0 in L ♯ , the sequence defined by y 0 = x 0 , y n+1 = y n ∇x n+1 converges in a finite number of iterations. Theorem 6.3 (Extrapolation). If ∇ is a widening operator on L ♯ and G is a continuous function on L ♯ , then the sequence y 0 = ⊥, y n+1 = y n ∇G(y n ) converges in a finite number of steps towards a limit ŷ that is a post-fixpoint of G (i.e., ŷ ⊒ ♯ G(ŷ)) and so, such that ŷ ⊒ ♯ lfp(G).

A trivial example of widening operator is: ∀x, y ∈ L ♯ , x∇y = ⊤.

The above theorem allows computing in every cases an approximation of the least fixed point of G. If the limit ŷ is not a fixed point, the solution is improved by the following theorem: Theorem 6.4 (Decreasing Sequence). If ŷ = G(ŷ), then the sequence z 0 = ŷ, z n+1 = G(z n ) is decreasing, and all its terms are overapproximations of lfp(G).

In short, to overapproximate the least fixed point of a function G, instead of computing the exact (and generally infinite) sequence

x 0 = ⊥, x n+1 = G(x n ).
we first compute an extrapolation sequence y 0 = ⊥, y n+1 = y n ∇G(y n ) that converges in a finite number of steps toward a limit ŷ, that is a correct overapproximation of lfp(G). Furthermore, if ŷ = G(ŷ), the result can be improved by computing the first terms (up to any finite given precision) of the decreasing sequence

z 0 = ŷ, z n+1 = G(z n )
that may be finite or infinite, but whose all terms are correct overapproximations of lfp(G).

Abstract Interpretations of Interpreted Automata

We have seen, in Section 5. 3.4, that the analysis of an interpreted automaton could lead to solving a system of fixed point equations, based on the control structure of the automaton, rather than a single equation. They are called partitioned systems. We focus on the forward system, the backward case being similar:

A k init = Val, ∀k = k init , A k = (k ′ ,g,a,k)∈Com a(A k ′ ∩ g)
This system is associated to the following abstract equation system:

A ♯ k init = ⊤ ♯ , ∀k = k init , A ♯ k = ♯ (k ′ ,g,a,k)∈Com a ♯ (A ♯ k ′ ⊓ ♯ g ♯ ).
This system can be solved iteratively, considering the equations one by one in any order, until stabilization. The convergence speed in influenced by the order in which the equations are considered, and the application of widening can be optimized:

• It makes sense to consider equations taking into consideration the control structure: when A ♯ k changes, the values A ♯ k ′ associated to the control points k ′ that are successors of k are likely to change too. Moreover, there is no use to consider control points located downstream a strongly connected component in the control graph, before stabilizing the values associated to control points in the component.

• It is unnecessary to use widening (that leads to information loss) in every equations of the system to ensure convergence: it is sufficient to choose a set K ∇ ⊆ K of widening control points such that every cycle in the graph contains a point of K ∇ , and to apply widening only in equations corresponding to widening points. As the choice of a minimal set of widening points is an NP-complete problem, we simply rely on an heuristic choice [Bou92].

Approximating Numeric Sets

0.5

We are interested by the verification of numerical programs, and we consider specifically the case where the state of the program is defined by a control point and a valuation of variables in a numerical set N (Z, Q or R). Using the semantics defined in Section 5. 3.4, we need to represent subsets of N n where n is the number of variables. Thus, the concrete lattice is:

P(N n ), ⊆, ∪, ∩, N n , ∅ .
In this section, we provide a brief overview of the main abstract lattices that have been proposed in the literature to abstract this domain.

The Sign Lattice

A rough abstraction consists in associating to each variable an abstract value in the lattice represented in Figure 6.2. For instance, the set of points shown in Figure 6.3(a) would be abstracted by the first, strict quadrant {x, y | x > 0, y > 0} (Figure 6.3(b)). In this lattice, abstract operations are straightforward, and, as the lattice is finite, there is no need for widening.

⊤ = 0 ≤ 0 ≥ 0 < 0 = 0 > 0 ⊥

The Lattice of Affine Equations

One of the first nontrivial abstractions, introduced by [START_REF] Karr | Affine Relationships among Variables of a Program[END_REF], consists in synthesizing invariant linear equations. In this analysis, a set of points in N n is abstracted by the smallest affine variety that contains it, that is to say by solution set of an affine system. Using this analysis, the set in Figure 6.3(a) would be abstracted into the whole space (no information, no equation), but the one in Figure 6.4(a) gives the equation x -2y = 1 (Figure 6.4(b)). The lattice of affine equations is infinite, but its depth is finite (there are at most n linearly independent varieties in N n ). Thus, there is no point in defining a widening operator.

The Lattice of Linear Congruences

In his thesis [START_REF] Granger | Analyses sémantiques de congruence[END_REF], Philippe Granger proposes, in particular, a lattice for abstracting parts of Z n as linear congruences, in the form a x ≡ b (mod c) (that is, ∃k ∈ Z, a x = kc+b). For example, the set in Figure 6.5(a) as abstracted as shown in Figure 6.5(b), that is x -2y ≡ 2 (mod 6).

The lattice is infinite, with infinite depth, but the increasing chain condition is satisfied: there is no infinite sequence of strictly increasing abstract values. Once again, widening is useless.

The Lattice of Intervals

The interval lattice [START_REF] Cousot | Static Determination of Dynamic Properties of Generalized Type Unions[END_REF] determines the variation interval of numerical variables. The abstraction of a set E consists to associate to each variable v the pair [min v , max v ] of minimal and maximal values (possibly -∞, +∞) of v in E. The abstraction by intervals of the set represented in Figure 6.6(a) is given in Figure 6.6(b). In the interval lattice, it is necessary to define a widening operator to ensure convergence. In [START_REF] Cousot | Static Determination of Dynamic Properties of Generalized Type Unions[END_REF], the following widening operator is proposed:

[a, b]∇[c, d] = [if c < a then -∞ else a, if d > b then + ∞ else b] x y (a)
x y (b) Figure 6.6: Abstraction Using Intervals

The Lattice of Convex Polyhedra

The lattice of convex polyhedra was first introduced by [CH78, Hal79] to analyse linear inequalities that are invariant amongst the variables of a program. The abstraction of a set of points in N n is the smallest convex polyhedron that contains it (Figure 6.7). Let us notice that while this definition makes sense for a finite set of points, this is not the case of an infinite set whose convex hull is not necessarily polyhedral (for instance, a sphere). This difficulty arises from the fact that the polyhedra lattice is not complete (the upper bound of a infinite set of polyhedra is not necessarily of polyhedron). The actual, complete lattice is the lattice of convex sets in N n , but as the computations are always finite (thanks, in particular, to widening), the computed solutions are always polyhedra. This lattice is used in this dissertation and is discussed in greater details in Chapter 7.

Particular Polyhedra

As the worst-case complexity of operations on polyhedra is exponential, some subclasses have been identified for specific applications.

Octagons

Areas have been generalized, without increasing complexity, by Antoine Miné [START_REF] Miné | The Octagon Abstract Domain[END_REF], that also allows constraints in the form x + y op c. The resulting polyhedra are called octagons (Figure 6.9).

Octahedra

Octagons have geen generalized too into octahedra [START_REF] Clarisó | The Octahedron Abstract Domain[END_REF], considering linear constraints whose coefficients are in {-1, 0, 1} (so, that may involve more than two variables). However, polynomial complexity is not preserved.

Semi-Linear Sets and Presburger Arithmetic

We finish this short review of numerical set representation techniques by citing works less directly related to abstract interpretation, that aim to analyze in an exact way the possible values of integer variables during a program execution: some authors [BW94, BGP97, CJ98, FS00] use the decidability of Presburger arithmetic to represent more accurately infinite sets of integer vectors. As long as the sets of states to be represented are semilinear, algorithms to handle them are available. For some programs, the set of accessible states can be computed with accelerations [START_REF] Gonnord | Accélération abstraite pour l'amélioration de la précision en analyse des relations linéaires[END_REF]. Acceleration is recognizing a fixed common pattern and exhibiting hardcoded solutions for this. It is used instead of widening in relevant cases and gives exact results.

Uses of Numerical Variables Analysis

There are various applications to the analysis of behavior of numerical variables. We classify them in two parts:

• Applications with a goal given a priori, in general a set of states we want to verify the inaccessibility of, and that can be taken into account by the analysis.

• Applications in which we simply want to get some information about, in general invariants, as precise as possible, to use it at a later stage.

Verifying Inaccessibility

We have seen in Section 5. 3.5 that all safety properties can be expressed as the inaccessibility of some states. The same applies to verify the absence of errors at runtime (dynamic consistency), for which there is no need to specify the property, and that often rely on numerical aspects (divide by zero, out-of-bounds array access [START_REF] Nguyen | Efficient Intraprocedural Array Bound Checking[END_REF], arithmetic overflow).

In all these applications, the set Q err of states we want to prove inaccessible is known, and we can proceed either forward -to prove acc(Q init ) ∩ Q err = ∅ -or backward -to prove

Q init ∩ coacc(Q err ) = ∅.
As it is an approximate analysis, we compute in fact either acc(Q init ) or coacc(Q err ), with

acc(Q init ) ⊆ acc(Q init ) and coacc(Q err ) ⊆ coacc(Q err ).
In terms of accuracy, it is beneficial to iterate forward and backward analyses, by computing successively

X 0 = acc(Q init ), Y n = coacc(Q err ∩ X n ), X n+1 = acc(Q init ∩ Y n )
that constitute a sequence increasingly precise, as long as the terms X n or Y n are not empty and do not stabilize.

Invariant Synthesis

A second class of applications is about invariant synthesis, with no preconceived objective, to be used later on. These applications are about:

• Compilation, with, for example, dead code elimination, finding expressions that can be computed at compile time or outside loops, the characterization of variation intervals of variables to optimize their memory implementation, the determination of the precise lifespan of array elements to improve memory management [LMC02]…

• Program parallelization [START_REF] Ancourt | A Linear Algebra Framework for Static HPF Code Distribution[END_REF] where an accurate analysis of inter-variable dependencies is required, especially between array elements.

• Interactive proof, where invariants found by analysis are used to make the human's work easier [BBC + 00, BLO98].

• In some applications, it is necessary to add auxiliary counters to the program, and in the analysis of their behavior. For example, [START_REF] Mercouroff | An Algorithm for Analyzing Communicating Processes[END_REF] proposes a method of communication analysis in a multiprocess system by adding counters at each meeting point, analyzing the code of each process one at a time and then, combining results while equalizing counters of corresponding meetings. We can also try to show loop termination (a liveness property), or try to determine the number of iterations [START_REF] Podelski | A Complete Method for the Synthesis of Linear Ranking Functions[END_REF], by adding a loop counter and then bounding it.

Chapter 7

The Lattice of Polyhedra details the translation from one representation to another, which is the key algorithm to perform all other operations, and most importantly to simplify representations. Section 7.3 describes the usual operations on convex polyhedra. An example of analysis is then detailed in Section 7. 4. Finally, the main existing polyhedral libraries are presented in Section 7. 5.

The general references of this chapter are [START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF], [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF] and [Jea00].

We consider a numerical set N and the affine space N n . The scalar product of two vectors X and Y in N n , equals to n i=1 X i Y i , is noted X • Y or simply XY .

Two Representations Definition (Convex polyhedron).

A closed, convex polyhedron (or simply polyhedron) is the intersection of a finite number m of closed half-spaces in N n .

Linear Constraints

Proposition 7.2 (Representation with affine constraints). Any convex polyhedron can be characterized as the solution set of a linear constraint system AX ≤ B, where A is an m × n matrix and B is a vector in N m . We call linear constraint an inequation aX ≤ b or an equation aX = b, where a ∈ N n is a line vector and b ∈ N is a scalar.

Example. Figure 7.1 represents the polyhedron corresponding to the constraint system In practice this is the representation we use, as it gives direct information on variables.

2x 2 -x 1 ≤ 7 -x 1 -x 2 ≤ -5 -x 2 ≤ -1 in other terms, AX ≤ B with A =   -1 2 -1 -1 0 -1   , X = x 1 x 2 and B =   7 -5 -1   . x 1 x 2 x 1 + 7 ≥ 2x 2 x 2 ≥ 5 -x 1 x 2 ≥ 1

Non-Canonicity of Constraint Systems

In this paragraph, we focus on the non-canonicity of the representation by linear constraints.

On the one hand, a constraint system may contain redundant constraints; on the other, when the polyhedron is contained within a strict vector subspace of N n (i.e., when it satisfies linear equations), it may be characterized by several minimal constraint systems.

Definition 7.3 (Redundant constraints).

A constraint aX ≤ b is said redundant in the system AX ≤ B if it can be removed while keeping unchanged the described polyhedron. In other words, the satisfaction of other constraints implies the satisfaction of aX ≤ b. A constraint system with no redundant constraint is said to be minimal.

The proliferation of redundant constraints obviously penalizes the operations on constraint systems. It is therefore essential to be able to minimize constraint systems.

Definition 7.4 (Dimension of a polyhedron).

The dimension of a polyhedron P is the dimension of the smallest affine subspace containing P . In other words, the dimension of P is equal to n, minus the number codim P of independent equations satisfied by P .

If the dimension of a polyhedron P is equal to n (i.e., codim P = 0), and if its constraint system does not contain any redundant constraint, then this constraint system is unique, up to a multiplication factor (since constraints multiple of each other are equivalent).

Definition 7.5 (Mutually redundant constraints).

Two constraints c 1 and c 2 are mutually redundant in a constraint system, if it it possible to remove either one of them -but not both -while leaving unchanged the described polyhedron. This situation may occur only if the polyhedron dimension is strictly inferior to n.

Example. In the system x 2 = 0, 0 ≤ x 1 ≤ 1, x 1 +x 2 ≥ 0, the constraints 0 ≤ x 1 and x 1 +x 2 ≥ 0 are mutually redundant (see Figure 7.2).

x 1 ≥ 0 As a consequence, when the dimension of a polyhedron if strictly inferior to n, it does not admit a unique, minimal representation as a constraint system.

x 1 + x 2 ≥ 0
If the origin is contained in the polyhedron P , it is possible to define a normal form using the scalar product. If P does not contain the origin, it can be reduced by translation. This representation will be used in Section 7. 1.5 to allow the translation from one representation to another.

Definition 7.6 (Normal form of constraints).

If the polyhedron P contains the origin, it is possible to obtain a normal form (A 1 , A 2 , A 3 ). These three sets separate the constraints according to their natures: an inequality with a positive constant, a zero constant or an equality to zero.

• The normal form of an inequation aX

≤ b with b > 0 is a b X ≤ 1 ( a b ∈ A 1 ). • The normal form of an inequation aX ≤ b with b = 0 is unchanged (a ∈ A 2 ).
• The normal form of an equation aX = 0 is unchanged (a ∈ A 3 ).

Generators

Definition 7.7 (Generating system). Any convex polyhedron P can be characterized by a generating system, constituted by three finite sets V , R and L of vectors in N n , called sets of vertices, rays and lines of P , respectively †1 . Then, P is the set of points that can be obtained as the sum of a convex combination of vertices, a positive combination of rays and a linear combination of lines:

P =    v i ∈V λ i v i + r j ∈R µ j r j + ℓ k ∈L ν k ℓ k λ i ∈ [0, 1], i λ i = 1, µ j ≥ 0, ν k ∈ N    .
Example. The polyhedron in Figure 7.1 can be represented in a dual way as a generating system (Figure 7.3), with

V = v 1 = 1 4 , v 2 = 4 1 , R = r 1 = 2 1 , r 2 = 1 0 , D = ∅.

Non-Canonicity of Generator Systems

As for constraints, the generator representation is not canonical in the general case: it may contain redundant elements (called non extremal in this context), and, when the polyhedron contains a line, it admits an infinite number of minimal generating systems.

x 1 x 2 v 1 r 1 v 2 r 2

(Extremal points and rays).

A point (respectively, a ray) in a polyhedron is extremal if it cannot be obtained by a convex combination (respectively, positive combination) of other points in P (respectively, of other rays that are not collinear to it).

The vertices v 1 and v 2 , the rays r 1 and r 2 of the example in Figure 7.3 are extremal. The point (4, 4) belongs to the polyhedron, but is not extremal, the vector (3, 1) is a non-extremal ray.

Definition 7.9 (Codimension of a polyhedron).

The codimension of a polyhedron P is the dimension of the biggest affine subspace contained in P . In other words, the codimension of P is equal to the number dim P of independent lines contained in P .

A polyhedron of codimension zero (no line) admits a unique, minimal generating system (ignoring a multiplication factor for the rays), whose vertices are the extremal points of the polyhedron, and whose rays are its extremal rays (normalized).

If the polyhedron contains a line, it has no extremal point, and no unique, minimal generating system. In this case, a minimal generating system is constituted by a minimal set of lines (a basis of the biggest affine subspace included in the polyhedron), and points and extremal rays of a section of the polyhedron: a section of P is the intersection of P with any supplementary subspace †2 of the biggest affine subspace included in the polyhedron.

In the example of Figure 7.4, the polyhedron contains the line ℓ of equation x 1 = x 2 and of direction vector (1, 1). The line of equation x 1 = 0 is a supplementary subspace to ℓ, so the segment -1 ≤ x 2 ≤ 0 is a section of the polyhedron, with vertices v 1 = (0, -1), v 2 = (0, 0) and with no ray. Thus, ({v 1 , v 2 }, ∅, {ℓ}) is a minimal generating system of the polyhedron. But we could have chosen another supplementary subspace (for instance x 2 = 0, or x 1 + x 2 = 0), so another section, and obtain another minimal generating system. If we have two lines or more, there is the choice for the basis.

Link between the Two Representations Satisfaction of Constraints by the Generator Elements

Definition 7.10. Let P be a polyhedron defined by the constraint system AX ≤ B. Then:

• Of course, a point v belongs to P if and only if Av ≤ B.

• A vector r is a ray of P if and only if Ar ≤ 0.

• A vector ℓ is a line of P if and only if Aℓ = 0. †2 Two affine subspaces E and E ′ are supplementary if their intersection is reduced to a point, and if any point of the whole space is a linear combination of points in E and E ′ .

Definition 7.11 (Polar polyhedron).

The polar polyhedron of P is the polyhedron P + defined by:

P + = {X ∈ N n | ∀Y ∈ P, XY ≤ 1}.

Properties.

• P + always contains the origin.

• If P contains the origin, then P ++ = P .

• If (V, R, L) is a generating system of P , then X ∈ P + if and only if:

∀v ∈ V, v • X ≤ 1, ∀r ∈ R, r • X ≤ 0, ∀ℓ ∈ L, ℓ • X = 0
• If P contains the origin and if (A 1 , A 2 , A 3 ) is the normal form of its constraint system (Definition 7.6), then (A 1 , A 2 , A 3 ) constitutes a generating system of P + . The converse also holds.

Considering a polyhedron representation and a characterization or an algorithm on this representation, the polar polyhedron allows to define a dual property or algorithm on the other representation. For instance, with an algorithm computing the intersection of two polyhedra given as constraint systems, we can deduce an algorithm to compute the convex hull (dual operation) of two polyhedra passed as generating systems.

Minimization of Representations

In this section, we show how the knowledge of both representations of a polyhedron can be used to minimize these representations.

Definition 7.12 (Saturation).

Let aX ≤ b be a linear inequation. Then a point v (respectively, a ray r) saturates the constraint, if and only if av = b (respectively, ar = 0). Let AX ≤ B a constraint system and (V, R, L) a generating system of the same polyhedron P . Then

• For any inequation c, Sat(c) denotes the subset of S∪R of vertices and rays that saturate c.

• For any vertex v ∈ V (respectively, for any ray r ∈ R), Sat(v) (respectively, Sat(r)) denotes the set of inequations, in AX ≤ B, saturated by v (respectively, by r).

Proposition 7.13 (Redundancy, extremality).

• A constraint c is redundant for P if and only if there exists another constraint c ′ in P such that Sat(c) ⊆ Sat(c ′ ).

• A vertex v (respectively, a ray r) is not extremal in P if and only if there exists another vertex v ′ (respectively, another ray r ′ ) in P , such that Sat(v) ⊆ Sat(v ′ ) (respectively, Sat(r) ⊆ Sat(r ′ )).

• Two constraints c and c ′ in P are mutually redundant if and only if Sat(c) = Sat(c ′ ).

The previous proposition allows to simplify the representations, and to identify mutually redundant constraints, at the cost of computing saturations (i.e., the product of the constraint matrix by the generator matrix). This computation can be complicated if there are many redundant elements, which increase the sizes of matrices. So we try to identify redundant elements earlier, when translating one representation into the other one (see Section 7.2).

In the next sections, we use the following definitions.

Definition 7.14 (Adjacency).

Let P be a polyhedron of dimension n -codim P and codimension dim P . Then

• Two vertices v and v ′ are adjacent if and only if they saturate (n -codim P -dim P -1) common inequations.

• A vertex v and a ray r are adjacent if and only if they saturate (n -codim P -dim P -1) common inequations.

• Two rays r and r ′ are adjacent if and only if they saturate (n -codim P -dim P -2) common inequations.

Computing the Dual Representation

Motzkin [START_REF] Motzkin | The Double Description Method[END_REF] has described an iterative algorithm that solves the problem of computing the dual form. A similar algorithm has been described by Chernikova [START_REF] Chernikoba | Algorithm for Discovering the Set of All the Solutions of a Linear Programming Problem[END_REF] and rediscovered in [START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF]. Thereafter, Le Verge [START_REF] Verge | A Note on Chernikova's Algorithm[END_REF] has proposed an optimization by improving the detection of redundant elements.

In practice, these algorithms do not manipulate polyhedra of dimension n but cones of dimension n+1, centered on the origin. A cone is represented by a vertex (the origin) and a set of rays and lines (seen as bidirectional rays). The polyhedron is then seen as the intersection of the cone with an hyperplane. This technical choice does not lessen the generality of the method. As for polyhedra, a cone can be described with constraints or with generators. Turning the polyhedron into a cone allows several simplifications in the process. First, by adding a dimension, vertices and rays are treated in the same way. Another advantage is that the cone always contains the origin (Section 7. 1.5) and we can settle for finding the positive solutions of a system, in which case simplex-based methods can be used.

Principle of Motzkin's Algorithm

Motzkin's algorithm can be used to compute the generating system of a polyhedron from its constraint system, as well as the inverse operation, by duality. We consider the first case.

The algorithm starts with a double, minimal description (constraints and generators) of a polyhedron that includes the polyhedron of which we want a double description. In practice, the universal polyhedron is used, whose a generating system is constituted by the origin vertex, an empty set of rays and a basis of the space as the set of lines. Initially, the set of constraints taken into account is empty.

Informally, each step of the algorithm consists in taking into account a new constraint c = (aX ≤ b) †3 . Then: †3 We present the case of an inequation, the case of equations is easily deduced.

1. If all the lines of the current generating system satisfy the constraint, the set of lines is left unchanged. Otherwise, a line ℓ not satisfying the constraint (i.e., such as aℓ = 0, see Definition 7.10) is chosen: it is composed linearly with all the other lines ℓ ′ not satisfying the constraint, to create new lines (ℓ ′′ = (aℓ ′ )ℓ -(aℓ)ℓ ′ ) that satisfy the constraint, and the vector ℓ or -ℓ, depending on the sign of aℓ, becomes a new ray r ℓ (Figure 7.5(a)).

2. For all pairs (r -, r + ) of adjacent rays in the current generating system, such that r - satisfies the constraint (ar < 0) and r + does not (ar + > 0), r + is replaced by r -= (ar + )r --(ar -)r + , a positive combination that saturates the constraint (Figure 7.5(b)).

3. For all pair (v -, r + ) (respectively, (v + , r -)) formed by a vertex strictly satisfying (respectively, not satisfying) the constraint (i.e., av -< b and av + > b, respectively) and a ray that does not satisfy it (respectively, that strictly satisfies it), (v -, r + ) (respectively,

(v + , r -)) being adjacent, r + (respectively, v + ) is replaced by the vertex v ′ = v -+ b-av - ar + •r + (respectively, v ′′ = v + + b-av +
ar -• r -) that saturates the constraint (Figure 7.5(c)). 4. For all pair of vertices (v -, v + ), adjacent and situated on either side of the constraint, v + is replaced by the vertex v ′ = λv -+ (1λ)v + , with λ = av + -b av + -av -∈ ]0, 1[, that saturates the constraint (Figure 7.5(d)). 5. Other elements in the current generating system that do not satisfy the constraint are removed. An invariant of this algorithm is that at each step, the generating system is minimal, thanks to the fact that only adjacent generating elements are combined. Therefore, the resulting generating system is minimal. Constraints are processed one by one until all constraints in the initial system have been dealt with.

r + v ′′ v + r - v ′ (c) ray + vertex → vertex c v + v ′ v - (d) vertex + vertex → vertex

Optimizations

As mentioned above, Motzkin's algorithm is, in general, implemented on a representation of polyhedra based on cones. Unfortunately, determining adjacencies is expensive, and the main optimizations have eliminated these computations:

• Wilde [START_REF] Wilde | A Library for Doing Polyhedral Operations[END_REF] and Le Verge [START_REF] Verge | A Note on Chernikova's Algorithm[END_REF] propose to combine all pairs of elements situated on either side of the constraint, without concern for adjacency. As the result of the combination is no longer necessarily extremal, its extremality is tested a posteriori with the following criterion:

Let n -codim P be the dimension of a polyhedron. Then a vertex (respectively, a ray) that does not saturate at least n -codim P (respectively, n -codim P -1) inequations is not extremal.

• The extremality condition of Le Verge is only necessary, and so does not implies the minimality of the result. It has to be ensured by a complete analysis of saturations (see Proposition 7.13). The Parma Polyhedral Library [START_REF] Bagnara | Possibly Not Closed Convex Polyhedra and the Parma Polyhedra Library[END_REF] relies on the following necessary and sufficient extremality condition:

Let n -codim P be the dimension of a polyhedron and dim P its codimension. Then a vertex (respectively, a ray) is extremal if and only if it saturates at least n -codim Pdim P (respectively, n -codim P -dim P -1) inequations.

Another optimization, used in the Linear library (Section 7. 5.4), arises from the remark that, in practice, matrices involved in representations are often very sparse, and consists in using classic data structures for sparse matrices.

Operations on Polyhedra

An abstract interpretation of interpreted automata (see Section 6.3, page 61), based on the lattice of polyhedra, requires defining the following operations:

• Lattice operations:

-Computation of the order relation (test for polyhedral inclusion).

-Comparison to ⊥, ⊤ (test for emptiness and fullness).

-Lower bound, upper bound (intersection, convex hull).

• Effect of an action, forward and/or backward.

• Widening.

These operations rely sometimes on one representation, sometimes on the other, or even both of them, confirming the interest of having both representations.

Tests

We need to perform the following tests:

Emptiness Test

A polyhedron is empty if its constraints are not satisfiable. This test can easily be performed on the minimized generating system, since a polyhedron is empty if and only if its vertex set is empty.

Fullness Test

A polyhedron is equivalent to the universe polyhedron if its constraints are always satisfied, in other words if its minimized constraint system is empty.

Inclusion Test

Here, both representations are used, since

P ⊆ P ′ ⇐⇒ ∀v ∈ V, A ′ v ≤ B ′ ∧ ∀r ∈ R, A ′ r ≤ 0 ∧ ∀ℓ ∈ L, A ′ ℓ = 0
where (V, R, L) is a generating system of P , and A ′ X ≤ B ′ is a constraint system of P ′ (see Definition 7.10).

Intersection

The intersection of two convex polyhedra is a convex polyhedron. The intersection satisfies the constraints of both polyhedra. Thus, trivially, the conjunction of constraint systems of both (or several) polyhedra is a constraint system of their intersection. This conjunction is not minimal in general and is often minimized immediately, to optimize future computations.

Convex Hull

The union of two convex polyhedra is not necessarily convex. The upper bound operation in the lattice of polyhedra is the convex hull, noted ⊔, that maps two (or several) polyhedra to the smallest polyhedron containing them (Figure 7.6).

P Q P ⊔ Q Figure 7.6: Convex hull of two polyhedra
This is the dual operation of intersection, and a generating system of P 1 ⊔ P 2 is formed by the union of vertices in P 1 and P 2 , the union of their rays, and the union of their lines. The resulting generating system is not minimal in general.

Applying Actions, Projections and Affine Transformations

The action of an interpreted automaton, as defined in Section 5. 3.1, page 52, is a relation between values of variables. Obviously, in all generality, it is impossible to define the effect of such an action on a polyhedron. Usually, it is simply defined in the case of affine relations, and roughly approximated in other cases, eliminating the involved variables (existential quantification).

Image by an Affine Relation

We only focus on the case of relations a in the form

a(X, X ′ ) ⇐⇒ X ′ = F X + G
where F is an n × n square matrix, and G is a vector of dimension n †4 . This case is interesting, from the point of view of the considered automata, since it corresponds to the actions defined by affine affectations to variables.

As in Section 5. 3.4, we note:

a(P ) = P [X := F X + G] = {F X + G | X ∈ P } and a -1 (P ) = P [X := F X + G] -1 = {X | F X + G ∈ P }.
Then, if AX ≤ B and (S, R, D) are the two representations of P :

• {F v + G | v ∈ V }, {F r | r ∈ R}, {F ℓ | ℓ ∈ L}
is a generating system of a(P ).

• AF X ≤ B -AG is a constraint system of a -1 (P ).

Existential Quantification

It is often useful to eliminate certain variables (i.e., to lose all information about them), for instance to overapproximate the effect of a non-affine affectation to these variables. This operation corresponds to an existential quantification on a polyhedron: ∃x, P . This can be done, using the generating system of P , by simply adjoining it a line which is the base vector corresponding to x, or, using the constraint system of P , by applying a Fourier-Motzkin elimination [Sch86].

Widening

The choice of a widening operator heuristic is not an absolute choice but a compromise, that requires experimental validation. The standard operator has been proposed by [START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF]. The initial idea, considering two polyhedra P and Q, with P ⊆ Q, is to build the constraint system of P ∇Q by preserving only the inequalities c = (aX ≤ b) of P verified by every point X in Q (see Figure 7.7).

x 1

x 2 P Q P ∇Q We should remember that widening is mainly used for fixed-point computation (Section 6.2.3), to guess the limit of an infinite sequence based on its first terms. Here, the idea is that a transformed inequality from P to Q (by a translation, a rotation…) may be transformed an infinite number of times and so, is removed. Polyhedra P and Q are included in P ∇Q.

As inequalities are removed from a finite system of inequalities, the algorithm stops in a finite number of iterations.

The problem with this definition is that it depends on the constraint system chosen for P : if this system is not canonical -and we stated in Section 7.1.2 that it happens when P satisfies equations -widening can unnecessarily lose information.

For instance, let us consider the following polyhedra (Figure 7.8):

• P = {x 1 = 0, 0 ≤ x 2 ≤ 1}. • Q = {0 ≤ x 1 ≤ 1, 0 ≤ x 2 ≤ x 1 + 1}.
If the widening operator is directly applied, the result is a quadrant (Figure 7.8(a)):

P ∇Q = {0 ≤ x 1 , 0 ≤ x 2 }.
x 1

x 2

P Q P ∇Q (a)
x 1

x 2

P Q P ∇Q (b)
P ∇Q = {0 ≤ x 1 , 0 ≤ x 2 ≤ x 1 + 1}.
To avoid this problem, the widening operator of [START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF] uses the notion of mutually redundant constraints (see Definition 7.5): the constraint system of P ∇Q is formed by the constraints of Q that are mutually redundant, in P , with a constraint in P . In the previous example, the constraint x 2 ≤ x 1 + 1 in Q is mutually redundant in P with the constraint x 2 ≤ 1, and the widening gives the result shown in Figure 7.8(b), whatever the shape chosen for the constraints in P .

There were subsequent work [START_REF] Bagnara | Precise Widening Operators for Convex Polyhedra[END_REF] to improve this widening operator, by distinguishing particular cases. But the standard operator, generally improved as follows, is the most commonly used.

Limited Widening

In [START_REF] Halbwachs | Verification of Real-Time Systems using Linear Relation Analysis[END_REF], a technique is proposed to improve the widening operator: when we know a set of constraints, U , that is likely to be invariant in a widening state, we can choose to keep for the widened polyhedron the constraints in U that are satisfied by all points in both widening operands. As U is finite, the fundamental properties of widening are not changed. However, this optimization often immediately obtains the results that the decreasing sequence would provide, and also gives more accurate results.

Analysis Example

Let us consider the car example first introduced in Section 5. 3.3, page 53, slightly simplified to show only numerical aspects. Each control point of the automaton is associated to a polyhedron (see Figure 7.9), as indicated in Section 6.3.

The abstract system of equations associated to this automaton is the following:

P init = ⊤ P = [d := 0][t := 0][v := 0]P init ⊔ [v := v + 1][d := d + 1](P ∩ (v ≤ 2 ∧ d ≤ 9)) ⊔ [t := t + 1][v := 0](P ∩ (t ≤ 3
t := t + 1 v ≤ 1, d < 8? → v := v + 1; d := d + 1
Figure 7.9: Car example with a polyhedron associated to each control point and P is chosen as widening point, which will be limited by the set of distinguished constraints {v ≤ 2, d ≤ 8, t ≤ 3} (see Section 7.3.6). Indeed, the control point P is attached to is exited (definitively) as soon as one of these constraints is violated, so they are a priori likely to be invariant in this control point.

The iterative resolution (Section 6.3) is the following. We only give the successive values of P until convergence, then the final value of P will allow us to compute P wall , P fast and P stop .

1. All polyhedra are initialized to the value ⊥. So, the first iteration gives:

P (1) = [d := 0][t := 0][v := 0]⊤ ⊔ ⊥ ⊔ ⊥ = (d = 0 ∧ t = 0 ∧ v = 0).

2.

The second iteration produces:

P (2) = P (1) ∇ [d := 0][t := 0][v := 0]⊤ ⊔ [v := v + 1][d := d + 1](P (1) ∩ (v ≤ 2 ∧ d ≤ 9)) ⊔ [t := t + 1][v := 0](P (1) ∩ (t ≤ 3)) = P (1) ∇ (d = 0 ∧ t = 0 ∧ v = 0) ⊔ (d = 1 ∧ t = 0 ∧ v = 1) ⊔ (d = 0 ∧ t = 1 ∧ v = 0) = P (1) ∇(0 ≤ d ≤ 1 ∧ 0 ≤ t ≤ 1 ∧ t + d ≤ 1 ∧ v = d).
The standard widening would give

P (2) = (v = d ∧ 0 ≤ t ∧ 0 ≤ d)
, but as the distinguished constraints are satisfied by both widening operands, the limited widening gives:

P (2) = (v = d ∧ 0 ≤ t ≤ 3 ∧ 0 ≤ d ≤ 2).
3. The third iteration of P is:

P (3) = P (2) ∇ (d = 0 ∧ t = 0 ∧ v = 0) ⊔ (v = d ∧ 0 ≤ t ≤ 3 ∧ 1 ≤ d ≤ 2) ⊔ (v = 0 ∧ 1 ≤ t ≤ 3 ∧ 0 ≤ d ≤ 2) = P (2) ∇(0 ≤ v ≤ d ≤ 2t + v ∧ t ≤ 3 ∧ d ≤ 2) = (0 ≤ v ≤ d ≤ 2t + v ∧ t ≤ 3 ∧ d ≤ 2).
4. The result for the fourth iteration is:

P (4) = P (3) ∇(0 ≤ v ≤ d ≤ 2t + v ∧ d -2 ≤ v ≤ 2 ∧ t ≤ 3 ∧ d ≤ 2) = (0 ≤ v ≤ d ≤ 2t + v ∧ v ≤ 2 ∧ t ≤ 3).
5. Finally, we find P (5) = P (4) ∇P (4) = P (4) .

The computation has converged towards a fixed point, that cannot be improved by a decreasing sequence. Now, we can compute the remaining polyhedra:

P wall = P (4) ∩ (d ≥ 10) = ⊥ P fast = P (4) ∩ (v ≥ 3) = (v = 3 ∧ 3 ≤ d ≤ 2t + 3 ∧ t ≤ 3) P stop = P (4) ∩ (t ≥ 4) = (t = 4 ∧ v = 0 ∧ 0 ≤ d ≤ 8)
and the result P wall = ⊥ shows that in the model, the car stops before the wall!

Existing Libraries

Several libraries, providing most of the polyhedral operators, are available; we present the main ones in this section. All libraries mentioned below operate with rational numbers (with limited or arbitrary precision), for the sake of computational accuracy.

PolyLib

Initiated at IRISA [START_REF] Wilde | A Library for Doing Polyhedral Operations[END_REF] in the early 1990s, the development of the PolyLib †5 has benefited from the cooperation of several teams, working on program parallelization and the design of systolic architectures.

The provided operations are the usual polyhedral operations, along with the possibility to manipulate explicit finite unions of polyhedra which, unlike convex hulls, do not lose information.

The PolyLib also offers sophisticated algorithms (for instance [START_REF] Nookala | A Library for Z-Polyhedral Operations[END_REF]) allowing to study integer solutions of linear systems.

The PolyLib is implemented in C.

NewPolka

The NBAC tool [START_REF] Jeannet | Partitionnement dynamique dans l'analyse de relation linéaire et application à la vérification de programmes synchrones[END_REF] uses the NewPolka library †6 , that relies on a preliminary version of the PolyLib and on optimizations proposed in [START_REF] Fukuda | Double Description Method Revisited[END_REF]. NewPolka is now integrated in the APRON library †7 , which offers access to several abstract domains, and also a higher-level interface.

Representation of coefficients

Computing with rational numbers, with reduction to a same denominator of constraint coefficients and vector components (reduction whose impact on algorithm efficiency is identified by [START_REF] Avis | Computational Experience with the Reverse Search Vertex Enumeration Algorithm[END_REF]) implies an increase of coefficient sizes, that, as soon as the dimension increases, causes arithmetic overflows. To solve this problem, NewPolka allows choosing between three integer representations: 32 bits, 64 bits and multiprecision integers using GMP †8 .

Lazy evaluation

The library keeps in general only one representation of the polyhedron. The conversion towards the other one as well as the minimization are delayed and performed only when necessary. When both representations are present, then they are both in a minimal form.

Sorting constraints and generators

Matrix lines are sorted (lazily) by lexicographic order for two reasons. In general, this sorting leads to an acceleration into the conversion algorithm. And, more specifically, it becomes easier to eliminate identical constraints, for instance in an intersection.

Strict inequalities

NewPolka allows to take into account strict inequalities, by adding an extra dimension [HPR97].

PPL

The University of Parma has also developed a polyhedral library, the Parma Polyhedral Library (or PPL) that is available on its website †9 . PPL is quite close to the library realized for NBAC.

The main difference is the programming language, C++ instead of C. Other differences derive more or less directly from the choice of C++ as programming language.

PPL is fully dynamic, no maximal size is specified during initiation.

The treatment of strict inequalities [START_REF] Bagnara | Possibly Not Closed Convex Polyhedra and the Parma Polyhedra Library[END_REF] is different from that of NewPolka. The PPL also provides alternative widening operators [BHRZ05].

PIPS/Linear

Linear is the linear library used in the PIPS compilation framework †10 [IJT91] developed at MINES ParisTech. This library uses sparse data structures to store the vectors, constraints and constraint systems, and dense data structures for matrices.

The C data type for coefficients can be set to 32 or 64-bit integers, but is incompatible with multiprecision numbers. However, intermediate computations may generate polyhedra with huge coefficients, leading to arithmetic overflows even with 64-bit integers, because constraint constants are transformed into coefficients by the convex hull operator. Eventually, when an overflow occurs, some constraints may be dropped, and the resulting invariant is less accurate than it should.

To address this problem, we added GMP support to some of the Linear polyhedral operators (Section 10. 3.2). This seemingly purely pratical implementation decision allows for a drastical simplification of the polyhedral algorithms because overflow exceptions no longer have to be handled. This is why we have developed ALICe, a toolset to compare automatic computation techniques suited for affine loop scalar invariants. It comes with a benchmark, that we built using 102 test cases found in the loop invariant bibliography, and interfaces three analysis programs, that rely on different techniques: Aspic, iscc and PIPS (see Section 8.2.5). Conversion tools are provided with ALICe to handle the various formats used by these programs.

Invariant Computation

The standard state-based model checking problem is to characterize the set of all reachable states of a transition system modeling some program. This information is generally used to check safety properties on the system, ensuring that bad configurations cannot be reached. The accuracy of computed invariants is very important, as it plays an essential role in the success of the program analysis.

Dealing with potentially infinite-state models requires to overapproximate invariants into a mathematical model (or abstract domain) such that invariant representation is finite and that allows to make the required computations. Many such domains exist in the literature, but we focus on the domain of affine invariants (Chapter 7), first introduced by N. Halbwachs [START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF][START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF], as it offers a good compromise between invariant accuracy and computational complexity.

Most of the usual analysis techniques consist in starting from a set of supposed predicates about a particular control position in the transition system, and then propagating it to other positions by evaluating the effect of each transition on the predicates. This is pretty straightforward, except in the case of loops, which requires a special treatment. This lead to intensive research, with many approaches based either on abstract interpretation [Jea00, Gon07, Mer05], that is, using widening operations, or on direct computation [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF]. The case of concurrent loops, i.e. when there are different possible transition cycles on the same control point, is particularly challenging.

As of today, several tools for linear relation analysis use a vast array of algorithms and heuristics to handle loops, aiming to maximum accuracy. We propose a toolset that compares these tools on a common set of small-scale, previously published test cases and to test the sensibility of these tools to different encoding schemes.

The ALICe Framework

The ALICe framework project aims to provide tools and a standardized set of test cases to compare as fairly as possible different polyhedral analysis techniques and programs. ALICe is a free software distributed under GPLv3, available at http://alice.cri.mines-paristech. fr/.

There are several motivations behind the ALICe project. First, we want to use it as a tool to compare polyhedral invariant computation tools on a common ground, using a set of published test cases issued from various sources, instead of evaluating their performances on ad hoc examples only, in the spirit of the TPTP library [START_REF] Sutcliffe | The TPTP Problem Library and Associated Infrastructure: The FOF and CNF Parts, v3.5.0[END_REF] and the CASC competition [START_REF] Sutcliffe | The State of CASC[END_REF][START_REF] Pelletier | The Development of CASC[END_REF].

ALICe also gives the opportunity to evaluate the benefits of model-to-model restructurations prior to analysis, and/or the sensitivity to encoding of the tools (Chapter 9).

Program Model

Test cases in ALICe -also called models -are interpreted automata as described in Section 5.3: they are transition systems constituted by a set K of control points (nodes), connected by guarded commands acting on integer variables (edges), and by a set of initial states Q init . Each model also comes with a set of error states Q err .

Definition 8.1 (Correctness).

A model m is correct if no error state is reachable from its initial states.

All models provided with ALICe are correct. An example of (correct

) model is shown in Fig- ure 8.1. k 1 k 2 t 1 : x ≥ 0 ? t 2 : x ≤ 0 ? x++ t 3 : x ≥ 1 ? x-- Initial states: Q init = {(k 1 , x) | x ∈ Val } Error states: Q err = {(k 2 , x) | x < 0}

Challenging a Tool

The purpose of error states is to challenge analysis tools. A tool successfully analyses a model if and only if the computed invariant, an overapproximation of the reachable states, does not intersect with the set of error states Q err .

More precisely, it is not possible in general to compute exactly the set of accessible states Acc(Q init ). Instead, analysis tools tested by ALICe compute supersets Q ′ of Acc(Q init ): if a given Q ′ does not intersect with Q err , then the tool that generated this Q ′ has verified that states in Q err are unreachable (Figure 8.2), thus the test case is considered a success for that tool. However, if the intersection Q ′ ∩ Q err is not empty, the tool user cannot conclude whether the property is violated or if the overapproximation is too inaccurate: this corresponds to a failure for the tool. But, as all analysis tools used within ALICe rely on forward analysis (see Section 8.2.5), this is not explored further.

Acc(Q init ) Q ′ Q err

Test Cases

The The list of models is given in Table 8. 1.

Test cases are usually relatively small, with typically 1 to 10 states and variables, and 2 to 15 transitions. Histograms showing distributions of test cases according to their sizes are shown in Figure 8. 3.

These test cases come in many forms and had to be encoded to a common format, described in Section 8.2. 4. The encoding was either done manually, from the model description, or automatically from C code using the C2fsm utility developed by P. Feautrier [Fea10, FG10].

Supported Tools

When dealing with a state transition relation, a transitive closure is useful to compute invariants [BGP97, ACI10], to check loop termination or derive loop complexities [START_REF] Gulwani | The Reachability-Bound Problem[END_REF], to build dependency constraints or to move convex array regions from a control point to another It is difficult to compare the algorithms and heuristics used by different tools designed for one of these goals because they require different inputs and produce incompatible outputs. Also, the encoding of the input often impacts the analysis.

As of now, ALICe is interfaced with three invariant computation tools: Aspic, iscc and PIPS.

• Aspic †1 [Gon13], a polyhedral invariant generator developed by L. Gonnord. Aspic relies on classic linear relation analysis, improved with abstract accelerations, identifying classes of loops in the abstract polyhedral domain whose effect can be computed directly instead of using widening operations, thus granting better accuracy.

• iscc is an interactive interface to the barvinok counting library †2 [VSB + 07, VWBC05]. barvinok is a library for counting the number of integer points in parametric and nonparametric polytopes based on isl (Section 7. 5.5) and PolyLib (Section 7.5.1). For a parametric polyhedron, the number of points is represented by either a piecewise quasipolynomial or an Ehrhart series.

• PIPS †3 [START_REF]PIPS: Automatic Parallelizer and Code Transformation Framework[END_REF], an interprocedural source-to-source compiler framework for C and Fortran programs, initiated at MINES ParisTech, that relies on a polyhedral abstraction of program behavior. Unlike other tools, PIPS performs a two-step analysis. First, the program is abstracted: each program command instruction is associated to an affine transformer representing its underlying transfer function. This is a bottom-up procedure, starting from elementary instructions, then working on compound statements and up to function definitions. Second, polyhedral invariants are propagated along instructions, using transformers previously computed. PIPS relies on the Linear library (Section

We do not consider the Omega+ [PRK + 14] library, as it appears to be superseded by isl [START_REF] Verdoolaege | isl: An Integer Set Library for the Polyhedral Model[END_REF].

It would be very interesting to use more tools, such as FASTer [Lab06, BFL04], NBAC [START_REF] Jeannet | The NBAC Verification/Slicing Tool[END_REF][START_REF] Jeannet | Dynamic Partitioning in Linear Relation Analysis: Application to the Verification of Reactive Systems[END_REF] or PAGAI [START_REF] Henry | PAGAI: A Path Sensitive Static Analyser[END_REF][START_REF] Henry | Succinct Representations for Abstract Interpretation[END_REF], but adding support is a time-consuming task, as explained below, and we have not been able to do it up to now.

Heterogeneity of Tools

ALICe test cases are written in the fsm format. This is a simple language that directly represents models, originally used by FAST [START_REF]Laboratoire Spécification et Vérification, ÉNS Cachan. FAST: Fast Acceleration of Symbolic Transition systems[END_REF] and then by Aspic. An example of fsm program for the model in Figure 8.1 is given in Listing 8. 1.

As we wish to analyze these test cases with different tools and compare results, we have to convert input and output formats. Basically, each tool uses different input and output formats:

• Aspic uses the fsm format as input and fsm expressions as output.

• iscc uses a custom format to describe both the input model as a relation on states with conditions on variables, and the output invariant, given as a map from states to polyhedral domains. An iscc relation corresponding to the model in Listing 8.1 is implemented in Listing 8.3, page 92.

• PIPS processes a structured C program according to a script written in a custom scripting language, called tpips, while resulting invariants are given as inserted before statements in the output C code. An example of PIPS output is given is Listing 10.3, page 103.

To check whether an analyzer works successfully on a test case, we follow these four steps:

1 convert a fsm model into C code or isl relation. There is also an export tool for fsm in dot format that allows visualization, fsm2dot. They are available within ALICe.

Conversely, it is possible to use ALICe to analyze simple models written in C by first translating them into fsm automata, using the C2fsm utility.

fsm2c

fsm2c is a conversion tool to generate C code from a fsm model. It is written in Python and distributed with ALICe. fsm2c works as follows:

1. The fsm code is parsed, using the PLY parsing library, and a finite state machine structure is generated.

2. A regular expression describing the possible paths in the automaton structure (regardless of transition guards) is generated, each character corresponding to a command. This regular expression represents the structure of the produced C code: unions are translated into conditional branches (if (rand()) { ... }), stars into loop statements (while (rand()) { ... }). Thus, it is important to generate a simple regular expression to avoid producing an unnecessarily complicated C code. This is done using Brzozowski's algorithm [START_REF] Brzozowski | Derivatives of Regular Expressions[END_REF], followed by a simplification step including the elimination of the ε-terms (the partial or final nature of an execution path is not relevant for invariant generation).

3. The C code is generated iteratively from the regular expression, as stated above. Each command g ? a is translated into a conditional statement with condition g and body executing the action in a. Temporary variables are introduced if needed, to handle concurrent modifications of variables in a, for instance in swaps.

Variables are renamed, to avoid any clash with C keywords. Dummy string instructions are also injected, to keep track of the variable names and control points of the original model in the produced C code. This allows to retrieve invariants at the model level when the C code is analyzed.

The C code produced from the fsm model in Listing 8.

Results for the Raw Test Cases

In this section, we present experimental results obtained with ALICe v 1.0. Benchmarks were run on a computer with a quad-core Intel i7-2600 processor clocked at 3.40 GHz with 16 GB of memory, using the following versions of analyzers (latest versions at the time of writing):

• Aspic version 3.

• iscc from barvinok version 0.37.

• PIPS revision 22 241 (September 2014).

The test suite was run 5 times and the average 3 runs were retained for each test case. Results obtained by these three tools are displayed in Table 8.2, along with the corresponding execution times. Aspic stands out as the winner in this comparison, with iscc coming second and PIPS, which is not primarily targeted at invariant analysis, being placed last, both in terms of success rate and of execution time. For analysis purposes, a model that fails to be checked can be restructured into another one, hopefully easier to analyze. Within ALICe, we have implemented two restructuring transformations on model states, to test analysis tools on a wider range of models and on specific model schemes, and explore the impact of model encoding. Model restructuring transformations are performed at the very beginning of ALICe execution, as an additional, preliminary stage, as shown in Figure 9.1.

Aspic iscc PIPS

Both restructuring transformations were proved sound in Coq, using a trace-equivalence scheme.

Definition 9.1 (Model restructuring).

A model restructuring is a sound abstraction on models, i.e. a function ρ : M → M that maps a model m 1 to a model m 2 such that: if m 2 is correct (i.e., its error region is not reachable), then m 1 is also correct:

∀m 1 , m 2 ∈ M, m 2 = ρ(m 1 ) ∧ correct(m 2 ) =⇒ correct(m 1 ). (9.1)
Thus, given a model m 1 and a restructuring ρ, it is sufficient to prove the correctness of m 2 = ρ(m 1 ) to deduce that m 1 is also correct. In addition, the restructuring ρ can be equivalent Considering an arbitrary, possible state trace with transitions

θ 1 = (k 0 , v 0 ) t 1 -→ (k 1 , v 1 ) t 2 -→ (k 2 , v 2 ) t 3 -→ • • •
in the original model m 1 , we show that for any corresponding trace θ 2 in the transformed model m 2 (corresponding to the same behavior, once the model has been transformed), then:

(for all state q 2 in θ 2 , q 2 ∈ Q err2 ) =⇒ (for all state q 1 in θ 1 , q 1 ∈ Q err1 ), thus ensuring Equation (9.1).

Control-Point Splitting Heuristic

The two main sources of imprecision that appear in Linear Relation Analysis are

• The computation of branches, performed with the convex union of invariants or transformers of each path (Sections 6.2.2 and 7. 3.3).

• The computation of loops, whether using widening (Sections 6.2.3 and 7. 3.5) or alternative techniques, as in PIPS (see Section 10.1.1).

These imprecisions are accumulated if there are structures in the automaton with two or more loops involved in the same control point. That is the case, for instance, of control k 2 in the model of Figure 8.1, page 84, that involves two self-loops t 2 and t 3 .

If such a structure is met during analysis, two ideas are used to improve accuracy. The first one is to refine invariants or transformers involved in loops, so that both the transitive closure approximation and the convex union might be more precise. The second is to reduce the number of parallel loops. Both can be achieved through a control restructuring of the analyzed automaton.

The first restructuring we use is a heuristic to split control nodes that contain several cycles. The global idea is to get rid of such nodes, which are usually the most difficult to automatically analyze, by splitting them with respect to the guards of the transitions, and adjusting the initial and error regions accordingly, as shown in Figure 9. 2. This heuristic was initially designed for PIPS and is presented in details in [START_REF] Maisonneuve | Convex Invariant Refinement by Control Node Splitting: A Heuristic Approach[END_REF].

Algorithm

The control-point splitting algorithm presented here is a restructuring aimed at achieving these two goals. For all i ∈ [1, m], add the transition (k P i , g i , a, k ′ ) with g i such that:

k 1 k ′ 2 (x ≤ 0) k ′′ 2 (x ≥ 1)
t 1 : x = 0 ? t ′ 1 : x ≥ 1 ? t 2 : x++ t 3 : x = 1 ? x-- t ′ 3 : x ≥ 2 ? x--
∀v ∈ Val, g i (v) ⇐⇒ g(v) ∧ v ∈ P i . (9.2)
3. Delete each transition (k ′ , g, a, k) that enters in k (k ′ = k). For all j ∈ [1, m], add the transition (k ′ , g, a j , k P j ) with a j such that:

∀v, v ′ ∈ Val, a j (v, v ′ ) ⇐⇒ a(v, v ′ ) ∧ v ′ ∈ P j . (9.3) 4.
Delete each transition (k, g, a, k) that loops on k. For all i, j ∈ [1, m], add the transition (k P i , g i , a j , k P j ) with g i and a j such that:

∀v ∈ Val, g i (v) ⇐⇒ g(v) ∧ v ∈ P i . (9.4) ∀v, v ′ ∈ Val, a j (v, v ′ ) ⇐⇒ a(v, v ′ ) ∧ v ′ ∈ P j . (9.5)
Controls that are not accessible nor coaccessible are not created, neither are the related transitions. Transitions whose guards are not satisfiable are not created either.

Interestingly, if guards g and actions a are polyhedral, and if every partition element P i of P is a polyhedron, then guards and actions in the resulting transitions are also polyhedral.

Correctness Theorems

We consider a general automaton m = (X, K, Trans, Q init , Q err ) and its image m P/k obtained by splitting the control node k ∈ K of m through the control-point splitting algorithm.

Theorem 9.2. For all j ∈ [1, m], for all v ∈ Val, if q = (k P j , v ′ ) is a reachable state of m P/k , then v ′ ∈ P j .

Proof. As the initial states of m P/k are the same as in m, k P j is not initial. As k P j is reachable, there exists a transition t of m P/k that enters in k P j . Let v, v ′ ∈ Val the valuations before and after entering in k P j through t. According to the algorithm above, two cases are possible:

• Either t = (k ′ , g, a j , k P j ), as described in Step 3.

• Either t = (k P i , g i , a j , k P j ), as described in Step 4.

Then a j (v, v ′ ) holds and hence, according to Equations ( 9.3) and (9.5), v ′ ∈ P j .

In other words, we have a guarantee that in every control k P i of m P/k , the invariant given by P i holds. These two theorems show that control-point splitting preserves reachable states: Theorem 9.5. Let q 1 a state of m and q 2 a state in m P/k such that q 1 ∼ St. q 2 . q 1 is reachable in m if and only if q 2 is reachable in m P/k .

Equivalence results 9.3, 9.4 and 9.5 give correspondences between automata m and m P/k in terms of traces and reachable states. Thus, safety and liveness properties on m can be translated into equivalent properties on m P/k and conversely. This allows to use the automaton m P/k instead of m to prove correctness of models, if it turns out to be easier.

Partition Choice

It is well known that the choice of a control structure affects result accuracy. Therefore, when dealing with a given model automaton m, an important question is to determine the control points that should be split as well as the partitions they should be split along, in order to make the analysis more accurate.

We have seen previously that most of the accuracy losses arise from the analysis of control points with parallel cycles. These control points are candidates to be split. Concerning the partition choices, a trade-off should be found between different criteria.

First, the automaton structure should be kept as simple as possible, in terms of control and, even more, of transition number. In the general case, partitioning a control k within m components not only adds m control states to the automaton, but also up to m 2 transitions between the newly created controls -from any control to any control -; the larger size and more complex structure will increase the analysis complexity. To avoid these issues, the number m of partition components must be bounded and the partition carefully chosen, so that some of the created states are not reachable or not co-reachable, or that some of the created transitions are not satisfiable, as they will not be present in m P/k . If possible, the priority is to eliminate transitions involved in loops or in cycles, as explained above in Section 9.1.

Also, the resulting transition guards should be as precise as possible, especially the ones involved in loops or parallel paths, to limit accuracy losses due to approximations.

Choosing the partition can be done manually, considering the system behavior. Since we want to use automatic invariant generation techniques, we also propose below an automatic heuristic technique.

Guard-Based Control Partitioning

Let m be an affine transformer automaton. On every control k with parallel cycles, let g 1 , . . . , g p be the guards of transitions looping on k.

Let G i = {v | g i (v)} and G i = Val \G i . If G i is a convex polyhedron, then G i can be partitioned into a finite set of convex polyhedrons G i,1 , . . . , G i,j i [Cha84]. So P i = {G i , G i,1 , . . . , G i,j i } is a partition of Val.
The partition taken on control point k is:

P k = P 1 ⊗ • • • ⊗ P p
where ⊗ is the product partition operation defined by:

∀E 1 , . . . , E n , ∀F 1 , . . . , F m , {E 1 , . . . , E n } ⊗ {F 1 , . . . , F m } = {E i ∩ F j | i ∈ [1, n], j ∈ [1, m]}.
In a nutshell, controls are created to explicitly allow or disallow each transition. The key idea behind this heuristic is that most of parallel loops have at least partly disjuncts guards. In this case, there are seldom controls with many parallel loops in the resulting model m P/k , if at all, and the memory state for these controls is well known (Theorem 9.2). This assumption proved reasonable in most of our test cases. Despite this, the main drawback of this technique is the important number of created controls and transitions, which limits its applicability to transition systems with a small number of transitions per control node.

In Figure 9.2, we show the effect of this heuristic applied on the model of Figure 8.

Reduction to a Unique Control Point

The other model restructuring reduces the set of controls to a unique control point ℓ, with all transitions turned into loops on ℓ, adding an extra boolean variable x k for each original control point k in both guards and actions (assuming symbols ℓ and {x k | k ∈ [1, . . . , n]} are not bound in the original model), to represent the corresponding control point of the original model. We ensure that, in every state of the system, exactly one x k is set to 1. Formally, a transition t between control states k i , k j , with guard g and action a t : k i g ? a ---→ k j is turned into the transition

t ′ : ℓ g∧x i =1∧ k =i x k =0 ? a∧x ′ j =1∧ k =j x ′ k =0
---------------------------→ ℓ.

Again, the initial and error set states are adjusted accordingly.

A more direct approach is to encode control information on a unique integer variable instead of several boolean variables, by numbering the control points in the initial model. We did not adopt it because it introduces encoding issues. Indeed, some relations on control points might or might not be expressed in terms of affine constraints, depending on the control-point encoding.

For instance, being in (k 1 or k ′ 2 ) in the model of Figure 8.1 cannot be expressed if the third control point k 2 is encoded with a value between those of k 1 and k ′ 2 . This issue is avoided with our boolean variable scheme.

The resulting model of this restructuring applied on the initial model of Figure 8.1 is shown in Figure 9. 3. The corresponding initial state set is:

Q init = {ℓ} × {b 1 = 1 ∧ b 2 = 0}. The error state set is: Q err = {ℓ} × {b 1 = 0 ∧ b 2 = 1 ∧ x < 0}. ℓ t 1 : b 1 = 1, b 2 = 0, x ≥ 0 ? b 1 = 0, b 2 = 1 t 2 : b 1 = 0, b 2 = 1, x ≤ 0 ? x++ t 3 : b 1 = 0, b 2 = 1, x ≥ 1 ? x--

Combining Restructuring Transformations

These model restructuring transformations can be used independently or together: first the model is reduced to a unique control point, then this control point is split, widening the scope of the splitting heuristic. Therefore, ALICe works with four versions of each initial model:

• The original version, with no restructuring (noted direct in the tables below);

• With all control points merged into a unique one, as described in Section 9.2 ( merged );

• Using the control-point splitting heuristic presented in Section 9.1 ( split );

• Combining both approaches ( merged-split ).

Impact of Restructuring Transformations on Experimental Results

The results obtained using these restructuring schemes are displayed in Table 9. 1.

We notice that the control state splitting heuristic leads to improved results for all tools, as shown by the comparison of Tables 9.1(a) and 9.1(c) with Tables 9.1(b) and 9.1(d), respectively. These results also confirm that iscc is significantly better in the treatment of concurrent loops, as previously noticed in Section 8.3: it outperforms the other tools on merged models

* * *

In classic Linear Relation Analysis, invariants are usually computed by propagating preconditions along the control paths of a program until stable preconditions are obtained for each control node, as seen in Section 7.4. To avoid infinite propagations, special approximation operators, e.g. widening operators, are used to guarantee the convergence within a finite number of steps.

Instead, PIPS relies on a modular alternative approach using affine transformers [IJT91], i.e., polyhedra representing transfer functions. In Section 10.1, we present how transformers then invariants are generated in PIPS. Section 10.2 gives some insight on the modularity and accuracy of the transformer-based approach. In sections 10.3 and 10.4, several refinements in transformer computation algorithms are introduced to improve invariant accuracy; they are experimentally evaluated in Section 10.5.

Generation of Invariants with Transformers

Transformers and preconditions are shown in Listings 10.2 and 10.3 for the example counter used by Halbwachs et al. [START_REF] Halbwachs | When the Decreasing Sequence Fails[END_REF], as comments just above the related statement. Transformers for elementary statements contain many trivial equations of the form x init = x because usually few variables are modified. These equations are kept implicit and instead the modified variables precondition for the next statement to improve the accuracy of T . This explains, for instance, why condition n ≤ 60 appears in the transformer of statement n++; in Listing 10.2.

Generation of Invariants

Invariants, also known as preconditions, are forward propagated from the initial state of the program using the transformers computed during the previous phase. However, accuracy is improved when some preconditions are recomputed directly for compound statements. For instance, the postcondition of a conditional if (c) T ⊤ else T ⊥ can be obtained either as the convex hull of the postconditions of the two branches:

Post = (T ⊤ • T c ) (Pre) ⊔ (T ⊥ • T ¬c ) (Pre),
or as the precondition transformed by the conditional transformer:

Post = (T ⊤ • T c ⊔ T ⊥ • T ¬c ) (Pre).
The first equation provides more accurate results at little cost, because the branch postconditions are computed anyway.

Modularity and Accuracy: Experimental Results

PIPS has been developed as a compilation framework, able to process large applications interprocedurally [START_REF] Nguyen | Efficient and Effective Array Bound Checking[END_REF]. It is important to check that the modularity provided by transformers results in the expected speed improvement and to measure the extent of its negative impact on accuracy. We show firstly that PIPS obtains accurate results in a small amount of time with respect to three other tools, Aspic, iscc and PAGAI, when dealing with loop nests and procedure calls.

We then recall previous experimental results showing a lack of accuracy and large execution times when dealing with the test cases shown in Chapter 9 to illustrate invariant generation algorithms.

Tools Used

As in Section 8.2.5, we choose to compare PIPS with Aspic and iscc using the ALICe framework.

Aspic and iscc's transitive closure functions use as input a state machine format, fsm, or a proprietary relational format. Unlike PIPS, these tools deal neither with the intricacies of C nor with their automatic abstraction. These differences in formats play an important role in the execution times obtained in Tables 8.2 and 9.1, as will be seen below.

New Improvements in Transformer Computation

We present two new improvements for transformer-based analyses. The first one deals with concurrent loops, and the second one with integer overflows.

Concurrent Loops

We are dealing with structured code. A control path set is built for each loop body. When a test or a loop is found, each pre-existing control path is duplicated to take into account the true and false branches, or the loop entrance and skip. This is not performed recursively down the branches or the inner loop bodies. Thus the total number of control path is at most 2 k , where k is the number of (possibly compound) statements in the loop body.

Control-Path Transformers in Loops

Let us assume that a loop contains several control paths, each defined by its transformer T i . Such loops are called concurrent loops above. The loop precondition, P * , can be decomposed into a set of preconditions, each obtained after a variable number of iterations and merged together by convex hull operators, ⊔: We propose to decompose the loop precondition P * as:

P * = P 0 ⊔ P 2 ⊔ P + ⊔ P 3+ (10.1)
P 0 is the precondition holding the first time the loop is entered. P 2 is obtained after two iterations with different control paths. P + corresponds to the case when only one control path is used for all iterations; it includes P 1 , the precondition of the second iteration. Finally, P 3+ is the loop precondition corresponding to all the longer control paths with at least three iterations and two different control paths.

P 2 = i j =i T i (T j (P 0 )) P + = i T + i (P 0 ) P 3+ = i j =i T + i (T j (C * (P 1 )))
C * is an over-approximation of the transitive closure of the convex hull C of the transitions, C = i T i and P 1 is the precondition after one iteration exactly, i T i (P 0 ).

When the default formula in PIPS, P * = P 0 ⊔ C + (P 0 ), is used, the convex hull operations are mostly performed before the transitive closure, in the transformer space, when computing C + . With Equation (10.1), they are executed later in the precondition space and each elementary transition T i is applied as last transition. Hence, the information brought by idempotent transformations is preserved.

The formula used to compute the loop invariant has been unrolled to take into account explicitly up to three iterations. It is possible to generalize this to k steps, with an exponential increase in the number of terms. But we are lacking experimental cases justifying such a development. [START_REF] Halbwachs | When the Decreasing Sequence Fails[END_REF] the example in Listing 10. 1. The reset to 0 is abstracted by a transformer that cannot be merged accurately with a transformer abstracting a conditional increment and no information is obtained with the equations presented in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF]. However, Equation ( 10.1) provides an accurate loop invariant, 0 ≤ n ≤ 60, by combining them, their transitive closures and the initial loop precondition.

Halbwachs et al. present in

Motivations for Equation (10.1)

Equation ( 10.1), combined with the definitions of P 2 , P + and P 3+ , was developped:

1. To transfer as much as possible convex hulls performed in the transformer space into convex hulls performed in the invariant (pre-and postcondition) space; for instance, the merge of an incrementation and a decrementation results in no information, regardless of the conditions used to control their executions.

2. C * , the transitive closure of C, the convex hull of the transformers related to each control path, is often imprecise. Information, especially idempotent components, can be restored using each control path transformer T i as last step for the iterations. This remark can be generalized to T + i . 3. In some cases, control path i can follow control path j, but j cannot follow i (T j • T i = ∅).

By considering the two last different kinds of iterations, some impossible multi-iteration paths are eliminated.

Arbitrary-Precision Numbers

With a few test cases, another issue appears with PIPS. Intermediate computations may generate polyhedra with huge coefficients, leading to arithmetic overflows even with 64-bit integers, because constraint constants are transformed into coefficients by the convex hull operator.

Eventually, when an overflow occurs, some constraints must be dropped, and the resulting invariant is less accurate than it should. To address this problem, we added GMP support to some of the PIPS polyhedral operators. This seemingly purely practical implementation decision allows for a drastical simplification of the polyhedral algorithms because overflow exceptions no longer have to be handled. It turns out that this positively impacts both the execution time (about 6 times faster now) and the comparison between transformer-and precondition-based analyses.

Other Improvements in Transformer Computation

We recall several improvements for transformer-based analyses. They can either be applied directly to the source codes, or to existing analyses. They have already been published in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF][START_REF] Maisonneuve | Convex Invariant Refinement by Control Node Splitting: A Heuristic Approach[END_REF] but they all still required implementation and experimentation.

Iterative Analysis

It is sometimes possible to improve the loop preconditions P * by recomputing the transformers a second time, using the first set of preconditions as input to limit their domains and ranges.

As explained in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF], the iterative relationship between transformers and preconditions is formalized by the two equations below where B stands for the loop body statement and the continuation condition, T for the function that converts a statement and its precondition into a convex transformer, P 0 represents the initial precondition entering the loop, and n is positive:

T n+1 = T (B, P * n ) ∩ P * n , P * 0 = ⊤, P * n+1 = P 0 ⊔ T n+1 T * n+1 (P 0 )
Note that the n-th loop precondition P * n , which can be computed more precisely using Equation (10.1), impacts the (n + 1)-th transformer T n+1 in two different ways. The affine abstrac-The example in Listing 10.7 was published in [START_REF] Dillig | Inductive Invariant Generation via Abductive Inference[END_REF]. Function foo has two different behaviors depending on the value of flag and it uses a non-affine condition in the loop, i % 2 == 0. Different behaviors that depend on a formal parameter can be separated by partitioning the function and values returned by the modulo operator can be analyzed precisely when information about the parity of i is known. When flag is not zero, it is possible to derive that a equals b. To obtain three equations in the function postcondition, transformers must be computed three times because each time a precondition makes the abstraction of a statement more precise (see Listing 10.8).

While-If to While-While Conversion

Note finally that an iterative analysis (see Section 10. 4.1) improves the invariant for the periodic function: the relationship between y and z is found (see Listing 10.10). As for the Dilig example in Listing 10.7, and unlike what is claimed in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF], the iterations are not linked only to non-polyhedral invariants, but also to polyhedral invariants and they may converge.

This optimization, the conversion of tests into while loops inside a while loop, was also used in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF] and proved correct in the extended version, but it has not been implemented, neither explicitly with a program transformation nor implicitly when computing loop invariants. In fact, it may have a detrimental effect with respect to using the convex hull to obtain a unique loop transformer, unless control path transformers are used. The example in Listing 10.11 [Mas14] shows how easier it is to prove the loop termination once the internal test has been changed into a pair of while loops.

Experimental Evaluation of the Improvements

Although PIPS has not been designed to analyze the small test cases used in academic papers dealing with loop invariants and termination and discussed in the previous chapters, it is yet interesting to use them to measure the impact of the improvements presented above and to compare PIPS to other tools able to compute invariants. Also, test cases left unsolved are good starting points for designing new improvements in invariant generation.

Impact of Improvements for PIPS

t naive = i∈ALICe time PIPS fsm2c(case i )
that is the sum of the PIPS processing time for all test cases in ALICe converted into C by fsm2c. This is highly detrimental to PIPS because of its large startup time, because number of passes, such as parsing, are not performed by other tools that are input with a simpler native format and because fsm2c has also negative impacts discussed in the next paragraph.

The C encoding chosen for fsm2c implies the parsing of the stdlib header, an interprocedural analysis by PIPS, lots of unreachable code and some exponential blow-ups in the generated code sizes. In order to reduce the overheads incurred by PIPS, we eliminate test cases which cannot be regenerated in less than 1000 lines, we replace the functional encoding of aleas by uninterpreted expressions, we use PIPS to eliminate unreachable statements and we process all test cases in one call to PIPS. Formally, the execution time is now:

t = time PIPS PIPS substitute i∈ALICe Csize(i)<1000 fsm2c(case i )
which reduces the PIPS execution time by a factor of 7 to 10.

Analysis of PIPS Failures

The ALICe benchmark contains 102 test cases, listed in Table 8.1, page 86. Out of them, there are 9 test cases that involve non-Presburger invariants and PIPS finds invariant for 82. This leaves 11 cases to investigate further, namely:

• halbwachs7.

• henry.

• metro and subway.

• microsoftex2.

• microsoftex5.

• popeea.

• realheapsort and realheapsort_step2.

• synergy_bad.

• ticket.

The next paragraphs discuss the reasons of the failures, and the improvements that can be introduced to solve them.

ALICe is based on the FAST format and the heuristics fsm2c used to generate structured C from FAST control flow graphs may blow up exponentially, up to 300 KLOC (metro, realheapsort, realheapsort_step2 and subway at the very least).

When some test cases are written in C, using while loops as they were published, the invariants for henry, halbwachs7 and synergy_bad [GHK + 06] are found by PIPS.

Case microsoftex2 is analyzed in [START_REF] Gulwani | The Reachability-Bound Problem[END_REF] to detect loop termination. The required information is not an invariant but a transformer, and PIPS computes the required transformer. Case microsoftex5 is also analyzed in [START_REF] Gulwani | The Reachability-Bound Problem[END_REF], using non-convex transformers. It contains two nested while loops, and the internal loop may or may not be the identity function. The information about the identity behavior is not preserved by PIPS because control paths are not built recursively going down loops to keep their number small. This test case requires a different algorithm to construct control paths in PIPS.

Case popeea [PC07, GZ10] has a non-convex invariant, which can only be found by a convex tool if new control points are added. The heuristics used by ALICe, fsmnodesplit, probably fails to discover a proper node splitting. Case ticket [BGP97] is interesting because it is easy to make the invariant convex. However, the number of control paths is large and the convex hulls used by PIPS lead to overflows when the while loop is unrolled. The algorithm is described as a transition system and the way it may be coded in C is critical to its analysis by PIPS.

*

Conclusion

(English version follows.)

Un nombre croissant de systèmes critiques, dont les dysfonctionnement peuvent avoir des conséquences désastreuses, sont aujourd'hui basés sur des composants informatiques. C'est le cas par exemple des commandes de vol des avions de ligne, des appareils médicaux ou des centrales nucléaires. À cause de la complexité toujours plus élevée des logiciels, la probabilité d'un vice de conception ou d'une erreur de programmation augmente. Aussi, les méthodes formelles sont de plus en plus souvent utilisées en complément des bancs de test et des normes de certification de code, pour fournir des garanties mathématiques de la correction des programmes, et ainsi atteindre une meilleure fiabilité. Ces méthodes proposent différents niveaux de précision et d'automatisation.

La vérification automatique repose en général sur le model checking et l'interprétation abstraite, une famille de méthodes efficaces pour générer automatiquement des invariants de programme, c'est-à-dire des propriétés mathématiques qui sont toujours vraies lors de l'exécution du programme, à différents endroits dans le code. Ces invariants sont ensuite utilisés pour assurer que le programme respecte sa spécification. Ainsi, il est crucial de générer des invariants précis, afin d'avoir des garanties sur les comportements complexes du programme. Les contributions présentées par ce manuscrit s'inscrivent dans ce cadre.

Contributions

Systèmes linéaires

La première partie de cette thèse concerne les sytèmes de contrôle, des dispositifs qui régulent le comportement d'autres systèmes.

* * *

Computers are involved in an ever-increasing number of critical systems, whose failures may have a disastrous impact. Examples include flight commands, medical devices or nuclear power plants. Due to the growing complexity of software, the probability of a misconception or programming error is even increasing. Hence, formal methods are involved more and more often in addition to extensive testing and code certification, to provide mathematical insurances of program correctness, and thus achieve greater reliability. These methods offer various levels of accuracy and automation.

Automatic verification generally relies on model checking and abstract interpretation, an efficient family of methods to generate automatically program invariants, i.e. mathematical properties that always hold at various locations in the code during the program execution. These invariants are then used to ensure that the program meets its specifications. Thus, it is crucial to obtain accurate invariants in order to have guarantees about complex program behaviors; this manuscript has contributed to this general subject.

Contributions Linear Systems

The first part of this thesis is about control systems, devices that regulates the behavior of other systems. Control systems are broadly utilized in industry, in a wide variety of applications. A usual topic in control theory is ensuring the stability of a control system, that is, the system started near an equilibrium point will stay forever in its neighborhood. Control theory states that a linear control system is stable if a distinctly shaped invariant, called Lyapunov invariant, holds, and many systems are designed and proved stable accordingly. But the stability proofs on these systems rely on exact arithmetic on real numbers, and do not necessarily apply when implementing the system on a device with machine-precision arithmetic.

We introduce a new framework to translate real-arithmetic proof invariants into invariants suited for machine arithmetic, taking into account number representation and rounding error issues. The generated invariants are guaranteed to be correct, but do not necessarily ensure the system stability. This generic translation scheme is checked using Coq, and is also implemented in the case of floating-point arithmetic, the most common standard for machine arithmetic. This implementation, LyaFloat, can check automatically system stability on floating-point platforms, and compute the number of bits necessary to ensure it, for both open-loop and closed-loop stabilities.

Linear Relation Analysis

The second part of this manuscript is about linear relation analysis, a classical abstract interpretation based on an over-approximation of reachable numerical states of a program by convex polyhedra. As of today, several tools for linear relation analysis exist. They rely on a large set of algorithms and heuristics to handle loops, aiming to maximum accuracy. We propose a new toolset, ALICe, to compare as fairly as possible three of these tools: Aspic, iscc and PIPS, on a common set of 102 small-scale, test cases previously published in the loop invariant and termination bibliography. Two main original results arise from the experimental measurements.

First, the choice of encoding is critical to the success of a model analysis. This issue is explored with two model-to-model restructuring algorithms that we have developed and proved correct with Coq. The first one reduces the encoding biases, by merging the model to a unique control point. The other one splits control nodes involving concurrent loops that are difficult to analyze. Used separately or together, they lead to significant improvements in analysis results and help to identify the strengths and weaknesses of each analysis tool.

The second experimental observation is that PIPS default algorithms have poor accuracy when dealing with concurrent loops. Unlike Aspic and iscc, invariant computation in PIPS relies on a first step that computes transformers, used later to propagate invariants. The greater scalability of this approach is useful to analyze large applications, at the cost of additional approximations.

We present new improvements in transformer computation to improve invariant accuracy in PIPS. Combined with other existing techniques, PIPS is able to perform at comparable levels with Aspic and iscc.

Perspectives Linear Systems

Regarding perspectives, several improvements could be introduced in our invariant translation tool LyaFloat. Currently, only the most common floating-point arithmetic is supported; yet it would be relatively easy to add support for fixed-point arithmetic, at least from a theoretical standpoint (the implementation raises technical difficulties, due to the lack of a suitable parametric fixed-point library for Python). In the same vein, it would be interesting to add support for more functions and propagation arguments, in order to extend the scope of the tool.

In its current design, LyaFloat generates invariants for machine arithmetic code following a formally proved pattern but there is no formal guarantee that LyaFloat itself is flawless. Thus, we would gain extra confidence if we could check the generated invariants by some other means. For this, LyaFloat could be rewritten and formalized in Coq rather than Python. This would require to prove (or at least formalize) the invariant propagation arguments needed for stability proofs, which could be a huge work. An alternative idea is to make LyaFloat generate invariants along with Coq scripts acting as certificates that the invariant propagation is correct, which would be sufficient to easily verify the generated invariants. This would require the development of a Coq library of theorems linked to stability theory.

From a general point of view, many steps are required to translate an idealized, proved linear controller model into the low-level C code that implements it, and shifting to the arithmetic representation is just one of them. Other outstanding issues include the discretization of continuous controllers, which poses the problem of guaranteeing a bounded computation time in the control loop, or the data acquisition, which is implemented on most controllers as interrupt handlers.

Linear Relation Analysis

In ALICe, dealing with C source code has proven problematic. Currently, all test cases are implemented as finite state machines in fsm format and compiled in C prior to being analyzed by PIPS. The translation from fsm to C is not straightforward and the generated C files tend to be overly long and complex: this is not optimal for analysis purposes, and is a disadvantage for PIPS in comparisons with Aspic and iscc. Even if the problem can be partially circumvented with an intermediate simplification step, we believe a better strategy should be employed to handle C code, either by storing a C version of each model alongside with the fsm version (which brings several other problems concerning representation choices, expressivity, consistency and restructuring), or by drastically improving the fsm2c tool to tackle these issues.

The scope of ALICe could be enhanced too, by collecting more models in the test suite on the one hand, and by supporting more analysis tools on the other. The later point would be timeconsuming because of the lack of standardization of formats used by analysis tools, requiring conversion routines for each added analysis tool. We should also check that the analysis tools generate realistic invariants, for instance shipping each model with minimal invariants. An analysis tool would be considered as correct only if the generated invariants are superset of the minimal ones. This would prevent a trivially faulty tool to trash the results.

Finally, even if the improvements made to transformer loop computations allowed to significantly improve the analysis results in PIPS, their computational complexity makes them unsuited for large code analysis. Accurate algorithms with a decent complexity class still remain to be found.

Analyse statique des systèmes de contrôle-commande : invariants entiers et flottants

Résumé : Un logiciel critique est un logiciel dont le bon fonctionnement a un impact important sur la sécurité ou la vie des personnes, des entreprises ou des biens. L'ingénierie logicielle pour les systèmes critiques doit donc répondre à des exigences drastiques de qualité, souvent réglementaires, permettant d'avoir un degré de confiance très élevé dans les programmes. C'est une tâche particulièrement difficile, qui combine différentes méthodes pour garantir la qualité des logiciels produits. Parmi celles-ci, les méthodes formelles peuvent être utilisées pour prouver qu'un logiciel respecte ses spécifications. Le travail décrit dans ce manuscrit s'inscrit dans le contexte de la validation de propriétés de sûreté de programmes critiques, et plus particulièrement des propriétés numériques de logiciels embarqués dans des systèmes de contrôle-commande. La première partie de cette thèse est consacrée aux preuves de stabilité au sens de Lyapunov. Ces preuves raisonnent sur des calculs en nombres réels ; or, ceci ne décrit pas correctement le comportement d'un programme exécuté sur une plateforme à arithmétique machine. Nous présentons un cadre théorique générique pour adapter les arguments des preuves de stabilité de Lyapunov aux arithmétiques machine. Un outil effectue automatiquement la transposition de ces preuves en nombres réels vers des preuves en nombres à virgule flottante.

La seconde partie du manuscrit porte sur l'analyse des relations affines, une interprétation abstraite qui approxime les valuations associées aux points de contrôle d'un programme par des polyèdres convexes. Nous présentons ALICe, un framework permettant de comparer différentes techniques de génération d'invariants. Il s'accompagne d'un benchmark composé d'une collection de cas de test tirés de publications sur l'analyse de programmes, et s'interface avec trois outils utilisant différents algorithmes de calcul d'invariants : Aspic, iscc et PIPS. Pour tester la sensibilité des outils à l'encodage, deux restructurations conservatives des cas de test sont proposées. Afin d'affiner les invariants générés par PIPS, nous introduisons dans un dernier chapitre plusieurs améliorations à ses algorithmes.

Mots clés : analyse statique, analyse des relations linéaires, calcul en virgule flottante, sécurité et sûreté, stabilité de Lyapunov

Static Analysis of Control-Command Systems:

Floating-Point and Integer Invariants Abstract: A critical software is a software whose proper functioning has an important impact on safety or life of people, companies and goods. Software engineering on critical systems must comply with drastic quality requirements, often regulatory, to achieve a higher level of trust in programs. This task is particularly difficult and combines different methods to ensure the quality of produced software. Among them, formal methods can be used to prove that a software obeys its specifications. The context of the work described in this manuscript is the validation of safety properties for critical software; specifically, of numerical properties for embedded software in control-command systems.

The first part of this thesis is dedicated to Lyapunov stability proofs. These proofs reason with exact real-number computations, but this is not an accurate description of how a program behaves on a machine-arithmetic platform. We present a generic, theoretical framework to adapt Lyapunov stability proof arguments to machine arithmetics. A tool automatically transposes these proofs on real numbers to proofs on floating-point numbers.

The second part of the manuscript is about linear relation analysis, an abstract interpretation that approximates valuations associated to a program's control points by convex polyhedrons. We present ALICe, a framework to compare different invariant generation techniques. It comes with a benchmark consisting of a collection of test cases drawn from the program analysis literature, and interfaces with three tools relying on different invariant computation algorithms: Aspic, iscc and PIPS. To evaluate the sensitivity of these tools to encoding choices, two techniques are proposed to conservatively restructure test cases. To refine invariants generated by PIPS, several improvements are made on its algorithms and introduced in the last chapter.

Figure 1 . 1 :

 11 Figure 1.1: Centrifugal governor at the Science Museum of London

Figure 1 . 3 :

 13 Figure 1.3: Closed-loop stability

Figure 2 . 1 :

 21 Figure 2.1: Mass-spring system

1 :

 1 Pseudocode of the controller

Figure 3 . 1 :

 31 Figure 3.1: Abstract scheme of invariant propagation on the original program

Figure 3 . 2 :

 32 Figure 3.2: General translation scheme

Figure 3 . 3 :

 33 Figure 3.3: Refined translation scheme

  Let d ∈ D F be the precondition corresponding to d after converting both constants and functions in the preceding instructions of the program. Let d′ = p( d, ĩ) ∈ D R . According to Equation (3.2), d′ is a valid postcondition for precondition d and instruction ĩ.

Now, let usTheorem 3 . 1 .

 31 consider an instruction ĩ in I R surrounded by a precondition d and a postcondition d′ = p( d, ĩ) in D R . We have: eval(ĩ, d) ⊂ d′ Let v err ∈ Val F be an upper bound of the arithmetic error induced by instruction ĩ on domain d. Then, the implication

#

  Definition of E P P = Rational("1e-3") * Matrix(rationals( ["0.6742 0.0428", "0.0428 2.4651"])) EP = Ellipsoid(P) # Definition of E Qµ mu = Rational("0.9991") Qmu = mu * P Qmu = Qmu.col_insert(2, zeros(2, 1)).

  constants( ["0.4990 -0.0500", "0.0100 1.0000"])) Bc = Matrix(constants(["1", "0"])) Cc = Matrix(constants(["564.48 0"])) Dc = Matrix(constants(["-1280"])) # Definition and verification of E R AcBc = Ac.col_insert(Ac.cols, Bc) R = (AcBc * Qmu.inv() * AcBc.T).inv() ER = Ellipsoid(R) print("ER included in EP :", ER <= EP) # Computation and verification of E R i = Instruction({Xc: Ac * Xc + Bc * Yc}, pre=[Zc in EQmu], post=[Xc in ER]) ERbar = i.post()[Xc] print("ERbar =", ERbar) print("ERbar included in EP :", ER <= EP) Listing 4.1: Computation of E R using LyaFloat

3 . 1 , 3 :

 313 the first step when translating Feron's program is to convert real constants into floating-point numbers. The first lines of code in the program: 1 Ac = [0.4990, -0.0500; 0.0100, 1.0000]; 2 Bc = [1; 0]; 3 Cc = [564.48, 0]; 4 Dc = -1280; Listing 4.2: Original constant definitions in the program become, assuming rounding to nearest value: Constant definitions, using floating-point values Theses matrices A c , B c , C c and D c will be used instead of the original matrices A c , B c , C c and D c in the sequel of the proof.

9991 10 uListing 4 . 4 :

 1044 = Cc*xc + Dc*yc; Beginning of the proof, unchanged

5 :

 5 Computation of u, original code First of all, matrices C c and D c must be replaced by their floating-point counterparts C c and D c both in the program instruction Line 10 and the ensuing invariant. This invariant relies only on algebraic arguments and does not depend on the values in the matrices, it still holds considering exact arithmetic operations. But this is not sufficient: indeed, this instruction is a sum of matrix multiplications, i.e. a set of additions and multiplications on floating-point numbers that yield rounding errors (see Section 3.3.2).

00247323x c 2 2 ≤

 2 1} we get the following values: α = 1.03246 • 10 -17 , β = 1.84829 • 10 -17 , γ = 7.17582 • 10 -32 . Our program finds δ = 5.35754 • 10 -16 ≫ γ and finally k = 1 + 5.35754 • 10 -16

c 2 :Figure 5 . 2 :

 252 Figure 5.2: An Interpreted Automaton

Figure 5 . 3 :

 53 Figure 5.3: Interpreted automaton with boolean input simulating a car controller

Figure 5 . 4 :

 54 Figure 5.4: The Car Example

Figure 5 . 5 :

 55 Figure 5.5: An Interpreted Automaton

  Figure 6.1: Galois Connection

Figure 6 Figure 6 . 3 :

 663 Figure 6.2: The Sign Lattice

Figure 6 . 4 :

 64 Figure 6.4: Abstraction Using Affine Equations

Figure 6 . 5 :

 65 Figure 6.5: Abstraction Using Linear Congruences

Figure 6 . 7 :Figure 6 . 8 :

 6768 Figure 6.7: Abstraction Using Polyhedra

Figure 6 . 9 :

 69 Figure 6.9: Abstraction with Octagons

Figure 7 . 1 :

 71 Figure 7.1: Representation of a polyhedron as a constraint system

Figure 7 . 2 :

 72 Figure 7.2: Mutually redundant constraints

Figure 7 . 3 :

 73 Figure 7.3: Representation as a generating system

2 Figure 7 . 4 :

 274 Figure 7.4: Polyhedron containing a line

Figure 7 . 5 :

 75 Figure 7.5: Steps of Motzkin's algorithm

Figure 7 . 7 :

 77 Figure 7.7: Standard widening for convex polyhedra

Figure 7 . 8 :

 78 Figure 7.8: Widening a polyhedron satisfying an equation

3 t ≤ 2 ?

 32 ))P wall = P ∩ (d ≥ 10) P fast = P ∩ (v ≥ 3) P stop = P ∩ (t ≥ 4) → v := 0;

Figure 8 . 1 :

 81 Figure 8.1: An example of model

Figure 8 . 2 :

 82 Figure 8.2: Model checking

  benchmark suite itself consists in 102 previously published test cases, including work by L. Gonnord [Gon07, ADFG10], S. Gulwani [CCG + 08, Gul09, GJK09, GMC08, GZ10], N. Halbwachs [Hal93, HPR97, Hal10, Hal79], B. Jeannet [Jea00] et al. The comprehensive list of model sources can be found in the bibliography [Bar05, Lee02, BHMR07, BGP97, CCG + 08, CS01, CS02, CPR06, Cou05, GR06, GR07, GHK + 06, Hen11, HJMS02, KEW + 85, Mer05, PR04, PC07], as well as on the ALICe website. They come mostly from works on loop invariant computation, loop bound analysis, loop termination and, to a lesser extent, protocol verification.

Figure 8 . 3 :

 83 Figure 8.3: Distribution of test cases

Figure 9 . 2 : 2 .

 922 Figure 9.2: Model of Figure 8.1 transformed by the splitting heuristic

Figure 9 . 3 :

 93 Figure 9.3: Model of Figure 8.1 reduced to a unique control point This restructuring has three purposes. First, it stresses the tool with more difficult test cases. It also reduces bias factors related to encoding choices. Finally, if used before the control-point splitting heuristic, merging increases the effect of splitting by generating the type of controlpoint structure to which it applies.

Keywords:

  Static Analysis, Linear Relation Analysis, Floating-Point Computation, Security and Safety, Lyapunov Stability

  Both properties are interesting: closed-loop stability offers guarantees on the plant state, whereas an open-loop stable controller is ensured not to diverge disastrously if there is a temporary malfunction in its input sensors. This is why both open-loop and closed-loop stability are often required simultaneously in critical systems. Both types of stability are considered in this dissertation.La stabilité constitue un attribut essentiel des systèmes de contrôle, en particulier lorsque la sécurité des personnes est en jeu, comme c'est le cas par exemple en médecine ou en aéronautique Cette partie de la thèse est organisée comme suit. Dans le reste de ce chapitre, nous décrivons l'exemple de système dynamique du second ordre utilisé par Feron dans[START_REF] Feron | From Control Systems to Control Software[END_REF], avec le contrôleur correspondant. Nous présentons ensuite l'analyse par Feron du contrôleur en boucle ouverte avec des nombres réels. Dans les chapitres qui suivent, nous présentons notre schéma générique de traduction (chapitre 3, page 29) et son implémentation dans une bibliothèque Python, et nous l'utilisons par la suite pour réécrire l'analyse de Feron en arithmétique à virgule flottante (chapitre 4, page 37). Enfin, on traitera le cas du système en boucle fermée. Cette partie se termine par une discussion sur la généralité de cette approche.Stability constitutes an essential attribute of control systems, especially when human safety is involved, as in medical or aeronautical domains (Chapter 1). Modern system developments, such as adaptive control technologies, rely on robust stability and performance criteria as the primary justification for their relevance to safety-critical control applications. Motivated by such applications, there exist many theorems that ensure system stability and performance under various assumptions and in various settings. Lyapunov's stability theory[START_REF] Lyapunov | The General Problem of the Stability of Motion[END_REF] plays a critical role in that regard.

	Chapter 2
	Lyapunov Stability Proofs * * *

(chapitre 1). Les développements récents, comme les technologies de contrôle adaptatif, utilisent principalement des critères de robustesse, de stabilité et de performance pour justifier leur pertinence dans des applications de contrôle critiques. Il existe de nombreux théorèmes, motivés par de telles applications, qui assurent la stabilité et la performance d'un système sous diverses hypothèses et dans divers contextes. La théorie de stabilité de Lyapunov

[START_REF] Lyapunov | The General Problem of the Stability of Motion[END_REF]

joue à cet égard un rôle essentiel. L'implémentation logicielle de bas niveau d'une loi de contrôle peut être inspectée par les outils d'analyse disponible pour aider au développement de programmes informatiques critiques. La technique d'analyse de programme la plus simple consiste à effectuer plusieurs simulations, en incluant dans certains cas une représentation logicielle ou matérielle du système asservi dans la boucle. Cependant, les simulations ne peuvent fournir d'information que sur un nombre limité de comportements du système. Les méthodes plus avancées comprennent le model checking et l'interprétation abstraite, par exemple en utilisant Astrée [CCF + 15]. Ces méthodes prennent en entrée les programmes informatiques eux-mêmes et renvoient en sortie des certificats garantissant le bon comportement des programmes par rapport aux spécifications choisies. Une autre possibilité est d'utiliser des techniques de démonstration assistée de théorèmes, proposées par des outils comme Coq, Isabelle ou PVS [INR15, PNW14, Com15]. Ces assistants de preuve peuvent être utilisés pour établir des propriétés de programmes et plus généralement des constructions mathématiques. Le model checking, l'interprétation abstraite et les assistants de preuve sont tous trois utilisés pour vérifier des applications critiques.

Dans

[START_REF] Feron | From Control Systems to Control Software[END_REF]

, Feron étudie comment des informations de domaine sur les système de contrôle -et, en particulier, les preuves de stabilité de Lyapunov sur une modélisation théorique de haut niveau d'un système -peuvent être transposées vers des certificats de comportement de logiciels de contrôle lisibles et vérifiables par un ordinateur, en s'appuyant sur le système de preuve de Floyd et Hoare

[START_REF] Peled | Software Reliability Methods[END_REF]

, appliqué à du pseudo-code MATLAB (voir les sections 2.2 et 2.3). Feron présente des implémentations en boucle ouverte et en boucle fermée de son contrôleur, selon que l'action du système sur son environnement et la rétroaction sont modélisées ou non. Mais les erreurs causées par l'arithmétique à virgule flottante ne sont pas considérées.

Dans cette partie, nous abordons cette question en présentant une approche automatique pour traduire des invariants de preuves de stabilité de Lyapunov sur du pseudo-code à arithmétique réelle, comme celui fourni dans l'article de Feron, vers des invariants similaires sur du code machine en prenant en compte les erreurs d'arrondi introduites par l'arithmétique à virgule flottante. Nous utilisons ces invariants pour vérifier si les conditions de stabilité sont toujours respectées, auquel cas la stabilité du système à virgule flottante est établie. Cette approche est basée sur un framework théorique générique, décrit plus loin dans ce document, dont la sûreté est prouvée en Coq ; elle est également implentée sous la forme d'un outil automatique.

The low-level software implementation of a control law can be inspected by analysis tools available to support the development of safety-critical computer programs. The simplest program analysis technique consists in performing several simulations, sometimes including a software or hardware representation of the controlled system in the loop. However, simulations provide information about only a finite number of system behaviors. More advanced methods include model checking and abstract interpretation, e.g. using Astrée [CCF + 15]. In these methods, inputs are the very computer programs and outputs are certificates of proper program behavior along the chosen criterion. Another possibility is to use theorem-proving techniques, supported by tools such as Coq, Isabelle or PVS [INR15, PNW14, Com15]. These proof assistants can be used to establish properties of programs and more general mathematical constructs. Model checking, abstract interpretation, and theorem-proving tools are all used to verify safety-critical applications.

zeros(2,1); % x c ∈ E P

	x c2
		20
	E P	10
	-40 -30 -20 -10	10	20	30
	-10
	-20
	Figure 2.2: The stability domain E P
	A stability proof of the controller code is provided in [Fer10], using Floyd-Hoare program annotation technique [Hoa69, Flo93, Pel01]: each program instruction comes with an invariant. The program annotated by Feron is reproduced below.
	6 receive(y, 2); receive(yd, 3); % x c ∈ E P 7 while (1) % x c ∈ E P 8 yc = max(min(y -yd, 1), -1); % x c ∈ E P , y 2 c ≤ 1
	9

5 xc =

  Domains are used to represent Floyd's and Hoare's invariants seen in previous sections. A domain d on K is a set of valuations on K, i.e. d ⊂ Val K . The set of domains on K is noted D K , with D F ⊂ D R . A real domain d ∈ D R can be associated to a limited-precision domain d ∈ D F defined by:

  cri.mines-paristech.fr/ people/maisonneuve/thesis/LyaFloat.v.

	Chapter 4
	Translating Stability Proof Scheme to Floating-Point Numbers

Nous voulons vérifier si la preuve de stabilité dans le programme de Feron est toujours valable sur un contrôleur implémenté avec des nombres flottants, en utilisant l'approche présentée dans le chapitre précédent. Dans ce chapitre, on utilisera le standard IEEE pour l'arithmétique à virgule flottante [IEE08], codée sur 64 bits †1 , comme c'est le cas dans la majorité des unités de calcul à virgule flottante actuelles. Dans cette norme, l'addition et la multiplication sont correctement arrondies en fonction du mode actif d'arrondi, ce qui permet de borner l'erreur en fonction de la valeur des opérandes. Le cas des représentations numériques alternatives est discuté dans la section 4.

5.

  Dans ce chapitre, nous donnons les principales définitions et notations relatives aux programmes, leur sémantique, et les propriétés que nous voulons vérifier les concernant. Tout d'abord, nous introduisons les systèmes de transition (section 5.1) qui constitueront le modèle sémantique standard des programmes ; nous définissons ensuite les propriétés de ces systèmes, et la problématique générale de la vérification (section 5.2). Les automates interprétés, un modèle de programme très général, sont introduits dans la section5.3.

	Chapter 5
	Models, Programs, Verification

  Comme indiqué en introduction, une partie de ce travail est consacrée à l'analyse des relations linéaires, basée sur le treillis des polyèdres convexes. Ce chapitre présente les principaux algorithmes sur les polyèdres. Les algorithmes utilisés reposent sur deux représentations des polyèdres -par système de contraintes et par système générateur -présentée dans la section 7.1.

La section 7.2 détaille le passage de l'une à l'autre des représentations, ce qui est l'algorithme clé pour effectuer toutes les autres opérations et surtout pour simplifier les représentations. La section 7.3 décrit les opérations usuelles sur les polyèdres convexes. Un exemple d'analyse est ensuite détaillé dans la section 7.

4

. Finalement, les principales bibliothèques polyédriques existantes sont présentées dans la section 7.

5

.

Les références générales de ce chapitre sont

[START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF]

,

[START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]

et

[START_REF] Jeannet | Partitionnement dynamique dans l'analyse de relation linéaire et application à la vérification de programmes synchrones[END_REF]

.

* * *

As mentioned in Introduction, part of this work is devoted to the analysis of linear relations. It is based on the lattice of convex polyhedra. This chapter presents the main algorithms relative to polyhedra. The algorithms we use are based on the two characterizations of polyhedra -by constraint systems and by generating systems -that are presented in Section 7.

1. Section 7.2

  Comme cela a déjà été souligné au chapitre 6, un point crucial de l'analyse de programme est le calcul des invariants de boucle. Il est nécessaire d'avoir des invariants précis pour prouver des propriétés sur un programme, mais de tels invariants sont difficiles à calculer. Des recherches approfondies ont été effectuées mais, à notre connaissance, aucun benchmark n'a été développé jusqu'à présent pour comparer les algorithmes et les outils.

	Chapter 8
	ALICe: A Framework to Improve Affine Loop Invariant Computation

C'est pourquoi nous avons développé ALICe, un ensemble d'outils pour comparer des techniques de calcul automatique d'invariants affines de boucle. ALICe s'accompagne d'un benchmark, que nous avons construit à partir de 102 cas de test trouvés dans la littérature sur les invariants de boucles, et s'interface avec trois programmes d'analyse, reposant sur différentes techniques : Aspic, iscc et PIPS (voir la section 8.2.5). Des outils de conversion sont fournis pour manipuler les différents formats utilisés par ces programmes.

* * *

As already underlined in Chapter 6, a crucial point in program analysis is the computation of loop invariants. Accurate invariants are required to prove properties on a program but they are difficult to compute. Extensive research has been carried out, but, to the best of our knowledge, no benchmark has ever been developed to compare algorithms and tools.

Table 8 .

 8

	aaron2 ackerman ax bakery bardin berkeley car car_simple counterex1 counters1 counters2 counters3 counters4 cousot9 determinant disj disjbnd1 disjbnd2 easy1 easy2 exmini forward gasburner gasburner_alt gcd gonnord1 gonnord2 gonnord3 gonnord4 gonnord5 gonnord6 gopan_reps gopan_reps_alt1 gopan_reps_alt2	gulwani1 gulwani1_alt gulwani2 gulwani2_alt gulwani3 gulwani4 halbwachs1 halbwachs2 halbwachs3 halbwachs4 halbwachs5 halbwachs6 halbwachs7 halbwachs8 halbwachs9 henry insertsort interleaving1 interleaving2 interleaving3 interleaving4 jeannet1 jeannet2 jeannet3 jeannet4 jeannet5 jeannet6 leeyannakis_bad loops maccarthy91 metro microsoftex2 microsoftex4 microsoftex5	microsoftex6 microsoftex7 multcounters1 multcounters2 multcounters3 multcounters4 nd_loop ndecr nested perfect popeea random1d random2d realbubble realheapsort realheapsort_step1 realheapsort_step2 realselect realshellsort relation1 rsd sipmabubble slam slam_bad speedometer subway synergy_bad terminate terminator ticket wcet1 wcet2 while2 wise

[START_REF] Khaldi | Parallélisation automatique et statique de tâches sous contraintes de ressources : une approche générique[END_REF]

. When dealing with a dependence relation, a transitive closure is useful to optimize or synthesize code [Won01, VCB11, BKK13]. 1: List of models currently provided with ALICe

Table 8 .

 8 2: Benchmark resultsThis ranking is not a total order because no tool is strictly better than another: for each tool, there is at least one model that is successfully analyzed only by this tool, as shown in Figure8.5; and it appears in Table8.3 that there is no clear trend as for the quality of generated invariants, in terms of invariant inclusion.

	Successes Time (s.)		75 10.9 35.5 43.3 63 43
	successes: 78 failures: 24		Aspic 12
			22	39	2
	iscc	1	1	1	PIPS
	Figure 8.5: Venn diagram of successes for each tool
	⊇ Aspic iscc PIPS	Aspic iscc PIPS -21 23 49 -54 32 23 -

Table 8 .

 8 3: Invariant inclusions, in terms of setsA closer analysis shows that iscc performs comparatively well on test cases encoded with concurrent loops (several loops on a single control point, similar to what is shown Figure9.3, page 98), unlike PIPS whose raw results appear to be the poorest. On the other hand, iscc can be quite slow on test cases that display a large, intricate control structure. Finally, despite its successes, Aspic has greater difficulty to deal with transitions featuring complex formulas, that it is not able to accelerate. This leads us to wonder if the tools are sensitive to the encoding, since a problem can be presented under many guises. This is the topic of Chapter 9.PIPS comparative slowness is primarily due to the fact that it uses C as input format, which is much more universal and thus, difficult to handle, than Aspic or iscc specific formats. This issue is revisited in Section 10.2.1.On mentionne au chapitre 8 la possibilité d'utiliser ALICe pour évaluer l'efficacité et la pertinence des restructurations de modèles. La motivation initiale était d'améliorer des résultats obtenus avec PIPS en jouant sur la structure de contrôle des modèles.À des fins d'analyse, un modèle que l'on ne parvient pas à vérifier peut être restructuré en un autre modèle, que l'on espère plus facile à analyser. Dans ALICe, nous avons implémenté deux restructurations sur les états des modèles, pour tester les outils d'analyse sur un plus grand échantillon de modèles et sur des schémas particuliers de modèles, et explorer l'impact de l'encodage des modèles. Les restructurations de modèles sont effectuées au tout début de l'exécution d'ALICe, lors d'une étape préliminaire additionnelle, représentée sur la figure 9.1.We mention in Chapter 8 the possibility to use ALICe to test the efficiency and relevance of model-to-model restructuring transformations. An initial motivation was to improve results in PIPS by playing on the model control structure.

	Chapter 9
	Model-to-Model Restructuring Transformations: Sensitivity to Encoding
	La sûreté de ces deux restructurations a été prouvée en Coq, en utilisant un raisonnement par équivalence de traces.
	* * *

  is correct if and only if m 2 is correct), although we are not interested in proving such properties in our toolchain. So, for instance, the restructuring that to any model associates the same trivial, inconsistent model is sound (but not very interesting).

			Aspic	out.asp	out.is	iscc
	in0.fsm restructuring in.fsm	in.is	iscc	out.is	iscc
		in.c	PIPS	out.c	out.is	iscc
	Figure 9.1: Analysis steps within ALICe, using restructuring transformations
	(m 1

  Given two controls k 1 in m, k 2 in m P/k , we introduce a state equivalence relation ∼ St. between two states of m and m P/k , defined by:(k 1 , v) ∼ St. (k 2 , v) if and only if • Either k 1 = k and k 2 = k P i where i ∈ [1, m] satisfies: v ∈ P i . • Either k 1 = k 2 = k.We also define a relation ∼ Tr. between traces of m and m P/k : T 1 ∼ Tr. T 2 if and only if states in T 1 and T 2 pairwise satisfy ∼ St. . Both relations ∼ St. and ∼ Tr. are bisimulations. They are one-to-one mappings, so notions of image and inverse image by ∼ St. and ∼ Tr. are defined. For all trace T 1 of m, there exists a trace T 2 of m P/k such that T 1 ∼ Tr. T 2 . For all trace T 2 of m P/k , there exists a trace T 1 of m such that T 1 ∼ Tr. T 2 .

	Then, the following two theorems hold:
	Theorem 9.3. Theorem 9.4.

  1, where the control point k 2 is split into two components with respect to the guards of transitions t 2 and t 3 . The initial state set is unchanged. Proving that the new error state set {k ′ 2 , k ′ 2 } × (x < 0) cannot be reached is easier than in the original model: it is true by construction for the control point k ′′ 2 , and can be easily deduced by looking at guards entering transitions for control point k ′ 2 .

  Dans l'analyse classique des relations linéaires, les invariants sont normalement calculés par propagation des préconditions le long des chemins de contrôle du programme, jusqu'à obtenir une précondition stable sur chaque point de contrôle, comme nous l'avons vu à la section7.4.

	Chapter 10
	Computing Invariants with Transformers: Experimental Scalability and Accuracy
	1 22 Aspic iscc PIPS 1 1 2 39 successes: 78 failures: 24 ⊇ Aspic -21 23 iscc 49 -54 PIPS 32 23 -(a) Direct ⊇ Aspic iscc PIPS Aspic -7 8 iscc 68 -54 PIPS 58 19 -(c) Merge	⊇ Aspic iscc PIPS ⊇ Aspic iscc PIPS	Aspic Aspic iscc PIPS iscc -21 27 54 -48 51 36 -successes: 83 (b) Split failures: 19 2 2 29 40 Aspic iscc PIPS -21 22 69 -58 68 33 -1 (d) Merge-split	8	1	PIPS
		(a) Direct Table 9.2: Invariant inclusions	(b) Split
	successes: 73 failures: 29		Aspic 1			successes: 88 failures: 14	Aspic 1
		21	36	1			20	47	2
	iscc	11	2	1	PIPS	iscc	4	12		2	PIPS
		(c) Merge			(d) Merge-split
		Figure 9.4: Venn diagrams of successes for each tool
	Globally, the merge-split restructuring leads to the best results with 88 out of 102 test cases correctly solved.

Pour éviter des propagations infinies, des opérateurs spéciaux d'approximation, comme l'élargissement, sont utilisés pour garantir la convergence en un nombre fini d'étapes.

Au lieu de cela, PIPS s'appuie sur une approche alternative, modulaire, en utilisant des transformeurs affines

[START_REF] Irigoin | Semantical Interprocedural Parallelization: An Overview of the PIPS Project[END_REF]

, c'est-à-dire des polyèdres représentant des fonctions de transfert. Dans la section 10.1, nous présentons comment PIPS calcule les transformeurs puis les invariants. La section 10.2 donne un aperçu de la modularité et de la précision de l'approche par transformeurs. Dans les sections

10.3 et 10.4

, plusieurs améliorations sont apportées aux algorithmes de calcul des transformeurs afin d'augmenter la précision des invariants ; ces améliorations sont évaluées expérimentalement dans la section

10.5.

Table 10 .

 10 For comparison purposes, we also use another polyhedral analyzer for C programs, PAGAI †1 , by Henry et al.[START_REF] Henry | PAGAI: A Path Sensitive Static Analyser[END_REF][START_REF] Henry | Succinct Representations for Abstract Interpretation[END_REF]. PAGAI -Path Analysis for invariant Generation by Abstract Interpretation -is a static analyzer working over the LLVM compiler infrastructure, which computes inductive invariants on the numerical variables of the analyzed program. PAGAI implements various state-of-the-art algorithms combining abstract interpretation and decision procedures (SMT-solving), focusing on distinction of paths inside the control flow graph while avoiding systematic exponential enumerations.Since the tools have very different structures and execution times, we report either directly the sum of the User and IO times reported by the time command for Aspic, iscc and PAGAI, or the times obtained using LOG_TIMINGS for the sole transformer and precondition passes of PIPS. In this way, the C parsing part of PIPS is eliminated, as it is for Aspic and iscc who use internal formats and for PAGAI who uses Clang for parsing, and we compare the execution time of the 2: Times in seconds for interprocedural analyses and intraprocedural analyses of inlined versions

	main-1 2 calls inlined 3 calls inlined 4 calls inlined 5 calls inlined 6 calls inlined main-2 main-3 main-4 main-5 -0.043 -0.061 -0.087 -0.108 -0.149 -261.810 -274.580 -370.960 -413.300 -456.360 PAGAI 0.980 Aspic iscc 1.417 1.383 5.680 2.030 14.677 2.990 30.007 4.467 53.247 PIPS 0.048 0.043 0.049 0.063 0.048 0.084 0.050 0.108 0.051 0.127

Table 10 .

 10 3 provides the numbers of successes and the total execution times obtained with the four different encodings of ALICe benchmark (Chapter 9) by Aspic, iscc, the legacy version of PIPS and PIPS with some of the improvements defined in Section 10.3. Accuracy figures for

	fsm Direct 5497 Split 9348 Merged 5579 Merged 6206 & Split	C Raw 15482 Simplified 4323 Raw 757222 Simplified 4199 Raw 15380 Simplified 5187 Raw 38941 Simplified 4549

Table 10 .

 10 4: Size of fsm and C files in number of lines Chapters 8 and 9. In previous chapters, the execution time is:

  Les systèmes de contrôle sont largement utilisés dans l'industrie, et interviennent dans un grande variété d'applications. Un sujet important en théorie du contrôle est de garantir la stabilité d'un système de contrôle, c'est-à-dire que le système démarré près d'un point d'équilibre restera dans un voisinage de ce point tout au long de l'exécution. La théorie du contrôle énonce qu'un système de contrôle linéaire est stable s'il existe un invariant caractéristique, appelé invariant de Lyapunov, et de nombreux systèmes sont conçus et prouvés en utilisant ce théorème. Mais les preuves de stabilité de ces systèmes reposent sur l'arithmétique des nombres réels, et ne s'appliquent pas nécessairement lorsque le système est implémenté sur une architecture avec une arithmétique machine.Nous introduisons un nouveau framework pour traduire les invariants de preuves en nombres réels vers des invariants en nombres machine, en prenant en compte la représentation des nombres et les problèmes d'arrondis. Les invariants traduits sont exacts, mais ne garantissent pas nécessairement la stabilité du système. Le schéma de traduction générique a été vérifié à l'aide de Coq, et est implémenté pour l'arithmétique à virgule flottante, la plus répandue des arithmétiques machine. Cette implémentation, LyaFloat, peut vérifier automatiquement la stabilité d'un système sur des architectures à virgule flottante, et calculer le nombre de bits nécessaires à la fois pour la stabilité en boucle ouverte et en boucle fermée. discrétisation d'un contrôleur continu, ce qui pose le problème de garantir des bornes de temps de calcul dans la boucle de contrôle, ou l'acquisition des données, qui est implémentée sur la plupart des contrôleurs sous la forme de routines de traitement d'interruptions.Dans ALICe, la manipulation de code source en C s'est avérée problématique. Actuellement, tous les cas de test sont implémentés sous forme de machines à états finis, dans le format fsm, et sont compilés en C avant d'être analysés par PIPS. La transformation de fsm en C n'est pas simple et les codes C générés tendent à être excessivement longs et complexes et donc difficiles à analyser, ce qui désavantage PIPS en comparaison avec Aspic et iscc. Même si le problème peut être en partie contourné avec une phase intermédiaire de simplification, nous pensons qu'il faudrait employer une meilleure stratégie pour les codes C, que ce soit en enregistrant une version C de chaque modèle en plus de la version fsm (ce qui amène d'autres problème : choix des représentations, expressivité, cohérence, encodage), ou bien en améliorant considérablement l'outil fsm2c afin de résoudre ces problèmes.Le domaine d'application d'ALICe pourrait également être étendu, d'une part en ajoutant de nouveaux modèles à la base de cas de test, et d'autre part en supportant davantage d'outils d'analyse. Cette dernière tâche demande du temps puisqu'à cause du manque de standardisation des formats utilisés par les outils d'analyse, il faut écrire des routines de conversion de formats pour chaque outil ajouté. Il faudrait aussi vérifier que les outils d'analyse génèrent des invariants réalistes, par exemple en associant à chque cas de test un invariant minimal. Un outil d'analyse ne serait considéré correct que si les invariants qu'il génère sont des sur-ensembles des invariants minimaux. Cela éviterait qu'un outil trivialement faux puisse bouleverser les résultats.Finalement, même si les améliorations apportées au calcul d'invariants par transformeurs ont permis d'améliorer significativement les résultats dans PIPS, les coûts de calcul les rendent impraticables pour l'analyse de gros programmes. Des algorithmes précis et de complexité raisonnable sont encore à trouver.

	Analyse des relations linéaires

send(u, 1); % x c ∈ E R

†1 La procédure qui suit serait exactement la même avec des nombres flottants 32 bits, les résultats numériques seraient simplement différents. La précision numérique est en fait un paramètre du système.†2 The procedure that follows would be exactly the same with 32-bit floating-point numbers, only with different numerical results. In fact, arithmetic precision is a parameter of the system.

†1 In the following, a straight line will be confounded with one of its direction vectors. Likewise, a ray will be confounded by one of its direction vectors.

†4 It would not be much more difficult to consider more general affine relations of the form F X + F ′ X ′ ≤ G.

†8 See http://gmplib.org/. †9 See http://www.cs.unipr.it/ppl/. †10 See http://pips4u.org/.

†1 See http://pagai.forge.imag.fr/.

Remerciements

, qui a été mon tuteur lors de mes études à

List of Figures

List of Tables

isl

The Integer Set Library †11 (or isl) [START_REF] Verdoolaege | isl: An Integer Set Library for the Polyhedral Model[END_REF] is a thread-safe C library for manipulating sets and relations of integer points bounded by affine constraints:

The descriptions of the sets and relations may involve both parameters (s in Equation (7.1)) and existentially quantified variables (z in Equation (7.1)). All computations are performed in exact integer arithmetic using GMP. The isl library offers functionalities that are similar to that offered by the Omega+ library †12 , but the underlying algorithms are in most cases completely different. 4. Use iscc to check whether the intersection of the model error region and the invariant computed by the analyzer is empty, i.e. whether the analyzer is able to solve the test case.

Those steps are illustrated in Figure 8. 4.

Translation Tools

Implementing these steps required to develop several translation programs, from and to Aspic, isl and PIPS formats. In particular, the tools fsm2c and fsm2iscc, presented below, respectively

Generation of Transformers by PIPS

The algorithms used in PIPS assume no cycles in the call graph and proceeds as follows.

First, each program command S, elementary or compound statement or procedure call, is over-approximated by an affine transformer T (S, P ), possibly using information about a precondition P of S. This is a bottom-up procedure, detailed in the next paragraphs, because a default value can be used for the precondition P when no information is available. Each function is analyzed once and its transformer is reused at each call site. Then, preconditions are propagated from the program starting point using the transformers.

This approach can also be used with unstructured programs: the control flow graphs are either turned into equivalent structured graphs [START_REF] Ammarguellat | A Control-Flow Normalization Algorithm and its Complexity[END_REF], or approximated and simplified by adding transitions, or analyzed using a decomposition into cycles [START_REF] Bourdoncle | Abstract Interpretation by Dynamic Partitioning[END_REF].

void foo(float x) { // P() {} int n = 0; // P(n) {n==0} while (1)

3: Preconditions for counter

The two-stage approach used by PIPS is the following. We suppose that elementary instructions have been turned into transformers, and show how control structures are handled.

Sequence

Let x i , x ′ i and x ′′ i denote values of variable x i at different states. A sequence of affine transformers T 1 followed by T 2 is overapproximated by the union of constraints in T 1 (on values

Choice

The effect of a choice T 1 or T 2 is the transformer T 1 ∪T 2 , which is not affine in the general case (the union of two convex polyhedra is not a convex polyhedron). The best convex approximation is the convex union T 1 ⊔ T 2 . This is a lossy operation in general.

Loop

Given an affine transformer T for a loop body, the affine transformer T * represents the effect of any number of iterations of T . It is computed by the Affine Derivative Closure algorithm [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF].

As noticed in Section 9.1, the two main sources of imprecision that appear are the abstractions of loops ( * ) and parallel paths (⊔). Their impact is cumulated when multiple control paths appear within loops and nested loops.

These imprecisions are accumulated if there are structures in the automaton with two or more loops involved in the same control. The state splitting heuristic described in Section 9.1 was developed with regards to this issue. Two alternative possibilities are explored in this chapter: transformer lists and iterative analysis.

Note that the pure bottom-up approach may be used or not. Since transformers are not computed concurrently but by traversing the AST, information gathered previously can be used right away. The range of a transformer or the condition of a test can be used as a passes of PIPS that might be replaced by new passes based upon the other tools. Furthermore, the evolution of the execution times for each tool is fully relevant and interesting.

Benchmarks were run on a computer powered by a quad-core Intel i7-2600 processor clocked at 3.40 GHz with 16 GB of memory, using Aspic version 3.3, iscc version 0.12.2, PAGAI version 14-04-07 and PIPS revision 22 241.

Impact of Cycle Nesting on Convergence

Nested loops are common in scientific codes and easy to analyze with transformers. However, a 2D loop nest is used by Halbwachs in [START_REF] Halbwachs | When the Decreasing Sequence Fails[END_REF] to show improvements in widening techniques. Let us consider the code in Listing 10.4, a matrix multiplication. To check that all array accesses are safe, invariants on i, j and k are needed, which requires the analysis of a 3D loop nest. for (i = 0; i <= n-1; i += 1) for( j = 0; j <= n-1; j += 1)

for (i = 0; i <= n-1; i += 1) for (j = 0; j <= n-1; j += 1)

Listing 10.6: Matrix exponentiation, inlined with PIPS twice as fast as the depth and because the polyhedral operators are known for their exponential worst case complexities. However, it does not increase very fast, especially for the nesting depths usual in functions when inlining is not used.

Interprocedural Analysis or Inlining

The code in Listing 10.5, a matrix exponentiation, mp, contains a call to a matrix multiplication, mm, within a loop. Function mp is called from a main function that reads a matrix and an exponent, and prints the resulting matrix. This code can be analyzed interprocedurally or intraprocedurally after inlining, as show in Listing 10.6. Table 10.2 contains execution time measurements for main programs calling mp from one to five times. Each measurement was performed 10 times and the median time is displayed. The code is either analyzed interprocedurally or intraprocedurally after the callees are inlined, which reduces the number of call sites to zero. Modularity becomes more and more useful when multiple call sites of the same function are present in the analyzed program. Also, inlining increases the loop nest depths, which is bad for the execution time as seen in the section above.

tion function T is sharpened and the domain of the resulting transformer also is restricted by the previous precondition, P * n . The iterative refinement process does not always converge and it may even lead to a precision loss, due to magnitude overflows. These can be addressed with arbitrary precision numbers (Section 10. 3.2).

printf("property verified\n"); } else { printf("property not found\n"); } } Listing 10.9: Periodic behavior As pointed out in [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF], there are different ways to encode the swap, and the above invariant may be more or less easy to generate. However, regardless of the encoding, the behavior is always periodic. More information is preserved if the transitive closure of the loop transformer is computed as a function of one of its powers [START_REF] Ancourt | A Modular Static Analysis Approach to Affine Loop Invariants Detection[END_REF]. For instance, the square is useful for idempotent functions and functions with a period of 2:

This can be generalized to any power of T . As no application has yet required a larger periodicity, our current implementation in PIPS is limited to As expected, the execution times increase and the benefit is sometimes quite small for the test cases in ALICe. However, using sets of control path transformers and a favorable encoding provides excellent results in the number of test cases solved per second, and future benchmarks for invariant generation may have different profiles. Note that periodic analysis (Section 10. 4.2), not shown in Table 10.3, does not improve any of the test cases currently in ALICe.

Sizes in number of lines are indicated in Table 10. 4 for FAST files and for C files generated by fsm2c. C file sizes are reduced by eliminating files greater than 1 KLOC and by a first PIPS pass (control simplification).

Note also that the execution times of PIPS in Table 10.

Perspectives

Systèmes linéaires

Concernant les perspectives, plusieurs améliorations peuvent être apportées dans l'outil de traduction d'invariants LyaFloat. Actuellement, seule l'arithmétique à virgule flottante, la plus répandue, est supportée ; il serait relativement facile de supporter l'arithmétique à virgule fixe, au moins en théorie (l'implémentation est rendue difficile par l'absence de bibliothèque de calcul en virgule-fixe et précision paramétrique en Python). Dans la même veine, il serait intéressant de supporter davantage de fonctions et d'arguments de propagation, afin d'étendre la portée de l'outil.

Dans sa forme actuelle, LyaFloat génère des invariants pour du code en arithmétique machine en suivant un modèle formellement prouvé, mais il n'y a pas de garantie formelle que LyaFloat lui-même n'est pas défectueux. On gagnerait un niveau d'assurance s'il était possible de vérifier par un autre moyen les invariants générés. Une solution serait de réécrire et de formaliser LyaFloat en Coq plutôt que d'utiliser Python. Cela nécessiterait de prouver (ou au moins de formaliser) les arguments de propagation d'invariants requis dans les preuves de stabilité, ce qui représenterait un travail conséquent. Une autre idée serait de générer avec les invariants des scripts Coq qui feraient office de certificat de correction pour la propagation d'invariants, ce qui serait suffisant pour vérifier facilement les invariants générés. Cela nécessiterait de développer une bibliothèque Coq de théorèmes autour de la théorie de la stabilité.

D'un point de vue général, de nombreuses étapes sont nécessaires pour traduire un modèle de contrôleur linéaire idéal et prouvé vers le code C de bas niveau qui l'implémente, et changer la représentation arithmétique n'est que l'une de ces étapes. Parmi les autres étapes, on trouve la

List of Tables