domain: spi phys math info

text:
Professor Tsiotras Holzinger

Dr Hasnaa Zidani

Dr Cesar Munoz

Elaud, Louis Helene Gabriel

Phillipe Trey

Thunde

Mark Kevin Corbin

Matthieu, Emanuele, Kevin Hugo Guillaume

The work presented in this thesis is at the border of multiple different areas, including optimal control, optimization theory and formal methods. Therefore having such experts on Control Systems, Optimal Control, Space Systems and Computer Science on my committee is a honor and I fully realize the chance I have. My advisor, Dr. Eric Feron, thank you for believing in me and for offering me to work with you as your PhD student. Working next to you brought me a lot, both on a

LIST OF SYMBOLS Symbol Description

R

The set of all real numbers.

R +

The set of all positive real numbers.

R *

+

The set of all stricly positive real numbers.

R n

The set of real vectors of length n.

R m×n

The set of real matrices of size m × n.

A F Frobenius norm of a matrix A.

A Two norm of a matrix A.

x Two norm of a vector x.

B n n-dimensionnal unit Eucliean ball. B n = {z ∈ R n : z ≤ 1}.

B r (x) Ball of radius r centered on x. B r (x) = {z ∈ R n : z -x ≤ r}.

Ell(B, c) Ellipsoid set defined by: Ell(B, c) = {Bu + c : u ∈ B n }.

Vol() Volume of a given set.

fl() Floating-point rounding to nearest of a given real number.

σ max (A) Largest singular value of a matrix A.

σ min (A) Smallest singular value of a matrix A.

k(A)

Condition number of a matrix A.

x Plant State Vector.

N

Model Predictive Control Horizon.

x Collection of state vectors to horizon: x = [x 1 . . . Constants Introduces the section where constants can be defined.

Variables

Defines the decision variable for the optimization problem.

Input

Optional Section defining the MPC input.

Output

Optional Section defining the MPC output.

Minimize

Introduces the cost to be minimized.

SubjectTo Introduces the constraints of the optimization problem.

Information Introduces the scalars needed for the verification.

Acronyms

SUMMARY

The efficiency of modern optimization methods, coupled with increasing computational resources, has led to the possibility of real-time optimization algorithms acting in safety critical roles. However, this cannot happen without addressing proper attention to the soundness of these algorithms.

This PhD thesis discusses the formal verification of convex optimization algorithms with a particular emphasis on receding-horizon controllers. Additionally, we demonstrate how theoretical proofs of real-time optimization algorithms can be used to describe functional properties at the code level, thereby making it accessible for the formal methods community. In seeking zero-bug software, we use the Credible Autocoding scheme. In this framework, with the use of a "Credible Autocoder", we are able to automatically generate C code implementation of Receding Horizon controllers along with its proof of soundness at code level.

We focused our attention on the Ellipsoid Method solving Second-Order Cone Programs (SOCP). Also, we present a modified version of the original Ellipsoid Method, in order to take into account and control its numerical error. Following this, a floatingpoint analysis of the algorithm and a framework is presented to numerically validate the method.

xvii Chapter I

INTRODUCTION

Cyber-physical systems (CPS) regroup all the mechanical systems that interact with computer-based algorithms. A cyber-physical system can be controlled by a human being or completely autonomous. Among those systems, we call safety-critical the ones for which a failure could cause human death. Such systems include planes, trains, nuclear power plants, human spaceflight vehicles, robotic surgery machines and many more.

Designing the embedded software carried by a safety-critical cyber-physical systems is a meticulous task and should be performed with caution. In most applications, safety-critical CPS are required to follow a given certification and requirements. In the case of Airborne Systems, the DO-178C, Software Considerations in Airborne Systems and Equipment Certification, give technical guidelines and requirements for developing avionics software systems. We discuss in more detail how formal methods fit into DO-178C and the scope of this PhD thesis in section 1.3.3.

The efficiency of modern optimization methods, coupled with increasing computational resources, has led to the possibility of real-time optimization algorithms acting in safety-critical roles. However, this cannot happen without addressing proper attention to the soundness of these algorithms. This PhD thesis discusses the formal verification of convex optimization algorithms with a particular emphasis on recedinghorizon controllers. We demonstrate how theoretical proofs of real-time optimization algorithms can be used to describe functional properties at the code level, thereby making it accessible for the formal methods community.

The need for enhanced safety and better performance is currently pushing for the introduction of advanced numerical methods into next generations of cyber-physical systems. While most of the algorithms described in this paper have been established

for a long time, their online use within embedded systems is relatively new and opens issues that have to be addressed. Among these methods, we are concerned with numerical optimization algorithms.

In this chapter, we recall preliminaries about formal verification, semantics of programs, and the Credible Autocoding framework. Finally, we develop the state of the art and explain how this thesis fits into the previous work. In Chapter 2, we give background about convex optimization, dynamical systems and model predictive control (MPC). Chapter 3 focuses on the axiomatization of a second-order cone program and the formal verification of the ellipsoid algorithm. A floating-point analysis of a modified version of the ellipsoid method is presented in Chapter 4, while Chapter 5 gives details about closed-loop management. This chapter presents how this framework can be automated and applied to a system as well. Chapter 6 concludes this thesis.

Credible Autocoding Framework

Credible autocoding, is a process by which an implementation of a certain input model in a given programming language is being generated along with formally verifiable evidence that the output source code correctly implements the input model.

Given that the mathematical proofs of high-level functional properties of convex optimization algorithms do exist, we want to translate, generate and carry them at code level. This is done by adding comments, which does not perturb the code compilation and execution. Furthermore, those comments are expressed in a formal specification language in order to be read and analyzed by other software. An illustration of this framework is given in Figure 1. In this thesis, we focus on automatically generating 3). In other word, we instantiate a given algorithm on a model predictive control (MPC) problem. This framework is similar to [START_REF] Mattingley | Cvxgen: A code generator for embedded convex optimization[END_REF].

For that, we build an Autocoder (see Chapter 5) that will automatically generate C code implementation of Receding Horizon Controllers from a text file containing the high-level MPC formulation. Along with generating the C code, this autocoder, that we call Credible Autocoder, also generates the semantics of the corresponding algorithms and its proof of soundness at code level. From a MPC formulation specified using a given programming language (detailed in Section 5.3) written by the user, the Credible Autocoder will generate C code algorithms implementing the original MPC formulation. As it was said before, along with generating those algorithms, it also generates the proof of soundness of those algorithms at code level. Once the semantics annotated C code generated, it is checked using a software analyzer and if it is proven correct, compiled then embedded at the heart of a system's feedback. Semantics of programs express their behavior. It gives a rigorous mathematical description of the meaning of a given program. For the same program, different means can be used to specify its semantics:

• a denotational semantics, expressing the program as a mathematical function,

• an operational semantics, expressing it as a sequence of basic computations, or

• an axiomatic semantics, as a set of observations.

In the latter case, the semantics can be defined in an incomplete way, as a set of projective statements, ie. observations. This idea was formalized by [START_REF] Floyd | Assigning meanings to programs[END_REF] and then [START_REF] Hoare | An axiomatic basis for computer programming[END_REF] as a way to specify the expected behavior of a program through pre-and post-condition, or assume-guarantee contracts.

Hoare Logic

A piece of code C is axiomatically described by a pair of formulas (P, Q) (see Fig- ures 3) that if P holds before executing C, then Q should be valid after its execution.

This pair acts as a contract for the function and (P, C, Q) is called a Hoare triple.

In most uses P and Q are expressed as first order formulas over the variables of the program. Depending on the level of precision of these annotations, the behavior can be fully or partially specified. In our case we are interested in specifying, at code level, algorithm specific properties such as the convergence of the analysis or preservation of feasibility for intermediate iterates. Software frameworks, such as the Frama-C plat- form [START_REF] Cuoq | Frama-c: a software analysis perspective[END_REF], provide means to annotate a source code with these contracts, and tools to reason about these formal specifications. For the C language, ACSL [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language. version 1.11[END_REF], (ANSI C Specification language) can be used as source comments to specify function contracts, or local annotations such as loop invariants. Local statements annotations act as cuts in proofs and are typically required when analyzing loops. Figure 4 shows an example of a Hoare triple. In this example, if the triple is indeed valid, we can conclude that assuming 0 ≤ x ≤ N is true, then after executing the command x := x + 1 we know that the property 1 ≤ x ≤ N + 1 will be true. Figure 5 displays another example of a Hoare triple. The two examples presented show correct Hoare triples and have been proven using Frama-C. ACSL also give means to introduce function contract. Those contracts consist in listing all the properties that we assume to be true before the execution of the function and all the properties that will be true after its execution.

For this, we use the keywords requires and ensures. For every function contract, we will then try to prove that, after the execution of a function, all the properties ACSL + C

1 // @ assert 0 <= x <= N ; 2 x = x + 1; 3 // @ assert 1 <= x <= N + 1; 1 // @ assert -2 <= x <= 2; 2 y = x * x ; 3 // @ assert 0 <= y <= 4; listed following a ensures keyword will be true assuming that all the properties listed within a requires were true before its execution. The part of the memory assigned by a function can also be specified using the keyword assigns. Checking the memory assignment is very important and could lead to checking errors if it is omitted.

Let us take a look at a simple and specific example. In Figure 6, we show an implementation of a function, "addInt", that takes two integers as inputs and store the addition of those two integers at a given memory address, referenced as an input pointer. Also, this function returns the corresponding addition. Let us assume the user only specifies that the output corresponds to the sum of the input integers. Then, while proving the second property at line 16, the SMT solvers would assume the value of a did not change and is correct. Therefore, missing an assigns clause could lead to logic issues and checking errors. Although Frama-C rises a flag reminding the user that a assign clause is missing, it still lacks the level of confidence we would like to achieve. In order to get around this, we make sure we write an assigns clause for every function contract.

Let us say that the user did write an assigns clause but a non correct one. Could this lead to logical issues? The answer is no. Indeed, the solvers might proved correctness for properties that are not correct but they will not be able to prove correctness for

7 ACSL + C 1 # include < stdio .h >

Weakest Precondition

In this section we give details about how this process can be automated at code level.

Meaning that, how can the correctness of a Hoare triple be checked automatically.

For this, we present how to prove a Hoare triple by weakest precondition. Let Q be a assertion and C a command. We define the weakest precondition of the couple (C, Q) as the weakest assertion P such that the triple {P } C {Q} is valid. For two assertions P 1 and P 2 , P 1 is weaker than P 2 means that P 2 =⇒ P 1 . The advantage of proving Hoare triple by weakest precondition lies in the fact that this latter step is easily automated. Indeed, the weakest precondition of a couple (C, Q), depending of the nature of the command C, follows a simple rule. Once this weakest precondition online with great success for the guidance of systems within safety-critical applications, including, autonomous cars [START_REF] Falcone | Predictive active steering control for autonomous vehicle systems[END_REF][START_REF] Raffo | A predictive controller for autonomous vehicle path tracking[END_REF][START_REF] Jerez | Embedded online optimization for model predictive control at megahertz rates[END_REF] and reusable rockets [START_REF] Ac ¸ikmese | Lossless convexification of nonconvex control bound and pointing constraints of the soft landing optimal control problem[END_REF][START_REF] Blackmore | Lossless convexification of control constraints for a class of nonlinear optimal control problems[END_REF]. The latter case has resulted in spectacular experiments, including landings of SpaceX's Falcon 9 and BlueOrigin's New Shepard. In order to perform the landing of the Falcon 9, the spaceX's teams used convex optimization technique [START_REF] Blackmore | Autonomous precision landing of space rockets[END_REF]. The landing problem was formulated as a Quadratic Programming (see Section 2) and ran during the descent phase of the Falcon 9 rocket. In addition, automatic code generation has been used for numerical optimization, with the use of the code generator for convex optimization CVXGEN [START_REF] Mattingley | Cvxgen: A code generator for embedded convex optimization[END_REF]. Thus, powerful algorithms solving optimization problems are already used online, have been embedded on board, and yet still lack the level of qualification required by civil aircraft or manned rocket flight. Automatic code generation for solving convex optimization problems has already been done [START_REF] Mattingley | Cvxgen: A code generator for embedded convex optimization[END_REF][START_REF] Mattingley | Code generation for receding horizon control[END_REF]. In these articles, the authors present a new code generator for embedded convex optimization. The generated code is said to be library-free, auto sufficient and without failure. Those algorithms are very complex and arguing that they do not contains any bug could be delicate, especially for safety-critical conditions. Also, this work does not include the use of formal methods.

Likewise, work within the field of model predictive control already exists where numerical properties of algorithms are being evaluated [START_REF] Patrinos | A dual gradientprojection algorithm for model predictive control in fixed-point arithmetic[END_REF]. This paper, [START_REF] Patrinos | A dual gradientprojection algorithm for model predictive control in fixed-point arithmetic[END_REF], proposes a dual gradient projection (DGP) algorithm implemented on fixed-point hardware.

A convergence analysis with taking into account round-off errors due to fixed-point arithmetic is presented. Following this, the minimum number of integer bits that guarantee convergence to a solution can be computed, minimizing the computational power. Nevertheless, this work is only valid for Quadratic Programming and using fixed-point numbers. Additionally, no formal verification was performed.

Similarly, as it is presented at [START_REF] Jerez | Embedded online optimization for model predictive control at megahertz rates[END_REF], MPC controllers have been successively implemented on FPGAs, letting the possibility to ran some problems at megahertz rate.

In this article, the authors focus on the alternating direction method of multipliers (ADMM) and perform a numerical analysis in fixed-pint arithmetic.

Formal Verification of Control Systems Software

Since the 60s, computer scientists were studying ways of analyzing programs. For this, different techniques were proposed based on mathematical supports. Those methods, developed for analyzing software, were based on mathematical descriptions of the programs behaviors (i.e. its semantics): a formal description. formal methods represents an advantage over testing because they guarantee the correction of properties and soundness of programs for a wide range of inputs and under certain known hypothesis. When testing programs, only certain test cases are executed, giving satisfactions when no run-time errors are detected. Unfortunately, the absence of run-time errors does not imply the correction of the program for other test cases, as close as they can be.

We shortly present in the following paragraphs, the different formal methods developed and successfully used to analyze embedded safety-critical software.

Abstract interpretation is one the most successful method developed towards the static analysis of programs. It was first proposed in the 70s. In practice, it is mainly used in order to compute numerical invariants over programs. For this, abstract interpretation uses results from set theory to compute over approximations of program behaviors. Abstract domains are used and depending on the programs nature, different shapes and geometries are considered, representing trade offs between accuracy and performance. Definitions and the use of the main abstract domains are detailed in [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF][START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Miné | A new numerical abstract domain based on difference-bound matrices[END_REF][START_REF] Miné | The octagon abstract domain[END_REF]. In 2010, a remarkable advanced was made towards formal verification of safety-critical control system using abstract interpretation by successfully proving the absence of runtime error of the flight control system of the Airbus A380 [START_REF] Kästner | Astrée: Proving the absence of runtime errors[END_REF].

On the other hand, another technique, called Satisfiability Modulo Theories (SMT), is growing in the area of formal verification for control system. SMT solvers are decision problems for logical formulas with respect to combinations of background theories expressed in first-order logic such as linear real/integer arithmetic. Roughly speaking, those SMT solvers are deciding procedures for the satisfiability of conjunctions of items, where an item is an atomic formula. SMT solvers handle sub-formulas by performing case analysis, which is the technique used in most automated deduction tools. SMT-solvers are used as back-end reasoning engines in a wide range of formal verification applications, such as model checking, test case generation, static analysis, etc. In this PhD thesis, we focus on the SMT solver Alt-Ergo connected to the software analyzer back-end Frama-C.

Contributions for higher-level properties have been made concerning formal verification of control systems [START_REF] Feron | From control systems to control software[END_REF][START_REF] Champion | Formal methods for the analysis of critical control systems models: Combining non-linear and linear analyses[END_REF][START_REF] Herencia-Zapana | Pvs linear algebra libraries for verification of control software algorithms in c/acsl[END_REF][START_REF] Roux | A generic ellipsoid abstract domain for linear time invariant systems[END_REF]. Those articles mainly focus on formal verification and code generation for linear control system and typical feedback control techniques. The authors show how theoretical system level properties (closed-loop Lyapunov stability) can be carried out from a high level programming language like Simulink all the way to code through automatic code generation. Following this, research has also been made toward the verification of numerical optimization algorithms [START_REF] Wang | Credible autocoding of convex optimization algorithms[END_REF], yet it remains purely theoretical and no proof was actually performed with the use of formal methods.

Formal Verification And Software Certification

In this section, we explain the differences between formal verification, testing and certification. The same way it is for Cyber-Physical System, we qualify a software of safety-critical if its failure can lead to catastrophic consequences. Therefore the correctness is a crucial issue in the design process of safety-critical software. For those systems, a higher authority is usually present and act as a supervisor checking software quality. For airborne systems many authority exists such as the FAA or EASA.

Certification is the legal recognition by a certification authority that a product fulfills specific requirements.

Software verification regroups all the computer science techniques that aim to show that a piece of software complies with predefined properties and that does not perform unexpected behavior. Until recently, the high majority of software verification performed by companies while developing software consisted of reviewing and testing.

As it is recalled in [START_REF] Gigante | Formal methods in avionic software certification: The do-178c perspective[END_REF], 30% to 50% of the software costs are dedicated to testing.

Unfortunately, detecting bugs by reviewing and testing is very limited, especially for safety-critical applications. Testing can only cover certain scenarios, and bugs can still be present even when a software passes all the test cases. Also, reviewing does not represent the perfect option because the reviewers might embrace the same logic as the programmer while reading the code and fall into the same mistakes. The reviewer's mind gets corrupted by adopting the programmer's point of view. For airborne software, the the legal authorities such as the FAA, EASA, follow the guidelines detailed in the document "DO-178C, Software Considerations in Airborne Systems and Equipment Certification" [START_REF] Rtca | DO-178C, software considerations in airborne systems and equipment certification[END_REF].

DO-178C Software Requirements: In the DO178C, it is explained that each part of the software is assigned a specific software level related to the criticality of the software. Software levels are assigned from levels A to E, A corresponding to the most critical piece of software. Thus, a level A software failure can lead to tremendous catastrophes. One on the requirement listed in the DO-178C requires that each line of code must be directly traced to a requirement. Also it specifies that this requirement needs to be traced to a corresponding test case.

Roughly speaking, the DO-178C decomposes the software design phase into three processes:

• the software planning,

• the software development,

• and the integral.

The software planning defines and manages all the software development-related activities. During the software development process, High-Level Requirements (HLR), software architecture and Low-Level Requirements (LLR) are defined. The software design, coding and integration are being performed during this process as well. The integral process addresses the verification. It aims to guarantee the correctness, control, and confidence in the safety and reliability of all previous lifecycle processes.

More details on DO-178B and DO178C can be found at [START_REF] Gigante | Formal methods in avionic software certification: The do-178c perspective[END_REF][START_REF]DO-178B: Software Considerations in Airborne Systems and Equipment Certification[END_REF][START_REF] Rtca | DO-178C, software considerations in airborne systems and equipment certification[END_REF] For each processes DO-178C gives requirements and guidelines on how formal methods can be incorporated. Also it specifies that formal methods cannot be used as a replacement for any existing requirements. The use of formal methods for verification purposes should be explained and fixed during the software planning process. All the assumption made during the formal modeling and verification should be specified.The properties aiming to be verified, the formal verification work plan and objectives should be detailed and listed beforehand. DO-178C specifies the repeatability characteristic that the formal verification should fulfill.

Contribution: This PhD thesis discusses the formal verification of convex optimization algorithms with a particular emphasis on receding-horizon controllers. We demonstrate how theoretical proofs of real-time optimization algorithms can be used to describe functional properties at the code level, thereby making it accessible for the formal methods community. This PhD thesis only dedicated work towards the generation of C code and checking its correctness. Indeed, the compiling process and its verification are not part of this PhD thesis scope. Work has previously been done towards those issues [START_REF] Leroy | Formal verification of a realistic compiler[END_REF]. In this PhD dissertation, we will present the following scientific contributions:

• the axiomatization of optimization problems using the specification language ACSL (mathematical properties of linear algebra, set theory and optimization theory were defined),

• the formalization of the algorithm's proof using this same language at code level, where ACSL annotations and Lemmas were formalized and proved,

• a modification of the original Ellipsoid Algorithm in order to account for numerical errors along with its complete numerical analysis.

• the development of a high-level optimization parser using a specific and predefined language (detailed in section 5).

• the development of an Autocoder (detailed in section 5) generating C code implementation of Convex Optimization Algorithms along with ACSL annotations [START_REF] Cuoq | Frama-c: a software analysis perspective[END_REF] Chapter II

OPTIMIZATION BASED CONTROL OF LTI SYSTEMS

Optimization based control deals with the problem of finding the control law for a given system such that a certain optimality criterion is achieved. In order to achieve this goal, two methods can be used: direct and indirect. In the latter case, we seek an analytical and closed form solution for the control input with the use of mathematical principles such as the Pontryagin maximum principle (PMP). The most famous result of this theory being the Riccatti equation, giving a closed form solution to a controller minimizing a quadratic cost on the inputs and states for a linear plant.

Further information about indirect methods can be found at [START_REF] Bonnans | Singular arcs in the generalized goddard's problem[END_REF][START_REF] Tsiotras | Goddard problem with constrained time of flight[END_REF]. In this thesis, we focus on direct methods, where we usually work with discretized systems and express the optimal behavior of a system by formulating an optimization problem. This optimization problem is then solved using numerical solvers and the optimal input sequence is found. In this chapter we recall first the setting of LTI theory and some stability results of LTI systems. Then, more details about convex programming and receding horizon controllers are given.

Linear Time Invariant (LTI) Systems

LTI Theory and State Space Realization

Given a system, a state-space representation is a mathematical model of a physical system as a set of inputs, outputs and state variables related by first-order differential equations.

In a general setting, a dynamical system can be represented by the system of equations 1. We write the state vector x, the input vector u and the output vector y.

       ẋ = f (x, u, t) y = h(x, u, t) (1)
Among dynamical systems, autonomous time-invariant systems are very important and represent a good trade-off between being simple enough to analyze but complicated enough to capture a wide variety of systems. For those systems, many stability results exist. In a general setting, an autonomous time-invariant system can be represented by the system of equations 2.

       ẋ = f (x) y = g(x) (2)
For LTI systems, the functions f and h are required to be linear in the input and state and to be independent of the time. Thus, LTI systems can be represented by the system of equations 3.

       ẋ = Ax + Bu y = Cx + Du (3)
Consequently, a LTI system that is additionally autonomous can be represented by equation 4.

ẋ = Ax (4)

Stability of LTI Systems

When analyzing dynamical systems or designing controllers, studying stability is a crucial characteristic. In most cases, we study the stability of a closed-loop system.

Meaning that the controller has already been designed and the loop closed (see Fig- equilibrium for the autonomous system 2 if ẋe = f (x e ) = 0. For autonomous LTI systems, the point x e = 0 is always an equilibrium. In term of stability, each equilibrium can be characterized as stable, asymptotically stable or unstable.

Definition 1 (Equilibrium Stability) The equilibrium point x e of the system 2 is:

• stable if, ∀ > 0, ∃ δ > 0 such that: x(0) -x e ≤ δ =⇒ x(t) -x e ≤ ∀t ≥ 0,
• unstable, if it is not stable,

• asymptotically stable if, We state now the Lyapunov Stability theorem which gives means to conclude on system stability under existence of a continuously differentiable function that we call a "Lyapunov function".

∀ > 0, ∃ δ > 0 such that: x(0) -x e ≤ δ =⇒ x(t) -x e ≤ ∀t ≥ 0,
Theorem 2 (Lyapunov Stability) Let x e be an equilibrium point for the system 2 and D ∈ R n be a domain containing x e . Let V : D → R be a continuously differentiable function such that:

V (x e ) = 0 , V (x) > 0 ∀x ∈ D\{x e },
and V (x) ≤ 0 ∀x ∈ D

Then x e is stable. Moreover, if

V (x) < 0 ∀x ∈ D\{x e },
Then x e is asymptotically stable.

For autonomous LTI systems, the results we have are even stronger than this last theorem. We know that we only have to look for quadratic Lyapunov functions.

Furthermore, we know that the solution of a LMI could give us a proof of stability.

Theorem 3 Suppose there exists a symmetric matrix P = P T ∈ R n×n such that P > 0

A T P + P A < 0
Then the system 4 is asymptotically stable.

In this section, we recalled definitions and results known about a subclass of dynamical systems. As it was said earlier, when designing a controller and closing the loop, studying the stability of the resulting system, is the first and most crucial job of the controls engineer. Further details and proof about control feedback system can be found at [START_REF] Ogata | System dynamics[END_REF][START_REF] Brogan | Modern control theory[END_REF]. In our case, because we focus on model predictive control (MPC), the controller carries an optimization algorithm and is therefore a complicated entity.

Nevertheless, having closed-loop stability guarantees cannot be avoided and details regarding this manner are discussed in Section 2.3.1. Additionally, we clarify the technical differences between MPC and Path-planning and explain why confusing those two techniques could be an issue for our application.

Convex Optimization Problems

Optimization algorithms solve a constrained optimization problem, defined by an objective function, the cost function, and a set of constraints to be satisfied:

min f o (x) s.t. f i (x) ≤ b i for i ∈ [1, m] (5)
This problem searches for x ∈ R n , the optimization variable, minimizing

f o ∈ R n →
R, the objective function, while satisfying constraints f i ∈ R n → R, with associated bound b i . An element of R n is feasible when it satisfies all the constraints f i . An optimal point is defined by the element having the smallest cost value among all feasible points. An optimization algorithm computes an exact or approximated estimate of the optimal cost value, together with one or more feasible points achieving this value. A subclass of these problems can be efficiently solved: convex problems. In these cases, the functions f o and f i are required to be convex [START_REF] Boyd | Convex optimization[END_REF]. When optimization algorithms are used offline, the soundness of their implementation and the feasibility of the computed optimizers is not as critical and solutions could be validated a posteriori [START_REF] Roux | Closed loop analysis of control command software[END_REF].

Here, we only present a specific subset of convex optimization problems: Second-Order Cone Programs. For x ∈ R n , a SOCP in standard form can be written as:

min f T x s.t. A i x + b i 2 ≤ c T i x + d i for i ∈ [1 , m] With: f ∈ R n , A i ∈ R n i ×n , b i ∈ R n i , c i ∈ R n , d i ∈ R. (6)
The focus of this PhD thesis is the online use of Convex Optimization algorithms.

Model predictive control or real-time based control is more general and the convex aspect is not necessary in general. When optimization problems that are not convex are being solved, in the general case, there is no guarantee that the global minimizer Frequently, optimization problems that are used online for control systems can be formulated as a SOCP. Therefore, only dealing with SOCP is legitimately not constraining for the application we are targeting. In control systems, SDP's are mostly used off-line, a priory checking system's stability [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF].

Real-time Convex Optimization Based Control and MPC

Model Predictive Control (also known as receding horizon control) is an optimal control strategy based on numerical optimization. In this technique, a dynamical model of the plant is being used to predict potential future trajectories. Also, a cost function J, that depends on the potential future control input and state, is being considered over the receding prediction horizon N and the objective here is to minimize this cost J. At each time t, a convex optimization problem where J has to be minimized, is being solved. From the solution of this problem, an optimal trajectory starting at x(t) is being calculated and the control input sent to the plant corresponds to the first input of this optimal trajectory. A time step later, at t + ∆t, the exact same process occurs and is repeated until final time. As an example, we present in problem 7 an example of a MPC formulation and an illustration of this process in Figure 13.

minimize x,u N k=1 J(x k , u k ) subject to x k+1 = Ax k + Bu k , ∀k ∈ [1, N ] x k ∈ X , u k ∈ U, ∀k ∈ [1, N ] x 1 = x(t) (7)
At each iteration, we aim to solve problem 7. For a given step k, we write

û 1,k û 2,k . . . ûN-1,k
the optimal input sequence found by the algorithm. The resulting input sent to the plant is therefore ûk = û1 (Figure 12).

When using online convex optimization techniques for the guidance of system, we usually account for the plants dynamics by constraining the trajectories. It is done with the help of equality constraints and therefore, only linear dynamics can be captured. For this reason, throughout this PhD thesis, we consider only linear dynamics.

Spring-Mass System

As an example, we develop a MPC controller to control a spring-mass system, shown in Figure 14. The state vector x regroups the position z

x k+1 = Ax k + Bu k C ûk x k
               ẋ =     0 1 -1 0     x +     0 1     u y = 1 0 x (8)
We want to apply the MPC controller detailed in problem 9. The initial point of the trajectory being the point

x o = 2 -1 T . minimize x,u N k=1 Q • x k subject to x k+1 = Ax k + Bu k , ∀k ∈ [1, N -1] u k ≤ 5 ∀k ∈ [1, N -1] x 1 = x(t) (9)
After discretizing the system at T = 0.01 sec (10 Hz) and for a MPC horizon of N = 10 we performed a simulation of the closed-loop system. Figures 15,16 and 17 show the resulting position, velocity and input force. We used the constants below:

Q =    5 0 0 1    ; A =    1 0.01 -0.01 1    ; B =    0 0.01    .
The matrix Q represents the weight associated with each components of the state vector. The main inconvenient of using indirect methods is the fact that a closed-loop solution is mathematically complex to found, particularly when having constraints on the control and the state. When using Receding Horizon Control one can use arbitrary cost function and constraints in the MPC formulation. When the cost and the constraints are convex, convex optimization solvers are available and a solution (when there is one) can be found rapidly, letting the possibility of using this technique online.

Difference between MPC and Path Planning

Model predictive control is a feedback technique used to control systems based on numerical optimization algorithms. Therefore, those controllers are meant to be used When implementing a path-planning scheme, the stability of the closed-loop system lies in the stability of the lower-level controllers and the correctness of the nominal trajectory generated beforehand. For receding horizon controllers, system stability is not that obvious. Using optimization techniques online, the closed-loop trajectory is not necessary close to the future potential trajectories computed at each iterations.

Because of this, we need more powerful mathematical tools in order to prove stability.

As it was explained in [START_REF] Jadbabaie | Stabilizing receding horizon control of nonlinear systems: a control lyapunov function approach[END_REF], Lyapunov functions could be used to prove stability of MPC controllers. In that paper, the authors show how end-point penalty can be used as Lyapunov functions and therefore prove closed-loop stability.

28

Chapter III

THE ELLIPSOID METHOD AND ITS SEMANTICS

In this chapter, we recall the last advances and steps in optimization theory and show its consequences in science and engineering. After this, we give details about the algorithm of interest in this thesis.

A Brief Recent History of Optimization Theory

Optimization has been a on going research area of mathematics for centuries. Many famous mathematicians were interested in it including Lagrange, Euler, Newton and many others. More recently, a revolutionary discovery was made in 1947 by George Dantzig with the elaboration of the Simplex Method solving Linear Programs. This method was a huge step forward in the area of linear optimization and represented at the time, an extremely efficient method although having poor theoretical complexity (exponential).

In 1970 the mathematicians Shor, Judin, and Nemirovski published a highly important result, the discovery of an algorithm, the Ellipsoid Method, solving in polynomial time any convex programs. Couple of years latter, in 1979 soviet mathematician Leonid Khachiyan applied this last method to the more specific setting of Linear Programs, showing for the first time the polynomial solvability of Linear Programs. At this point, although polynomial time algorithm for LP was found using the Ellipsoid Method, the Simplex Method was still widely used because of better running time.

Indeed, theoretically very slow, the Simplex Method remains in practice extremely fast.

Following this, in 1984, Indian mathematician Narendra Karmarkar published the first interior-point method solving linear programs in polynomial time. This drastically increased the interest for interior-point methods to solve linear optimization problems, which was until there only applied for non-linear optimization. This discovery was quite a revolution, unifying linear and nonlinear optimization techniques.

Up until 1984, linear programs and non-linear programs where seen as two different problems that should be treated with different types of methods. With Karmarkar's first interior-point method solving LP in linear time, the community realized that linear and non-linear optimization could be treated the same way and shared interesting properties. This last algorithm being also efficient on both complexity and practice point of view.

Following the increasing interest for interior-point methods, mathematician Nemirowski, in the 90s came up with a marvelous result, the extension of interior-point methods to solve semi-definite programs (SDP). Indeed, Arkadi Nemirovski showed and explained how interior-point methods could be expanded to solve semi-definite optimization problems. This result is quite remarkable, allowing the community to expand the known results of Linear Programs to semi-definite optimization problems.

The Ellipsoid Method

As it was recalled in the introduction on this chapter, the Ellipsoid Method was first published by the mathematicians Shor, Judin, and Nemirovski in 1970. This algorithm was said to solve convex problems in polynomial time. One of the most famous consequence of this discovery is the elaboration of the Khachiyan method, which is named after its author Leonid Khachiyan and consist of an application of this method to linear programs.

Despite its relative efficiency with respect to interior point methods, the Ellipsoid Method benefits from concrete proof elements and could be considered a viable option for critical embedded systems where safety is more important than performance.

For the Ellipsoid Method, several algorithms that are mathematically equivalent exist.

All algorithms implement the same successive ellipsoids and present a mathematically equivalent update. Because there are several ways of encoding and defining an Ellipsoid, several algorithms exist. The point here being the fact that those algorithms are not numerically equivalent and we need to focus on the most accurate one. An

Ellipsoid can be defined by a positive-definite matrix P > 0. Equivalently, because for all positive-definite matrix P > 0, there exists a square-root matrix B such that B 2 = P , one can choose to implement the same algorithm but propagating the successive square root matrices B. The algorithm presented in [START_REF] Boyd | Linear controller design: limits of performance[END_REF] show a implementation of the Ellipsoid Method propagating the positive-definite matrices P and the one developed in [START_REF] Nemirovski | Introduction to Linear Optimization[END_REF] (The first one published) show an implementation propagating the matrices B. The same way it is for the Kalman Filter, we experienced that the algorithm implementing the square-root matrices is more numerically accurate. For that reason, we chose to focus ourselves on this implementation. This chapter presents a way to annotate a C code implementation of the Ellipsoid Method to ensure that the code implements the method and therefore shares its properties (convergence, soundness). Before recalling the main steps of the algorithm, the needed elements will be presented.

Ellipsoids in R n . An ellipsoid can be characterized as an affine transformation of an Euclidean Ball. Before defining an Ellipsoid set, we first recall the definition for an Euclidean ball.

Definition 2 (Euclidean balls) Let n ∈ N we denote B n the unit Euclidean ball in

B n = {z ∈ R n : z ≤ 1}
Also, Vol(B n ) denotes its volume. Also, we define B r (x) as the ball of radius r centered

on x i.e {z ∈ R n : z -x ≤ r} . Definition 3 (Ellipsoid Sets) Let c ∈ R n and B ∈ R n×n a non-singular matrix (det(B) = 0). The Ellipsoid Ell(B, c) is the set : Ell(B, c) = {Bu + c : u T u ≤ 1} (10)
Also, because the convergence of the method is a consequence of a volume decreasing property, we define below the volume of an Ellipsoid.

Definition 4 (Volume of Ellipsoids) Let Ell(B, c) be an ellipsoid set in R n . We denote by V ol(Ell(B, c)) its volume defined as :

Vol(Ell(B, c)) = |det(B)| • Vol(V n ) (11)
Algorithm. Let us now recall the main steps of the algorithm detailed in [START_REF] Bland | The ellipsoid method: A survey[END_REF][START_REF] Nemirovski | Introduction to Linear Optimization[END_REF][START_REF] Boyd | Linear controller design: limits of performance[END_REF].

In the following, we denote E k = Ell(B k , c k ), the ellipsoid computed by the algorithm at the k -th iteration.

Ellipsoid cut. We start the algorithm with an ellipsoid containing the feasible set X, and therefore the optimal point x * . We iterate by transforming the current ellipsoid E k into a smaller volume ellipsoid E k+1 that also contains x * . Given an ellipsoid E k of center c k , we find a hyperplane containing c k that cuts E k in half, such that one half is known not to contain x * . Finding such a hyperplane is called the oracle separation step, cf. [START_REF] Nemirovski | Introduction to Linear Optimization[END_REF]. In our SOCP setting, this cutting hyperplane is obtained by taking the gradient of either a violated constraint or the cost function. e, one can update the ellipsoid E k to its next iterate E k+1 according to equations [START_REF] Brogan | Modern control theory[END_REF] and [START_REF] Champion | Formal methods for the analysis of critical control systems models: Combining non-linear and linear analyses[END_REF]. In addition to that, we know an upper bound, γ, of the ratio of Vol(E k+1 )

to Vol(E k ) (see Property 1).

c k+1 = c k - 1 (n + 1) • B k p , (12)
B k+1 = n √ n 2 -1 B k + n n + 1 - n √ n 2 -1 (B k p)p T (13)
with:

p = B T k e e T B k B T k e . ( 14
)
Termination. The search points are the successive centers of the ellipsoids. Throughout the execution of the algorithm, we keep track of the best point so far, x. A point

x is better than a point y if it is feasible and have a smaller cost. When the program reaches the number of iterations needed, the best point so far, x, which is known to be feasible and -optimal, is returned by the algorithm. We state a volume related property, at the origin of the algorithm convergence, then state the main theorem of the method.

Property 1 [Reduction ratio.] Let k ≥ 0, by construction:

Vol(E k+1 ) ≤ exp -1 2 • (n + 1) • Vol(E k ) ( 15
)
Proof We give a proof to show that the successive ellipsoids computed by the method are actually decreasing by a ratio γ. First, let us put the update formula 13 into the form:

B k+1 = αB k + β(B k p)p T with: α = n √ n 2 -1 and β = n n + 1 - n √ n 2 -1 .
Let us now take the determinant of both sides.

det(B k+1 ) = det αB k + β(B k p)p T = det B k • αI n + βpp T = det B k det αI n + βpp T = det B k α n det I n + β α pp T
Using Sylvester's determinant identity:

det(I n + AB) = det(I m + BA) ∀A ∈ R n×m , B ∈ R m×n
the determinant on the right side of the equality can be express as:

det(B k+1 ) = α n det B k • 1 + β α p
But, From equation 14, we can see that p = 1. Therefore,

det(B k+1 ) det(B k ) = α n • 1 + β α ≤ exp -1 2(n + 1)
Hypotheses. In order to characterize the number of steps required for the algorithm to return an -optimal solution, three scalars and a point x c ∈ R n are needed:

• a radius R such that: • a scalar r such that:

X ⊂ B R (x c ) (16)
B r (x c ) ⊂ X (17)
• and another scalar V such that:

max x∈X f o -min x∈X f o ≤ V. ( 18
)
Those hypotheses are illustrated in Figure 21. This example corresponds to a Linear

Program setting where the feasible set is a bounded and not flat polyhedral set. The different balls are shown in shades of gray. More generally, the existence of those assumptions implies that the feasible set needs to be both bounded and not flat.

Unfortunately, when implementing MPC controllers, equality constraints are present, implying that the feasible set is flat on some dimensions. Therefore, we performed an initial equality constraint elimination (see Chapter 5.1.1).

The main result can be stated as:

Theorem 4 Let us assume that X is bounded, not empty and such that R, r and V are known. Then, for all ∈ R * + , the algorithm, using N iterations, will return x, satisfying:

f o (x) ≤ f o (x * ) + and x ∈ X ( -solution)
Furthermore, if we write n the dimension of the optimization problem, the number of steps, N , is polynomial in n and is of the form:

N = 2n • (n + 1) log RV r .
This result, when applied to LP, is historically at the origin of the proof of the polynomial solvability of linear programs. Its proof can be found at [START_REF] Khachiyan | Polynomial algorithms in linear programming[END_REF][START_REF] Nemirovski | Introduction to Linear Optimization[END_REF]. In order to produce formally verifiable code, we want to generate annotated code with Hoare triples, including function contracts expressing pre and post conditions. To do so, we present the development of ACSL theories related to optimization problems.

Building ACSL Theory Related to the Ellipsoid Method

Linear Algebra Axiomatization

In this section we give details about the ACSL theories we had to build in order to prove mathematical properties at code level. Indeed, the software analyzer takes as an input the annotated C code augmented with ACSL theories that define new abstract types, functions but also axioms, lemmas and theorems. The lemmas and theorems need to be proven but the axioms are always assumed to be true. For the SMT solver, properties at code level are usually harder to prove than lemmas within ACSL theories. Thus, our approach here was to develop the needed ACSL theories enough to be able to express and prove the main results used by the algorithm within the Linear Algebra Based ACSL Theory. In this ACSL theory, we defined new abstract types for vectors and matrices. We also defined functions that allow us to create a vector and a matrix from a C code pointer. Additionally, all the very well known operations have also been axiomatized (give a mathematical description to it) such as matrix multiplication, matrix addition, vector-scalar multiplication, scalar product, norm. This ACSL theory is automatically generated during the autocoding process of the project, thus, all the sizes of the vectors and matrices are known. Within this theory, we only defined functions that will create, from a C code pointer, objects of appropriate sizes (as illustrated in figure 22). ACSL code is printed in green and its keywords in red. The C code keywords are printed in blue and the actual C code is printed in black. Figure 22 presents the definition of two abstract types matrix and vector, the ACSL constructors for those types and the axiomatization of vectors addition. Ellipsoid Method Based ACSL Theory. In addition to defining new types for optimization problem, we also axiomatize the calculation of the vector constraint, feasibility, epsilon optimality, etc, . . . . Also, Ellipsoids and related properties are axiomatized, as presented in Figure 23. Within this theory, all the axioms and lemmas ACSL 1 # include " axiom_linalg . h " 2 /* @ axiomatic Ellipsoid { required for the main proof will be autocoded and proved. Before stating the main result, some preliminary results needs to be obtained. For instance, volume related properties are defined and proved (Figure 24). The proof of the main lemma for the method is stated in Figure 25. The first assumption listed on this lemma is the fact that the reals r, , V are strictly positive. It also assumes that the variable V satisfies the property [START_REF] Falcone | Predictive active steering control for autonomous vehicle systems[END_REF]. The lines 13 and 14 express the fact that we are assuming the ellipsoid Ell(P, x) to be a minimum localizer (meaning that any point outside of this ellipsoid cannot be better than x best, in term of cost or feasibility. Lines 15 to 17 specify the assumption for the scalar r to fulfill, corresponding to equation [START_REF] El Ghaoui | Inversion error, condition number, and approximate inverse of structured matrices[END_REF]. Finally, at line 18, we express the fact the the ellipsoid Ell(P, x) has a volume less than ( * r/V ) n . The conclusion of the lemma being the fact that x best represents an -solution to the optimization problem. In Figure 24, we show a axiomatization of a set. We wrote two axioms. One specifying that every set have a positive volume, and another one specifying that if a set B contains all the elements of a set A then the volume of B is at least equal to the volume of A. Then, we state a lemma expressing ACSL # include " axiom_linalg . h " /* @ axiomatic Optim { ...

Optimization Theory Axiomatization

To axiomatize an optimization problem, we intend to see it, independently of the method used to solve it, as a pure mathematical object. Our goal is to axiomatize it with enough properties, allowing us to state all the needed optimization-level properties at code level for the proof. Let us consider the second-order cone program, described in Eq. 6.

Encoding an SOCP. In order to fully describe an SOCP, we use the variables:

f ∈ R n , A =     A 1 . . . A m     , b =     b 1 . . . b m     , C =     c T 1 . . . c T m     , d =     d 1 . . . d m    
And also the vector m = n 1 . . . n m collecting the sizes of the vectors

A i • x + b i .
Furthermore, if na = m i=1 n i = 0, then, the SOCP we are considering is actually an LP. Using ACSL, we define a new type and a high level function, providing the possibility to create objects of the type "optim" (Figure 26).

When applying a method to solve an optimization problem, many concepts are important. The work here is to highlight those concepts and write a library translating those concepts into a formal specification language. That way, formally verifiable code could be produced, independently of the method and implementation used. The concepts of feasibility and optimality are being axiomatized. For this, given a second-order cone program, we gave an axiomatic definition for the vector constraint, the gradient of a constraint, the cost, optimal point (making the assumption that it exists and is unique), etc. For instance, Figure 27 trix, the size of the considered object needs to be know since it is hard-coded in the ACSL axiomatization. This does not represent an issue at this time since we already know all the sizes of the variable used (from the autocoder, see Section 5). Also, working with predefined and hard-coded size objects will help the analyzers proving the goals. The work presented here is generic and the code can be generated for any size of matrices and vectors.

Annotating a C code implementation of the Ellipsoid Method

We now give details about how we annotated the C code and the type of Hoare triples present in the code. For this, we adopted a specific technique. Every C code function will be implemented in a separated file. That way for every function, a corresponding C code body (.c) file and header file (.h) will be automatically generated. The body file contains the implementation of the function along with annotations and loop invariants. The header file contains the declaration of the function with its ACSL contract.

The first kind of Hoare triples and function contract we added to the code was to check the basic mathematical operations. For instance, Figures 29 and28 present the C code body and header files of the function computing the two norm of a vector of size two. Thanks to this contract, we can prove that the value returned by the function is indeed the two norm of the vector of size two associated with the input pointer.

Furthermore, we proved that the result is always positive or null, and that, assuming the corresponding vector is not equal to zero, the output is necessary strictly greater than zero. This last property being interesting when proving there are no division by zero (normalizing vectors). Then, once all the functions implementing elementary mathematical operations have been annotated and proven, we annotate the higherlevel C functions such as constraint calculations, gradient calculations, matrix and vector update, . . . . The Figures 30 and31 show the annotated C function and contract for the function "getp" that computes the vector p as described in equation ( 14), needed to perform the ellipsoid update.

C Code + ACSL The requires clause at line 17, expresses the fact that we are assuming the matrix of size (2, 2) defined by the pointer P minus to be invertible. Thanks to the assigns clause, we point out that this function only affects the variables p, temp2 and temp matrix. Finally, this triple lists three post execution properties that are supposed to be true. One expresses that the two norm of the vector of size two defined by # ifndef getp_lib # define getp_lib # include " a xi o m_ d ef _ li n _a l g . h " # include " socp / sizes . h " # include " getNorm_2 . h " # include " changeAxis . h " # include " scaleAxis . h " # include " getTranspose . the pointer p is equal to one. Another specifies that the matrix of size (2, 2) defined by the pointer P minus is invertible. The last expresses that the vector of size two defined by the pointer p is equal to the normalized multiplication of the transpose of the matrix of size (2, 2) defined by the pointer P minus times the vector of size two defined by the pointer grad. This latter property expresses equation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF]. In Figure 30, we show the implementation of the corresponding C code function and all the annotations required for the proof. For this, we decomposed the function into elementary mathematical operations performed via function called. That way, the code is traceable and easier to prove for the SMT solvers. Both function contracts have been successfully proved. One can wonder why we have the specification for the matrix P minus to be invertible after the execution of the function. Indeed, we assumed this same property to be true beforehand and specify that the variable it applies for would not change. These two properties making it trivial. In order to accelerate the verification process, SMT solvers need to be guided. Therefore, every needed property that is required somewhere else in the project should be listed, even if, logically speaking, it does not give any additional information. Thus, in order to guide the SMT solvers and to accelerate the verification, we added the ensures property at line 42.

Chapter IV FLOATING POINTS ANALYSIS

So far we presented how one can formalize existing results about the Ellipsoid Method and prove them at code level. We show now how the original algorithm can be modified in order to manage numerical errors. In this chapter we present a way of analyzing and controlling the numerical errors of a modified version of the Ellipsoid Algorithm.

Numerical errors of embedded systems can lead to tremendous catastrophes, therefore analyzing the numerical property of an algorithm that is meant to be embedded on a safety-critical CPS is an inevitable task. After recalling some past system failures we present in this chapter the work that has been done towards those issues. It turned out the issue was coming from internal variable used by the system to compute the current time in sec. The time from the internal clock of the system was stored in a variable t clock , encoded on a 24 bit fixed-point number and represented the time since system boot in decasecond. When current time in sec was needed, the program was performing equation [START_REF] Feron | From control systems to control software[END_REF], multiply as expected the internal variable by 0.1.

Past System Failures and Motivation

t (sec) = t clock × fl(0.1) (19)
Unfortunately, the constant 0.1 cannot be represented using 24 bit fixed-point numbers (it is also not the case for floating-point numbers). Therefore, this constant was chopped and induced an error on the time computation that was not expected by the engineers. The small rounding error, when multiplied by the large number giving the 19)). Equivalently, because a Scud travels at about 3, 750 mph, it corresponds to an error of more than 1, 600 f t.

Ariane 5 Rocket

This error is known as one of the most expensive floating-point error. It caused a damage worth half a billion dollars. The horizontal velocity of the rocket encoded on a 64 bit floating-point number was converted to a 16 bit signed integer. The problem here being the fact that this number was larger than the largest integer representable on a 16 bit signed integer. The conversion failed and the software ended up triggering a SAFETY mode, switching to a backup computer. Unfortunately, the same error happened and it was misinterpreted as a scenario where aggressive control input from the motor was needed. Only 40 seconds after ignition, at an altitude of about 12000 f t, the rocket went out of its nominal trajectory, broke up and exploded.

Motivations

As seen from the past sections, it is undeniable that floating-points errors could lead to system failure. We now argue that, in addition to that, the nature of the algorithm we are developing are particularly sensitive to numerical errors and thus analyzing its numerical properties is even more important. A large majority of optimization algorithms are iterative algorithms. Meaning that within its structure, lies a main loop in which the n -th approximation is derived from the previous ones. Each iterate have given properties, which are necessary for the algorithm convergence. If we consider the operations performed by the computer to have errors, those errors could potentially be added-up at each iteration and could lead to a software error.

Therefore, studying numerical stability is indeed a need and would give guarantees that are highly appreciated.

Controlling the Condition Number

In this section we explain how the condition number of the successive Ellipsoids of the algorithm can be bounded through the execution of the code and why it is important.

We recall below the definition of the condition number of a non-singular matrix.

Definition 5 (Condition Number of a Matrix)

Let A ∈ R n×n a non-singular matrix. We define the condition number of the matrix A the scalar k(A) such that:

k(A) = A • A -1
By extension, we talk about the condition number of an Ellipsoid Ell(B, c) by taking the condition number of the matrix B, k(B).

Bounding the condition number of the matrix B is fundamental and represents the main argument of the algorithm numerical stability. Unfortunately, for the original algorithm, no reasonable bound on k(B) can be found. Therefore, we slightly modified the ellipsoid algorithm to make it able to correct the current ellipsoid E i in the case where its condition number had become too high (ellipsoid too flat). That way we can control the condition number of B. Additionally, we made sure that this correcting step, when it occurs, does not break the convergence of the algorithm and its semantics described in Section 3.2. In this section, we give details about the modification performed on the original algorithm.

Bounding the Singular Values

When updating the matrix B i by the usual formulas of the Ellipsoid Algorithm, B i evolves according to

B i+1 = B i • D i , (20)
where n -1 singular values of D i are n/ √ n 2 -1, and one singular value is n/(n + 1).

It follows that at a single step the largest and the smallest singular values of B i can change by a factor from [1/2, 2].

Proof From the update equation (20), we have the property:

σ max (B i+1 ) = σ max (B i • D i ) = B i • D i 2 ≤ B i 2 • D i 2
because the two norm is a consistent norm and

A 2 = σ max (A) ∀A ∈ R n×n Therefore, σ max (B i+1 ) ≤ σ max (B i ) • σ max (D i ) = σ max (B i ) • n √ n 2 -1 .
Which implies that:

σ max (B i+1 ) ≤ 2 • σ max (B i ) ∀n ∈ N , n ≥ 2
For σ min , let us take first the inverse of the update equation (20) (we know that all the matrices are indeed invertible). We have:

B -1 i+1 = D -1 i • B -1 i . Thus, σ max (B -1 i+1 ) = σ max (D -1 i • B -1 i ) ≤ σ max (D -1 i ) • σ max (B -1 i )
Using the fact that for an invertible matrix A, we have:

σ max (A -1 ) = 1 σ min (A) we conclude that: 1 σ min (B i+1 ) ≤ 1 σ min (D i ) • 1 σ min (B i )
Rearranging this last equation, we end up with:

σ min (B i+1 ) ≥ σ min (D i ) • σ min (B i ) = σ min (B i ) • n n + 1 .
Again, if we consider all possible sizes, we end up with:

σ min (B i+1 ) ≥ 1 2 • σ min (B i ) ∀n ∈ N , n ≥ 2
Let us argue now that one can bound the singular values of the matrix B i throughout the execution of the program.

Minimum Half Axis: First, we claim that if σ min (B) is less than r /V then the algorithm has already found an -solution. The scalar being the wanted precision and the scalars r and V being defined in Section 3.2.

Let us assume σ min (B) < r /V . In this case, E i is contained in the stripe between two parallel hyperplanes, the width of the stripe being strictly less than 2 • r /V

and consequently E i does not contain X θ (defined in equation ( 21)), where x * is the minimizer of f o and θ = /V . This argument being a consequence of the fact that X θ contains a ball of radius r /V (by definition).

X θ = θX + (1 -θ)x * = {θz + (1 -θ)x * , z ∈ X} (21)
Consequently, there exists z ∈ X such that y

= θz + (1 -θ)x * ∈ X but y / ∈ E i ,
implying by the standard argument that the best value f + of f processed so far for feasible solutions satisfies

f + ≤ f (y) ≤ f (x * ) + θf (z -x * ) which implies that f + ≤ f * + .
We can thus stop the algorithm and return the current best point found (feasible and smallest cost).

Maximum Half Axis:

We argue in the section that one can modify the original ellipsoid algorithm in order to bound the value of the maximum singular value of B i . When the largest singular value of B i is less than, say, 2R √ n + 1, we carry out a step as in the basic ellipsoid method. When this singular value is greater than 2R √ n + 1, we take some time to "correct "B i , namely, to pass from

E i = B i X to E + i = B + i X in
such a way that E + i is a localizer along with E i , meaning that E i ∩ X ⊂ E + i ∩ X. In addition to that, we have the following properties:

(a) The volume of E + i is at most γ times the volume of E i ;

(b) The largest singular value of

B + i is at most 2R √ n + 1.
For this, let us define σ = σ max (B i ) > 2R √ n + 1 and let e o being corresponding direction and index. We then consider the matrix G such that:

G = diag n/(n + 1), √ n + 1/σ, . . . , √ n + 1/σ
We conclude then this case by performing the below update on B i and c i :

B i+1 = B i • G and c i+1 = c i -(e T o c i ) • e o (22)
Figure 34 shows an illustration of such a correction. The unit ball being the feasible

σ min (B i ) ≥ 1 2 r V = r 2V and σ max (B i ) ≤ 2 × 2R √ n + 1 = 4R √ n + 1
We have the following element of proof. First, let us compute the volume of the ellipsoid obtained using this correction step.

Vol(

E + i ) = det(G) • Vol(E i ) So, Vol(E + i ) = (1 + 1/n) (n-1)/2 √ n + 1/σ • Vol(E i ) ≤ exp(1/2) 2R • Vol(E i ).
Usually, R being quite large, we have exp(1/2)/(2R) ≤ γ. If we look at this equation more carefully, we can see that having R ≥ exp(1/2) = 1.65 implies that exp(1/2)/(2R) ≤ γ for all dimension n (we recall that γ depends on n) greater than 2. Hence, as R is usually a substantial scalar (radius of a ball including the feasible set), it will most likely be true. If it is not the case, then we could just make R bigger until we have the wanted property.

Corresponding Condition Number

Let us now see how this process impact the condition number of B. First, from the definition of the condition number we have:

k(B) = B • B -1 = σ max (B) σ min (B)
Because of the very well known property of the two norm: A = σ max (A). and for A non-singular:

A -1 = σ max (A -1 ) = 1 σ min (A)
.

Thus, by bounding the singular values of B, we concluded on a bound on the condition number of the matrix B.

k(B) ≤ 2 1/2 • 2R √ n + 1 r /V = 8R √ n + 1 r /V and B = σ max (B) ≤ 4R √ n + 1

Corresponding norm on c

At each iteration we know that we have:

x * ∈ Ell(B, c) Thus, x * -c = Bu ≤ B • u ≤ B , for some u ∈ B 1 (0) Finally, c ≤ R + c + B

Consequences on Code

In this section, we explain how to implement this correcting step and give the tools to verify it. In order to detect ellipsoids with large semi-major axes, we need to compute the largest singular value σ max of the current matrix B k . However, performing a singular value decomposition would be way too expensive and slow (this decomposition being performed at each iteration). On the other hand, because σ max (P ) is equal to the two norm, we compute an over approximation of the two norm, the Frobenius norm. The latter one being an over approximation of σ max (P ) and extremely fast to compute, it represents the best option. Because in practice, the successive ellipsoids stay within acceptable bounds, no semi-major axis should be detected. Thus, if we know that the Frobenius norm of the matrix B i is less than the maximum value acceptable for semi-major axis, we can conclude that the ellipsoid is well conditioned.

The equations ( 23), ( 24) and ( 25) show more details about matrix norm equivalence and semi-major axis. We have the following very well known properties corresponding to matrix norms:

A 2 ≤ A F ≤ √ n A 2 ∀A ∈ R n×n (23)
with:

σ max (A) = A 2 ∀A ∈ R n×n . ( 24
)
Thus, we have the property needed below:

σ max (A) ≤ A F ∀A ∈ R n×n . ( 25
)
The mathematical definition of the Frobenius norm of a matrix is presented in definition 6.

Definition 6

The Frobenius norm of a matrix A ∈ R n×n is:

A F = n i=1 n j=1 a 2 i,j
We axiomatized the Frobenius norm of a vector by being equal to the vector two norm of the "vectorized" matrix. We "vectorize" a matrix by concatenating all its rows in 1 # include " a x i o m _ d e f _ l i n _ a l g . h " 2 # include " getNorm_256 . h " 3 # include " socp / sizes . h " In the case where a too large semi-major axis is detected, The direction e in which this axis lies is also needed. Therefore, we perform a power's iterative algorithm in order to compute this information. Other algorithms could be used in order to get the full decomposition of the current matrix, but does not represent useful information for our use. Work already have been made concerning the formal verification of matrix decomposition algorithms [START_REF] Roux | Formal proofs of rounding error bounds[END_REF].

Propagating the rounding errors through the algorithm

Let F denotes the set of all floating-point numbers and R the set of reals. We use standard notation for rounding error analysis [START_REF] Roux | Closed loop analysis of control command software[END_REF][START_REF] Rump | Verification of positive definiteness[END_REF][START_REF] Rump | Error estimation of floating-point summation and dot product[END_REF], fl(•) being the result of the expression within the parenthesis computed in rounding to nearest. We write the relative rounding error unit u and the underflow unit eta. For IEEE 754 double precision (binary64) we have u= 2 -53 and eta= 2 -1074 . We present in this section an analysis targeting the numerical properties of the ellipsoid algorithm. Contributions already have been made concerning finite-precision calculations within the ellipsoid method [START_REF] Khachiyan | Polynomial algorithms in linear programming[END_REF]. However, this work only shows that it is possible to compute approximate solutions without giving exact bounds; it remains very theoretical and only applied to Linear Programming (LP). Also, the analysis performed considers abstract finite-precision numbers and floating-points are not mentioned. Thanks to the analysis performed in this section, using the IEEE standard for floating-point arithmetic and knowing exactly how the errors are being propagated, we would be able to check a posteriori the correctness of the analysis using static analyzers [START_REF] Goubault | Static Analyses of the Precision of Floating-Point Operations[END_REF][START_REF] Goubault | Static Analysis of Finite Precision Computations[END_REF][START_REF] Putot | Static Analysis-Based Validation of Floating-Point Computations[END_REF].

Preliminaries

Within this algorithm, we focus our attention on the update formulas ( 12), ( 13) and ( 14), allowing us to update the current ellipsoid into the next one. This program dealing with ellipsoids, when investigating its numerical property, the condition number is a decisive information.

Theorem 5 [Matrix perturbations and Inverse] Let A be a non-singular matrix of R n×n and ∆A a small perturbation of A. Then, from [START_REF] El Ghaoui | Inversion error, condition number, and approximate inverse of structured matrices[END_REF], we know that,

(A + ∆A) -1 -A -1 A -1 ≤ k(A) ∆A A (26)

Norms and Bounds

To successfully perform the numerical analysis of the algorithm, we need to know how "big" the variables can grow within the execution of the algorithm. Indeed, for a given instruction, the errors due to floating-point arithmetic are usually proportional to the value of the variables.

Bound on variables c, B, k(B) and p.

For the variable p, we have:

p = B T e √ e T BB T e = 1 ( 27
)
After modifying the original algorithm in the way described in Section 4.2 we have the following results:

B ≤ 4R √ n + 1 (28) k(B) ≤ 8RV √ n + 1 r (29) c ≤ R + x c + B ( 30
)
Where n, R, r, V, x c and are the variables described in Section 3.2.

Floating-Point Rounding of Elementary Transformations

In this section, we express the floating-point errors taking place when performing the update formulas ( 12) and [START_REF] Champion | Formal methods for the analysis of critical control systems models: Combining non-linear and linear analyses[END_REF]. For this, we present first the error analysis for basic operations appearing in the algorithm.

Rounding of a Real.

Let z ∈ R z = fl(z) = z + δ + η with |δ| < u and |η| < eta/2
Product and Addition of Floating-Points. Let a, b ∈ F.

fl(a × b) = (a × b)(1 + 2 ) + η 2 fl(a + b) = (a + b)(1 + 1 ) with: | 1 | < u, | 2 | < u, |η 2 | < eta and 2 η 2 = 0 Reals-Floats Product. Let z ∈ R and a ∈ F, |fl fl(z) • a -z • a| ≤ |z||a| • u + |a| • 2u(1 + u) Scalar Product. Let a, b ∈ F n . We define, a, b = n i=1 a i b i and |a, b| = n i=1 |a i b i |.
We have then:

|fl a, b -a, b | ≤ A n |a, b| + Γ n
With:

A n = n • u 1 -n • u and Γ n = A 2n eta u = 2n • u 1 -2n • u eta u

Multiplication of Reals and Floating Scalar Product.

Let z ∈ R and a, b ∈ F n . We have the following property:

|fl fl(z) • fl a, b -z × a, b | ≤ |a, b| • |z| • 2u(1 + n) + |a, b| • 4u (31)
Now that the propagation of the numerical errors through the elementary transformations have been presented. In the same fashion, we identify the kind of operations performed to update the current ellipsoid. For this, we define ∆ B , ∆ B -1 and ∆ c representing the floating-point errors, such that:

∆ c = fl(c + ) -c + (32) ∆ B = fl(B + ) -B + (33) ∆ B -1 = fl(B + ) -1 -B + -1 (34)
and assume that after performing the floating-point analysis we found E B and E c such that:

|(∆ B ) i,j | ≤ E B ∀i, j ∈ [1, n] (35)
and

|(∆ c ) i | ≤ E c ∀i ∈ [1, n]. (36)
We dedicated Section 4.3.3 to the computation of E c and E B .

Rounding Error on c + . Knowing how the errors are being propagated through elementary transformations, we want to compute the error for a transformation similar to what happen for the vector c update. For each component of c, we have:

c + i = c i -1/(n + 1) • Row i (B), p .
Therefore, the operation performed, in floating-point arithmetic is:

fl c + fl fl(z) • fl a, b .
Ell B + , c + ⊂ Ell λ • fl(B + ), fl(c + ) ⇐⇒ fl B + -1 • B + u + c + -fl(c + ) ≤ λ ∀u, u ≤ 1 Proof Let us denote B 1 (0) the unit Euclidean ball. (i.e. B 1 (0) = {x ∈ R n : x T x ≤ 1}
). Using definition [START_REF] Boyd | Convex optimization[END_REF], we have the following equivalent statement:

∀ u 1 ∈ B 1 (0) , z = B + u 1 + c + → ∃ u 2 ∈ B 1 (0) , z = λ • fl(B + )u 2 + fl(c + )
Let us now reformulate the second part of the statement.

∀ u 1 ∈ B 1 (0) , z = B + u 1 + c + → ∃ u 2 ∈ B 1 (0) , λ • fl(B + ) -1 • z -fl(c + ) = u 2
which is equivalent to:

∀ u 1 ∈ B 1 (0) , z = B + u 1 + c + → λ • fl(B + ) -1 • z -fl(c + ) ≤ 1
Let us now replace z by its formula in the second part of the statement.

∀ u 1 ∈ B 1 (0) , λ • fl(B + ) -1 • B + u 1 + c + -fl(c + ) ≤ 1
Putting on the other side of the inequality λ, we end up with the wanted property.

∀ u ∈ B 1 (0) , fl(B + ) -1 • B + u + c + -fl(c + ) ≤ λ

Analytical Sufficient Conditions for Covering

From Lemma 2, one can see that the calculation of a widening coefficient λ highly depends on the accuracy of the matrix fl(B + ) -1 . Therefore, we also need to compute a number E B -1 such that:

|(∆ B -1 ) i,j | ≤ E B -1 ∀i, j ∈ [1, n]. (43)
The quantity (B + ) -1 is not used explicitly in the algorithm and its floating-point error could not be evaluated by numerically analyzing the method. Instead, we will use perturbation matrix theory [START_REF] El Ghaoui | Inversion error, condition number, and approximate inverse of structured matrices[END_REF] and Theorem 5. This will give us an upper bound on E B -1 given E B , the norm of B and its condition number. The result is stated in the following lemma.

Lemma 3 [Widening -Analytical Sufficient Condition] λ ≥ 1 + k(B) B √ n • √ n • k(B)E B + E c + k(B) B nE B E c =⇒ Ell B + , c + ⊂ Ell λ • fl(B + ), fl(c + ) .
Proof First, in order to prove this theorem, we need to evaluate ∆ B -1 . Using equation [START_REF] Jadbabaie | Stabilizing receding horizon control of nonlinear systems: a control lyapunov function approach[END_REF],

∆ B -1 ≤ k(B) B -1 B ∆ B = k 2 (B) B 2 ∆ B But, ∆ B ≤ ∆ B F ≤ nE B
So, we have:

∆ B -1 ≤ k 2 (B) B 2 nE B
let us now define the three constants below:

I = fl B + -1 B + ,
In that case, if N denotes the original number of iteration needed, the algorithm implementing the widened ellipsoids will require:

N λ = N 1 -n(n + 1) log(λ) iterations ( 45
)
Proof We recall that the algorithm converges if and only if the volumes of the successive ellipsoids are decreasing.

Vol Ell(B k+1 , c k+1 ) < Vol Ell(B k , c k )
But we know when using the original update process of the Ellipsoid Method we have:

Vol Ell(B k+1 , c k+1 ) ≤ γ • Vol Ell(B k , c k )
With γ = exp(-1/(2(n + 1))). Also, thanks to equation [START_REF] Boyd | Linear controller design: limits of performance[END_REF], we know that

Vol Ell(λ • B k , c k ) = λ n/2 • Vol Ell(B k , c k )
Therefore, the corrected algorithm, with widening coefficient λ will converge if and only if:

λ n/2 • γ < 1
Which is equivalent to:

λ < exp 1 n(n + 1)
Therefore using equation [START_REF] Boyd | Linear controller design: limits of performance[END_REF],

Vol(E k+1 ) = λ n/2 • Vol(E k+1 )
So, we have:

γ λ = Vol(E k+1 ) Vol(E k ) = Vol(E k+1 ) Vol(E k+1 ) Vol(E k+1 ) Vol(E k ) = λ n/2 • γ
In order to end up at the final step with an ellipsoid of the same volume, we need:

Vol(E o ) • γ N λ λ = Vol(E o ) • γ N
Which implies:

N λ • log(γ) + n 2 log(λ) = N • log(γ)
Replacing γ by its value, using property 1:

N λ • -1 2(n + 1) + n 2 log(λ) = N • -1 2(n + 1)
And thus, finding the wanted formula:

N λ = N 1 -n(n + 1) log(λ)
.

An application of this framework is presented in Section 5.4.

Chapter V AUTOCODER AND CLOSED LOOP MANAGEMENT

In this Chapter, we give details about the closed-loop management of optimization algorithms. We show how to extract convergence guarantees for parameterized optimization problems. Following this, we present how to use the autocoder GENEMO that we built and how to modify the output generated code. Finally, real-time simulations using generated code for different systems are presented.

Closed Loop Management -Sequential Optimization Problems

In this section, we are interested in extracting convergence guarantees for a class of optimization problems that are going to be used online. For this, considering the MPC controller input as a parameter, we end up with an optimization problem with parameterized constraints and cost. The work presented in this section is about studying how this parameter affects the constraints and the cost at each iteration and trying to find hypothesis for this parameter to fulfill in order to conclude on the convergence of the algorithm for every points along the trajectory.

First, we focus on the special case where only linear constraints are present. Following this, another section will be dedicated to the more general setting of SOCP constraints.

Parameterized Linear Constraints

Let us write the optimization problem we are aiming to solve in real-time. The vector X ∈ R nx denotes the decision vector. minimize

X f o (X) subject to AX ≤ b (46)
Let us assume the initialization of this optimization problem is being done such as equation [START_REF] Tsiotras | Goddard problem with constrained time of flight[END_REF], where S is a full rank matrix. Usually, S represents a selector matrix and is of the form:

[I p O p×(nx-p) ].
In that case it is obviously full rank. x o denotes the input of the controller. We write xo to account for the fact that x o will change from one optimization problem to another.

S • X = xo (47)
We decompose and separate the equality and inequality constraints hidden behind the original matrix A and vector b. That way, Problem 46 can be written as: minimize

X f o (X) subject to A eq • X = b eq A ineq • X ≤ b ineq S • X = xo (48)
The idea here is to project all the equality constraints in order to eliminate it while keeping track of the variable parameter, xo . We know there exist matrices M , A 1 and A 2 such that for all vector X ∈ R nx satisfying the equality constraints of Problem 48, there exists a vector Z ∈ R nz such that:

X = A 1 • b eq + A 2 • xo + M • Z (49)
M being a matrix formed by an orthonormal basis of the null space of the matrix   

A eq S    . d denotes the total number of equality constraints (number of rows of A eq

A f =                 -1 1 1 1 1 -0.5 0 1 -1 0 0 -1                 ; φ(A f,l ) =                 1 1 1 1 1 1                 ; b f,o =                 1 2 1 1.5 0.5 0.5                 ; r o = 0.5.
We state an important Lemma, illustrating the fact that for every point along the

∀x o ∈ R n s.t. x o 2 ≤ r o , P l ⊂ P xo ⊂ P u Proof Let us take x o such that x o 2 ≤ r o .
First, let us establish a very simple inequality using Cauchy-Schwarz inequality:

|(A f,l • x o )(i)| = |row(A f,l , i) T • x o | ≤ row(A f,l , i) 2 • x o 2 ≤ φ(A f,l )(i) • r o ∀i (55)
Therefore,

-φ(A f,l ) • r o ≤ A f,l • x o ≤ φ(A f,l ) • r o Then, if x ∈ P l we have A f • x ≤ b f,o -r o • φ(A f,l
) . Using the inequality below, it is clear that it implies x ∈ P xo . The same way, assuming that x ∈ P xo and using the inequality below, it is clear that we have x ∈ P u .

The work of the control engineer consists in finding three scalars r, R and V such that:

∃ z1 such that B( z1 , r) ⊂ P l , (56)
∃ z2 such that P u ⊂ B( z2 , R) , (57)
V ≥ max

z∈P xo f o (x sol + M z) -min z∈P xo f o (x sol + M z) , ∀ xo ≤ r o . (58)
For the first scalar, r, one can compute a numerical value by running an off-line optimization problem finding the largest ball inside P l . If no solution can be found, the value of r o need to be decreased, and we repeat the process until finding an acceptable r o and radius r. Further information about finding the largest ball in a polytope can be found at [START_REF] Boyd | Convex optimization[END_REF].

As a consequence of equality constraint elimination, and recalled in Eq. 49, we have the relation below between the original decision vector X and the projected one Z:

X = A proj •    b eq xo    + M Z (59)
We decompose now the decision vector X into two parts, x and u. The part u being bounded due to constraints in the original optimization problem. If no constraints on u were originally present, one can add some making the problem bounded and simpler to analyze. The point here being that from bounded variables within the vector X, one can conclude on bounds on the projected vector Z. There is no need to have original bounds particularly on the collection of future inputs u. Rewriting equation 59 yield:

   x u    =    A 11 A 12 A 21 A 22       b eq xo    +    M 1 M 2    Z (60)
Following this, we would have:

Z = M -1 2 • u -A 21 b eq -A 22 xo (61)
and therefore assuming again that xo ≤ r o , one can compute a value of R such that:

Z ≤ M -1 2 • u + A 21 b eq + A 22 r o = R (62)
On the other hand, from the physical meaning of the variables and the constraints of the optimization problem, one can conclude on bounds in which the variables should live in, and therefore find a lower bound for V .

Additionally, when the original problem is too hard to analyze one can add constraints in order to make the analysis easier. Adding constraints in order to find acceptable values for the scalars R and V more easily is a very efficient technique but can have a negative effect on the value of r. When adding constraint to a given optimization problem, its feasible set get smaller and it is harder to find a inscribed ball that lies inside.

We now state the general lemma, giving an upper bound on the number of iteration throughout a trajectory.

Fact 2 [MPC Ellipsoid Method Convergence]

Let us assume we want to run the problem 48 on-line in order to implement a receding horizon control.

That way, using the Ellipsoid Method and initializing the first Ellipsoid by B( z2 , R), the method will find an -solution using N it iteration for all x o such that x o ≤ r o , with:

N it = 2 • n z (n z + 1) • log R r V and n z = dim Null A eq S T Proof
Let us assume, we are running this problem online and we are currently trying to solve this problem for a given point. We assume x o such that x o ≤ r o . Let us write X f to denote the feasible set of the current optimization problem to solve. that way, we know that we have Eq. 56 Thus, by definition of r, we know that:

∃ z1 such that B( z1 , r) ⊂ P l ⊂ X f
The same way, because we assume R satisfies Eq. 57, we also have:

∃ z2 such that X f ⊂ B( z1 , r)
We also have the scalar V meeting the wanted requirement. Finally, we see that all the needed hypothesis are fulfilled and we can conclude that the returned point will indeed be -optimal using N it iteration thanks to Theorem 4.

Linear Programs: In the case of Linear Programs, the method developed is identical and having a linear cost, the scalar V is easily found by running the two optimization problems below.

minimize xo,z c T • (A 1 b eq + A 2 xo + M z) A f • z ≤ b f,o + A f,l • xo xo 2 ≤ r o maximize xo,z c T • (A 1 b eq + A 2 x o + M z) A f • z ≤ b f,o + A f,l • x o x o 2 ≤ r o Having f o (x) = c T x.

Second-Order Conic Constraints

In the previous section, we used the fact that we had linear constraints to construct extreme polyhedral sets and conclude on geometric characteristics on feasible set. In this section we aim at giving ways of computing the same needed constants for the general setting of second-order constraints. Let us now consider we want to solve a model predictive control optimization problem of the form: minimize X f o (X)

A i • X + b i 2 ≤ c T i X + d i , i = 1 . . . m A eq X = b eq S • X = xo (63)
First, let us assume that no equality constraint are hidden in the second-order constraints (if it is not the case, a very simple analysis will confirm it and one can extract those equality constraints and put it into the couple (A eq , b eq )). Unfortunately, even after eliminating the equality constraints of problem 63, in order to find a value for the scalar r, we need to find the largest balls inside a second-order cone, which is not an easy task. To get around this, we use the equivalence of norm in finite dimensions (see Eq. 64, 65 and Figure 38), noting that • 1 and • ∞ can be expressed with linear constraints. We perform a linear relaxation, by replacing the second-order constraints with stronger constraints that can be expressed with linear constraints. That way, the feasible set becomes a polyhedron and the previous tools could be applied. For instance, problem 66 represents a linear relaxation of the original problem 63.

x ∞ ≤ x 2 ≤ √ n • x ∞ (64) n -1/2 • x 1 ≤ x 2 ≤ x 1 (65) minimize X f o (X) √ n • A i • X + b i ∞ ≤ c T i X + d i , i = 1 . . . m A eq X = b eq S • X = xo (66)
Figure 39 shows an illustration of a second-order linear relaxation, along with the largest ball inside the relaxed polyhedral set.

Running Time Evaluation

In this chapter, we analyze the online use of optimization algorithms in order to extract convergence guarantees for every point along the trajectory. In Section 5.1, we showed how one can decomposed the online optimization problem as a parametrized one and compute an a priori upper bound on the number of iterations needed to be performed for every points along the trajectory. Those algorithms are run online and have very strict time constraint on their execution time. Therefore, the key information here is the running time of those algorithms, which could be directly evaluated from the number of iterations. From the floating-point operations per with the use of the keyword "subjectTo". A constraint is of the form:

"constraintj: ... ; ".

The same way we used a receding index variable to implement the sum function, we also used such index to define sliding constraints. Figure 42 presents an example of constraints written in GENEMO. Also, with the use of the keyword "Information", one can incorporate the a priori information needed described in Chapter 3 (see Fig- the GENEMO file could implement an MPC controller but also a single optimization problem. In the latter case, the goal being to solve and give the solution of a given single optimization problem and the sections described below (Input/Output sections) should be ignored. On the other hand, if the purpose of using this autocoder was to automatically generate C code implementation of MPC controller, the user has to complete these sections in order to give information about the Input/Output aspect of the MPC controller. Figure 44 presents the input and output sections completed for the spring-mass example. Comments can be added along the code thanks to the sign #. Only single line comments are managed.

Internal Aspect of The Autocoder

In this section we give details about the internal aspect of the autocoder. The user first has to formalize the MPC controller desired and write the corresponding input file detailed in Section 5.3. In addition to this, the user has to specify the implementation of the method used to solve optimization problems. By default, this autocoder has been created in order to generate C code implementation of the Ellipsoid Algorithm.

Nevertheless, the C code generator is generic and could be applied to any algorithms.

In order to generate C code implementation of another solving method, the user has to give and formalize the augmented AST (Abstract Syntax Tree) using a python library wrote beforehand. This library implements all the necessary properties relative to trees, AST and C code generation. Additionally, this library has been written in a way to facilitate the augmentation for semantics and Hoare Triples to those AST.

By default, the usual control structures will be autocoded along with annotations.

For instance, when creating a node relative to a for loop using the AST library, the autocoder will by default generate the C code implementation of the for loop structure along with trivial loop invariants, loop variant and assigns clauses. That way, a considerable amount of annotations, required for the proof, are handled directly by the autocoder which makes the development process faster. We present in this section an application of this framework for the spring-mass system described in Section 2.3, Figure 14. For this, we use again a discretization period of T = 0.1 sec and an horizon of N = 10.

We recall the discrete state-space realization matrices A and B of this system below:

A =    1 T -T 1    ; B =    0 T    .
Let us assume the initial point at t = 0 is:

x o = 2 -1 T ,
and the matrix Q equal to:

Q =    10 0 0 1    .
We want to generate C code in order to perform the below MPC for this system. The goal is to stabilize the origin as fast as possible (minimizing the norm of the state along the future trajectories).

minimize

X=[x,u] N k=1 Q • x k x k+1 = Ax k + Bu k , k = 1..N -1 -5 ≤ u k (1) ≤ 5 , k = 1..N -1 -10 ≤ x k (1) ≤ 10 , k = 1..N -10 ≤ x k (2) ≤ 10 , k = 1..N x 1 = x i (68)
We use an accuracy of = 0.1 and using the method developed in the previous sections, we found: r = 4.55, R = 18.19, V = 141.42.

From that, we conclude that the program will need N it = 1556 iterations in order to compute an -solution. Using double precision floating points, we can also conclude on a widening coefficient λ = 1 + 5. system with a inner feedback controller, and a discretization step of T = 0.5 sec. The resulting system is a stable linear system that we control using MPC. The problem we are trying to solve is supposed to model a landing of the 3 DOF helicopter. We initialize the system with an angle of 25 deg in elevation and 15 deg travel and we want to make the system go back and land at the origin while avoiding the ground.

The ground is the area below λ = 0. Thus, the constraint in order to avoid the ground is a combination of the elevation and pitch angle and can be formulated as:

h sin(λ) ± d sin(φ) >= 0
We linearized those constraints, assuming small angles perturbations, and naturally end up with linear inequalities of the form: A obs • x ≤ b obs . The MPC controller that we want to implement is formalized in problem 69.

minimize

X=[x,u] N k=1 x k x k+1 = Ax k + Bu k , k = 1..N -1 u k ≤ 50 , k = 1..N -1 x k ≤ 200 , k = 2..N A obs • x k ≤ b obs , k = 2..N -40 ≤ x k (2) ≤ 40 , k = 1..N -1 x 1 = xo (69)
Assuming that x o ≤ 30, we found, using the method developed in section 5.1 a radius 8.0612 (running an off-line optimization problem that finds the largest ball inside P l ). For R, using the method described in 5.1 and equation 62 we found R = 322. From the problem formulation, one can see that if the problem is feasible, the sum of the norm of the successive x is supposed to be minimized. That way, the worst case is when x o had the largest norm, and when the system stay at this point throughout the whole trajectory (or keep the same norm at least). That way, we should have a V of:

V = N • x o ≤ N • 27 = 162 (70)
For all those cases, we can also conclude on a widening coefficient lambda to apply in order to control floating-point errors. The method we use to control the floating points error is affecting the volume ratio of the method. Therefore, in order to account for it, we need to do more iteration depending on the value of the widening coefficient We present in this section an application of this framework for Quadcopter Drones.

For this we used a linear discrete-time representation of the drone of the form:

x k+1 = Ax k + B 1 u k + B 2 u k+1 + f d. (71)
The vector input u collects the vector thrust and the vector x regroups the position and velocities and the drone along each axes. The matrix A, B 1 and B 2 are described below. The mass of the system is written m , T s denotes the sampling period and g, the gravity constant.

In this section, we implement a path-planner based on the autocoder we built. For this, we generate beforehand a nominal trajectory using convex optimization techniques. We formulate an optimization problem, to compute an optimal trajectory for the drone to go from a point A to a point B, while satisfying constraints.

One of the constraints to respect is the fact that the input of this system, the thrust the fact that we do not want the drone to rotate too much because non-linearities phenomenon will rise and it will be difficult to control it. One of the problem, was the fact that this constraint is not convex. For that, we perform a lossless convexification step adding slack variables in order to encode this constraint in a convex manner.

More details about lossless convexification can be found at [START_REF] Ac ¸ıkmes ¸e | Lossless convexification of a class of optimal control problems with non-convex control constraints[END_REF]. We end up with the optimization problem presented in equation 72. Following the method described in Section 5.1, we found the values for the variable r, R and V listed in Table 3 52.

A =                 1 T s 0 0 0 0 0 1 0 0 0 0 0 0 1 T s 0 0 0 0 0 1 0 0 0 0 0 0 1 T s 0 0 0 0 0 1                 ; f d =                 0 0 0 0 0.5 • T s • T s • g T s • g                 B 1 =                 (1/3) • T s • T s /m 0 0 (1/2) • T s /m 0 0 0 (1/3) • T s • T s /m 0 0 (1/2) • T s /m 0 0 0 (1/3) • T s • T s /m 0 0 (1/2) • T s /m                 B 2 =                 (1/6) • T s • T s /m 0 0 (1/2) • T s /m 0 0 0 (1/6) • T s • T s /m 0 0 (1/2) • T s /m 0 0 0 (1/6) • T s • T s /m 0 0 (1/2) • T s /m                 minimize X=[x,s,u] s x k+1 = Ax k + B 1 u k + B 2 u k+1 + f d. , k = 1..N -1 u k ≤ s(k) , k = 1..N T min ≤ s(k) ≤ T max , k = 1..N -u k (3) ≥ s(k) • cos(φ max )
x k (5) ≥ 0 In Chapter 1, we introduced the subject, gave details about the needed elements and presented a literature review. Following this, we recalled some important results of control theory for understanding the application we are targeting in this research in Chapter 2. Chapter 3 presents an axiomatization of second-order cone programs, a subset of convex optimization problem, using the specification language ACSL. Additionally, annotations for numerical algorithms solving these problems were proposed. We focused our attention on the Ellipsoid Method. A numerical analysis of the method has been presented in Chapter 4. We first show how to modify the original Ellipsoid Algorithm to properly bound the programs variables. Then the propagation of numerical errors due to floating-point calculations through the operations performed by the program is presented. Finally, Chapter 5 gives mathematical evidence of the online use of these optimization algorithms. Closed-loop system is handled and receding-horizon controllers can be autocoded with soundness guarantees.

x 1 = x i ; u 1 = u i ; x N = x f ; u N = u f (72)
As it was recalled earlier, one of the purposes of this research was to improve V&V techniques while developing critical software. Nevertheless, this work should not be seen as a replacement for unit or integration testing. Instead, we intend to see it as a process that takes places throughout the development of the embedded software and a way to force control engineers and programmers to keep in mind that algorithms should follow mathematical evidence and need to be implemented with caution and rigor.

Distributed systems and architectures are growing in interest within Aerospace Industries. It is especially the case for engine manufacturers. Formal verification or certification of an engine controller running on a multi-core architecture is very challenging.

The work developed in this thesis naturally fits the use of multi-core architecture in the sense that it is extremely promising in terms of computational resources and power, which are crucial features when using real-time optimization based control.

Work towards the parallelization of optimization algorithms and its formal verification would therefore be interesting and might need to be addressed.

u

  Collection of input vectors to horizon: u = [u 1 . . . u N -1 ]. xiv X Original Decision Vector for an Optimization problem.ZProjected Decision Vector for an Optimization problem.X fFeasible set of an optimization problem.XEpsilon optimal set of an optimization problem.ACSLrequires Introduces a precondition in a local contract or a function contract. ensures Introduces a postcondition in a local contract or a function contract. assigns Listing the memory being assigned in a function contract. \result Referring to the output of a C function. loop invariant Introduces a loop invariant, a property true at each iteration of the loop. loop variant Introduces a strictly positive and decreasing quantity in a denumerable set, usual N. Needed to prove termination. loop assigns Listing the memory being assigned in a loop. axiomatic Defining a new ACSL Theory. Containings new types, functions, axioms, lemmas, theorems. lemma Defining a new lemma. Need to be proven. theorem Defining a new theorem. Need to be proven. axiom Defining a new axiom. Assumed to be true. logic Defining a new ACSL function. predicate Defining and naming a given ACSL property. xv genemo Files

Figure 1 :

 1 Figure 1: V & V Cycle for Credible Autocoding

Figure 2

 2 Figure 2 gives an illustration of the corresponding toolchain.

Figure 2 :

 2 Figure 2: MPC Integration of Credible Autocoded Algorithms

Figure 3 :

 3 Figure 3: Hoare Triple

Figure 4 :

 4 Figure 4: Hoare Triple Expressed in ACSL Example 1

Figure 5 :

 5 Figure 5: Hoare Triple Expressed in ACSL Example 2

2 3 /* @ 4 @ 7 intFigure 6 : 2 3 /* @ 4 @ 5 @ 6 @ 10 *Figure 7 :Figure 8 :

 247624561078 Figure 6: Memory Assignment Issue Example in ACSL

Figure 9 : 9 . 1

 991 Figure 9: Weakest Precondition Illustration

ure 10 )Figure 10 :

 1010 Figure 10: Closed-Loop System

  and lim t→∞ x(t) = x e .

Figure 11 :

 11 Figure 11: Classification of Some Convex Optimization Problems

Figure 12 :Figure 13 :

 1213 Figure 12: MPC Closed-Loop System

Figure 14 :

 14 Figure 14: Spring-Mass System

Figure 15 :

 15 Figure 15: Closed-Loop Position Versus Time

Figure 16 :

 16 Figure 16: Closed-Loop Velocity Versus Time

Figure 17 :

 17 Figure 17: Input Force Versus Time

Figure 18 :

 18 Figure 18: Ellispoid Method Trade off

  Then, we define the ellipsoid E k+1 by the minimal volume ellipsoid containing the half ellipsoid Êk that is known to contain x * . The Figures19 and 20illustrate such ellipsoids cuts.

Figure 19 :

 19 Figure 19: Ellipsoid Cut

Figure 20 :

 20 Figure 20: Ellipsoid Cut In an LP Settings

Figure 21 :

 21 Figure 21: Included and Including Balls

ACSL 1 / 4 logic

 14 vector v e c _ o f _ 1 6_ s c a l a r ( double * x ) reads x [0..15]; 5 logic vector v e c _ o f _ 3 6_ s c a l a r ( double * x ) reads x [0..35]; 6 ...

7Figure 22 :

 22 Figure 22: ACSL Linear Algebra Theory

  Figure 22 represents an extract from the autocoded ACSL linear algebra theory (which can be found in appendix A).

Figure 23 :

 23 Figure 23: Ellipsoid Type Definition

Figure 24 :Figure 25 :

 2425 Figure 24: Existence vector ACSL Lemma

ACSL 1 Figure 26 :

 126 Figure 26: ACSL Optim Type Definition

  illustrates the axiomatization of a constraint calculation and the feasibility predicate definition. When instantiating a object of type vector or ma-ACSL 1 /* @ 2 axiom c o n s t r a i n t _ l i n e a r _ a x i o m : 3 \forall optim OPT , vector x , integer i ; 4 getm ( OPT )[ i ] == 0 == > 5 constraint ( OPT , x , i ) == 6 -scalarProduct ( getci ( OPT , i ) ,x , size_n ( OPT )) -getdi ( OPT , i ); 7 axiom c o n s t r a i n t _ s o c p _ a x i o m : 8 \forall optim OPT , vector x , integer i ; 9 getm ( OPT )[ i ] != 0 == > 10 constraint ( OPT , x , i ) == 11 twoNorm ( vector_affine ( getAi ( OPT , i ) ,x , getbi ( OPT , i ))) -12 scalarProduct ( getci ( OPT , i ) ,x , size_n ( OPT )) -getdi ( OPT , i ); 13 ...

14 predicate 15 isFeasibleFigure 27 :

 141527 Figure 27: ACSL Feasible Predicate Definition

Figure 28 :

 28 Figure 28: getNorm 2 Header C Code File

1Figure 29 :Figure 30 :

 2930 Figure 29: getNorm 2 Body C Code File

Figure 31 :

 31 Figure 31: getp.h Header File

4. 1 . 1

 11 US Patriot MissileDuring the Gulf war in Saudi Arabia, 1991, a US patriot missile failed to intercept an incoming Iraqi Scud missile. It ended up killing 28 US soldiers and injuring around 100 other people.

Figure 32 :

 32 Figure 32: US Patriot Missile

Figure 33 :

 33 Figure 33: European Ariane 5 Rocket

Figure 34 :

 34 Figure 34: Corrected Ellipsoid

4 5 6 7 /* @ 8 @ 9 @Figure 35 :

 468935 Figure 35: getFrobeniusNorm Annotated C Function

Figure 36 :

 36 Figure 36: Ellipsoid Widening

Figure 37 :Fact 1 [

 371 Figure 37: P u and P l Polyhedral Sets

Figure 38 :

 38 Figure 38: Unit balls of R 2

Figure 39 :

 39 Figure 39: Linear Relaxation of a Second Order Cone

1 Figure 41 :

 141 Figure 41: Objective Section For The Spring-Mass System

ure 43 ) 1 Figure 42 :

 43142 Figure 42: Constraint Section For The Spring-Mass System

Figure 43 :

 43 Figure 43: Information Section For The Spring-Mass System

5 #

 5 This is a comment

Figure 44 :

 44 Figure 44: Input/Output Sections For The Spring-Mass System

Figure 45 :

 45 Figure 45: Internal Aspect of The Autocoder

5. 4 . 2

 42 Figure 46: 3 DOF Helicopter

λ.

  Using double precision floating points and a accuracy of = 0.25 we have the following results: λ = 1.0063428 An original number of iteration equal to N = 5528 and therefore an updated number of steps of: N λ = 6817. Doing the simulation on a Intel Core i5-3450 CPU @ 3.10GHz × 4 processor we have a running time of approximately 0.2 sec for a single point and therefore it took the computer 4 sec to simulate 10 sec (because we discretized the system at 2Hz). The results of the simulation are presented in Figures 49 and 50. We used the text file presented in Figure 68 in appendix C to generate the C code in order to performed the simulation. The whole simulation was performed using the generated C code and the Simulink model presented in Figure 48. The block "MPC CONTROLLER", with the help of the Matlab coder, compiles and calls the generated and annotated C code.

Figure 48 :

 48 Figure 48: Simulink File Used for the Quanser Simulation

Figure 50 :

 50 Figure 50: Simulation of Lowest Altitude versus time

Figure 51 :

 51 Figure 51: Cone constraint on drone

Figure 52 :

 52 Figure 52: MPC Simulation for a Quadcopter

  > pow (a , n ) == a * pow (a ,n -1); axiom pow_mult : \forall real a ,b , integer n ; pow (a , n ) * pow (b , n ) == pow ( a *b , n ); axiom pow_r_eps_N : \forall real epsilon , V , r ; r == 0.2 == > epsilon == 0.1 == > V == 110 == > pow ( epsilon / V * r , 3) >= 6.0105004 e -12; axiom n o n _ s i n g u l a r _ m a t _ m u l t _ s c a l a r : \forall matrix A , real a ; invertible ( A ) == 1 == > a != 0 == > invertible ( m a t_ mu lt _ sc al ar (A , a )) == 1; axiom n o n _ s i n g u l a t I _ A B : \forall vector x , real a ; twoNorm ( x ) == 1 == > a != -1 == > invertible ( mat_add ( ident (3) , m at _m u lt _s ca l ar ( ketbra ( x ) , a ))) == 1; axiom det_non_null : \forall matrix A ; invertible ( A ) == 1 == > det ( A ) != 0; axiom de t_ m at ri x_ m ul t : \forall matrix A , B ; det ( mat_mult (A , B )) == det ( A )* det ( B ); axiom d e t _ m a t r i x _ m u l t _ s c a l a r : \forall matrix A , real a ; mat_col ( A ) == 3 == > mat_row ( A ) == 3 == > det ( ma t_ m ul t_ sc a la r (A , a )) == pow (a ,3)* det ( A ); axiom det_I_AB : \forall vector x , real a ; twoNorm ( x ) == 1.0 == > det ( mat_add ( ident (3) , ma t _m ul t_ s ca la r ( ketbra ( x ) , a ))) == 1.0 + a ; axiom n o r m a l i z e _ v e c t o r : \forall vector x ; twoNorm ( v ec _ mu lt _s c al ar (x ,1/ twoNorm ( x ))) == 1; axiom k e t b r a _ a x i o m _ s e l e c t : \forall vector p , integer i , j ; 0 <= i < 3 == > 0 <= j < 3 == > mat_select ( ketbra ( p ) ,i , j ) == vector_select (p , i ) * vector_select (p , j ); axiom k e t b r a _ a x i o m _ s i z e : \forall vector p ; mat_row ( ketbra ( p )) == vector_length ( p ) && mat_col ( ketbra ( p )) == vector_length ( p ); axiom a b s o l u t e v a l u e _ p o s i t i v i t y : \forall real a ; absolutevalue ( a ) >= 0; axiom a b s o l u t e v a l u e _ p o s i t i v e : \forall real a ; a >= 0 == > absolutevalue ( a ) == a ; axiom a b s o l u t e v a l u e _ n e g a t i v e : \forall real a ; a <= 0 == > absolutevalue ( a ) == -a ; axiom t r a n s p o s e _ s e l e c t : \forall matrix A , integer i , j ; 0 <= i < mat_col ( A ) == > 0 <= j < mat_row ( A ) == > mat_select ( transpose ( A ) , i , j ) == mat_select (A , j , i ); ( A )) == mat_row ( A );axiom v e c _ m u l t _ s c a l a r _ l e n g t h : \forall vector A , real x ; vector_length ( v ec _mul t _s ca la r (A , x )) == vector_length ( A );

Figure 54 :Figure 55 :Figure 56 :Figure 57 :Figure 58 :Figure 62 :Figure 64 :

 54555657586264 Figure 54: Linear Algebra ACSL Theory (Part 2)

Table 1 :

 1 2 • 10 -5 to apply in order to account for the rounding errors. As stated in Section 4.3, a new number of iterations of N λ = 1563. is now needed. In order to automatically generate the C code used for the control, we used to GENEMO input file presented in Appendix C.3. We repeated the same analysis for different horizon values. The results are collected in Tables1 and 2. Performances For Different Problem Sizes

	N	r	R	V	n z N it f * (Ell) RunTime(s) f * (CV X)
	5 4.79 12.24 70.71 4 300 5.1565	0.001	5.15648
	7 4.69 14.90 98.99 6 677 6.9293	0.003	6.92928
	10 4.54 18.19 141.42 9 1556 9.3116	0.010	9.31154
	15 4.30 22.67 212.13 14 3916 12.6910	0.050	12.6910
	20 4.05 26.47 282.84 19 7466 15.6975	0.169	15.6974
				N N it λ (double)	N λ (double)
				5 300 1 + 5.9e -7	301
				7 677 1 + 5.2e -6	678
				10 1556 1 + 5.2e -5	1563
				15 3916 1.000709	4548
				20 7466 1.004669 Not Converging

Table 2 :

 2 Floating Point Consideration

Table 3 :

 3 Performances For different Horizon -Drone of the MPC simulation applied for the quadricopter using the autocoded C code is shown in figure

	. A plot

ACKNOWLEDGEMENTS

I would like to thank everyone who made this thesis possible and an unforgettable adventure for me. First, I would like to thank my committee members, Professor

As usual, none of this would have been possible without funding, therefore I want to thank the sponsor, the National Science Foundation (NSF) and the French Aerospace Laboratory (ONERA).

Using the type of floating-point transformation [START_REF] Rtca | DO-178C, software considerations in airborne systems and equipment certification[END_REF] and neglecting all terms in eta and powers of u greater than two we get:

Error on B + . Similarly, For each component of B, we have: Hence, we found:

4.4 Necessary conditions for numerical stability

Problem Formulation

In order to take into account the uncertainties on the variables due to floating-point rounding, we want to modify the algorithm to make it more robust. For this, we choose to evaluate those uncertainties and conclude on a coefficient λ that represents by how much we are going to widen the ellipsoid E k at each iteration (see Figure 36).

Let us assume we have B ∈ F n×n , p ∈ F n , c ∈ F n . We want to find λ ≥ 1 ∈ R such that:

Equivalent and Sufficient Conditions for Covering

In this section, we state two lemmas that give an equivalent and a sufficient condition for the ellipsoid Ell λ • fl(B + ), fl(c + ) to include the ellipsoid Ell B + , c + .

We now state lemma 2, which gives a sufficient condition for the coefficient λ to have Ell λ • fl(B), fl(c) including Ell(B, c). If the calculations of B + and c + were perfect, using lemma 2, we would have λ = 1 working; no correction is indeed necessary.

Lemma 2 [Widening -Sufficient Condition]

Proof let us assume that

The goal here is to end up this the inclusion property [START_REF] Patrinos | A dual gradientprojection algorithm for model predictive control in fixed-point arithmetic[END_REF]. Using the fact that the two norm is a consistent norm, we have:

and:

Therefore using equations (41), ( 42) and (40) we have:

Then, using the triangle Inequality ( x + y ≤ x + y ), we finally have that:

So, using lemma 1 we have:

Using equation (57):

So, if:

Using lemma 2 we conclude then on the inclusion property [START_REF] Patrinos | A dual gradientprojection algorithm for model predictive control in fixed-point arithmetic[END_REF].

Thus, due to this lemma, following the floating-point analysis of the algorithm, we are now able to compute a coefficient λ such that equation ( 39) is valid. After founding such a λ, we would like to know whether the algorithm is still converging. On the other hand, because the method's proof lies in the fact that the final ellipsoid has a small enough volume, this correction will have an impact of the number of iterations.

Lemma 4 addresses those issues.

The algorithm implementing the widened ellipsoids, with coefficient λ, will converge if:

+ number of rows of S). Then Z ∈ R nz with n z = n x -d. Thus, the original optimization problem 48 is equivalent to the below projected problem.

minimize

With:

More details and references about equality constraints elimination can be found at [START_REF] Boyd | Convex optimization[END_REF]. In order to solve this optimization problem with the Ellipsoid Method and have convergence guarantees for every possible x o , we need to compute geometric characteristics on the feasible set of this parametric optimization problem. For that reason, for now, let us assume:

We define the parameterized polyhedral set:

Having this collection of polyhedral sets, we want to compute ball radii that will tell us about the volume of the feasible set that would be true for every initialization point x o . Beforehand, let us define the operator φ(.) that returns for a matrix A the vector φ(A) whose coordinates are the two norm of the rows of A. Let us now consider the two extreme polyhedral sets below:

In order to give an example of this concept, we show in Figure 37 an illustration of such polyhedral sets (the illustrated sets have no physical meaning and do not represent any MPC problem. We used the values below: time for a given problem size and number of iterations. In this thesis, we only highlight that it is a needed work in order to predict feasibility of the real-time use of those algorithms but do not present such analysis.

GENEMO Programming Language Syntax

The user can access the norm of a vector x by writing ||x||. On the other hand, if

x is a matrix, x(:, k) denotes the k-th column of the matrix x. A ,b ,C ,d ,f ,m )) == 2; axiom equalityOpt : \forall matrix A1 , A2 , C1 , C2 , vector b1 , b2 , d1 , d2 , f1 , f2 , int * m1 ,* m2 ; \forall optim OPT , vector x , real epsilon , V ;

in ( e p s i l o n _ o p t i ma l _ s e t ( OPT , epsilon , V ) , x ) == > cost ( OPT , x ) <= epsilon + minimum ( OPT ); lemma e p s i l o n _ o p t i m a l _ s e t _ c o n s t r a i n t s : \forall optim OPT , vector x , real epsilon , V ; 0 < epsilon / V < 1 == > V > 0 == > in ( e p s i l o n _ o p t i ma l _ s e t ( OPT , epsilon , V ) , x ) == > ( \exists vector y ; isNegatif ( constraints ( OPT , y )) && isNegatif ( constraints ( OPT , Optimal ( OPT ))) &&

x == vector_add ( vec_ mu l t_ sc al a r (y , epsilon / V ) , vec _mul t_s ca la r ( Optimal ( OPT ) , 1 -epsilon / V )) && isNegatif ( constraints ( OPT , x )) ); lemma e p s i l o n _ o p t i m a l _ s e t _ c o n s t r a i n t s 2 : \forall optim OPT , vector x , real epsilon , V ;

in ( e p s i l o n _ o p t i ma l _ s e t ( OPT , epsilon , V ) , x ) == > isFeasible ( OPT , x ); lemma e p s i lo n _ o p t i m a l _ s e t : \forall optim OPT , vector x , real epsilon , V ;

in ( e p s i l o n _ o p t i ma l _ s e t ( OPT , epsilon , V ) , x ) == > i s Ep s il o nS ol u ti o n ( OPT , x , epsilon ); lemma b e t t e r _ e p s i l o n _ o p t i m a l _ s e t : \forall optim OPT , vector y , x , real epsilon , V ;