2.1 Isabelle/HOL

Isabelle is an interactive higher-order-logic (HOL) theorem prover, a kind of proof assistant, based on automated reasoning techniques and logical calculus. It is written using ML [Reference Milner37], a functional programming language, in which rules are presented as propositions and the proofs are structured by combining rules using the $\lambda$ -calculus. Isabelle provides the ability to express mathematical theorems of control engineering in formulae of formal logic and prove them using the automated derivation tools provided. When a proof is completed in Isabelle, one can be certain that the mathematical theorem at hand is valid.

Isabelle/HOL is one of Isabelle’s platforms with quantifiers over domains of variables and formal semantics. Isabelle/HOL has a proof-language called Isar, which structures the proofs in a human readable and understandable mathematical formal text for both users and computers. The mathematical formulae can be formalised and proven in Isar with the aid of Isabelle’s logical tools.

Examples of such tools are the simplifier, which performs operation and reasoning on equations; the classical reasoner that carries out long chains of reasoning procedures; and algebraic decision procedures using advanced pattern matching by the package sledgehammer for automatically finding the proofs based on already-proven theorems in Isabelle’s library. The latter also calls external FOL (first-order-logic) provers (ATPs) and Satisfiability Modulo Theories (SMT) solvers [Reference Barrett and Tinelli38]. Isabelle has been chosen in this work for its rich logical and automation techniques since it has a large library, which contains most of the mathematical theories that are useful in control systems theory proofs.

The use of Isabelle is justified as advanced robust and adaptive control theory heavily uses concepts of convergence in limits, differentiation and integration, complex analysis, abstract vector spaces, spaces of functionals and nonlinear operators, stability theory, and formal measures of control performance. Isabelle also supports other logical systems such as Hoare logics for systems verification [Reference Nipkow39] in addition to the ability to generate executable codes of the proven statements. There is a wide range of syntax in Isabelle/HOL, the most common and useful Isabelle/HOL expressions and symbols are summarised in Table 1.

Table 1. Isabelle/HOL symbols and expressions

2.2 MetiTarski

MetiTarski is a FOL automated theorem prover, which is designed as a combination of the resolution prover Metis [Reference Hurd40] and a decision procedure tool QEPCAD [Reference Brown41] to conduct proofs of theorems stated algebraically over the field of real numbers. It can be utilised to prove universally quantified inequalities of transcendental and other functions as well, for instance ln, log, exp, etc. MetiTarski reduces problems to decidable ones by replacing some functions with their upper and lower bounds. It can be made to perform proofs automatically with the aid of three reasoners: QEPCAD, Z3 [Reference De Moura and Bjørner42] and Mathematica. It can, for instance, be used to check the stability of an unmanned aerial vehicle (UAV) during flight, due to its ability to prove inequalities, which occur in the theory of dynamical stability, at a reasonable speed.

4.1 Verification using the Isabelle/HOL prover

The attitude controller of a generic multi-rotor UAV, as proposed in [Reference Jasim and Veres31], is considered here. The aircraft’s attitude dynamics are controlled by a robust nonlinear controller that takes into account the dynamical modelling uncertainties and payload disturbances. Figure 2 illustrates the attitude controller.

Figure 2. The inner control loop of the multi-rotor controller [Reference Jasim and Veres31] for attitude control. q is the measured and $q_r$ is the desired attitude, $q_e$ is its error in quaternions, $\xi$ is angular error vector, $\omega, \omega_r$ are the angular velocity and reference, $K_q, K_\omega$ are control gains, $I,\tilde{I}$ are the inertia matrix and its estimate. $\Gamma(\omega)$ is a cross-product matrix for the Coriolis forces such that $\Gamma(\omega) I \omega = \omega \times I \omega$ is satisfied and $u_d$ is torque disturbance of the attitude.

Simulation cannot be used to fully prove that the UAV’s control system is robust for all possible flight conditions within its flight envelope. To ensure correctness of the controller robustness stated in [Reference Jasim and Veres31], the controller design’s derivations and stability analysis are to be verified using the Isabelle/HOL prover. The Isabelle/HOL software system is applied for this purpose due to its rich library of mathematical theorems, which are required to perform the UAV’s control system verification.

The verification process using Isabelle/HOL is illustrated in Fig. 3, which consists of two parts: formalisation (HOL Platform) and proof procedures (Proof Tools). The first part starts by formalising the multi-rotor UAV system into the Isabelle/HOL syntax such as the coordinate system, rotational dynamics, time-domain functions, proposed assumptions and the aircraft’s stability analysis. The implementation of the control design and aircraft dynamics includes a series of Isabelle items such as definition, lemma and theorem.

Figure 3. Formalising and proving UAV’s controller in Isabelle/HOL theorem prover.

In this demo some new lemmas were needed to be stated and proven, which did not yet exist in Isabelle as that library is still under development. The formalisation also required importing some proven mathematical theories and lemmas from the prover’s library, which were used in formalising the control system equations with the controller’s assumptions and definitions.

The theories implemented under an HOL formalism will be described here to illustrate the proof procedures. The $HOL.thy$ is at the core of HOL formalism, which includes definitions of real numbers (real.thy), functions (fun.thy), sets (set.thy), etc. These are necessary in most formal procedures. The $Quaternions.thy$ can be used for quaternion definitions, operations and also to represent aircraft attitude. A multi-variable analysis package $Analysis.thy$ can be exploited for function operations over real numbers. The Finite_Cartesian_Product.thy and Inner_Product.thy can be utilised for definitions and operations over real vectors, L2_norm.thy and Norm_Arth.thy are used for real vector norms and their operations, respectively.

A time sub-domain has been defined as $T= \{t | \ t \in [0,\infty\}\}$ and relied on in definitions of sets of vector signals. For instance, the concept of a function of time (i.e. signal) over $[0,\infty)$ and its mapping to another signal is an operator. Formal nonlinear operators have been defined within Isabelle to express unknown dynamics, also used in a previous work [Reference Jasim and Veres1], to map signals to signals by dynamical operators. These nonlinear operators are very important in descriptions and proofs of controller robustness. It is worth mentioning that in the Isabelle/HOL it was possible to define such operations using the “locale” syntax [Reference Nipkow, Paulson and Wenzel30].

The mentioned theory $Analysis.thy$ includes definitions of real vectors (Finite_Cartesian_Product.thy), vector norms (L2_norm.thy, Norm_Arth.thy) and their operations. These theories enable the definitions of the aircraft’s three-dimensional rotation vectors and their norms including the torque, angular velocity and acceleration vectors, where each component of a vector is represented as a continuous time-domain function for signals. The continuous functions defined in $Fun.thy$ , $Function\_Algebras.thy$ and $Topological\_Spaces.thy$ theories have been heavily relied on.

Matrices are formalised in $Matrix.thy$ , and their operations and derivations performed in $Analysis.thy$ , $Finite\_Cartesian\_Product.thy$ and $Real\_$ $Vector\_Spaces.thy$ . The rate change of the multi-rotor’s attitudes, i.e. velocities and accelerations, are formalised by time derivatives in the theories $Deriv.thy$ and $derivative.thy$ .

Multi-rotor controller design requires several robustness assumptions, which need inequalities over the field of real-numbers. Fortunately, such inequalities have been defined in the Isabelle prover in the theory $Orderings.thy$ . This is an important feature for any robust control design.

The second part “Theorem Proving” in Fig. 1 is an interactive process between the designer/engineer and the automated proving tools that the Isabelle prover supports. The role of designer/engineer is to exploit the ITP step-by-step the achieve a proof, where the system is not able to find a proof automatically, by simplifying the statement into several steps. Each step should be proven using the provided automated tools before moving to the next one, otherwise the prover will not pass the statement. Examples of the automated tools supported in Isabelle are: CVC4, Z3, SPASS, E prover, Remote_Vampire and SMT solvers. In addition, Isabelle has its own automatic proving tools such as auto, simp, blast, etc. Most of the control systems verified in this work required an interaction with the prover due to the design complexity making it not possible for the prover to solve them automatically.

The Isabelle code used in the case study is too long to be stated here, only the relevant definitions and proofs are shown, while the complete code can be found in an online repository [Reference Jasim and Veres43]. Formalisation and proofs are illustrated in the remainder of this section.

The multi-rotor attitude dynamics in [Reference Jasim and Veres31],

(1) \begin{equation} I \dot{\boldsymbol{\omega}} + \Gamma(\boldsymbol{\omega}) I \boldsymbol{\omega} + \boldsymbol{\tau}_d = \boldsymbol{\tau},\end{equation}

where $I \in \Re^{3 \times 3}$ is a symmetric and positive-definite inertia matrix of the craft about its mass centre. $\boldsymbol{\tau}_d(t)=[\tau_{d \phi}(t) \ \tau_{d \theta}(t) \ \tau_{d \psi}(t)]^{T} \in \Re^{3}$ are unknown disturbances torques with $\phi,\theta$ and $ \psi$ being the roll, pitch and yaw angles, respectively. $\boldsymbol{\tau}(t) = [\tau_{\phi}(t) \ \tau_{\theta}(t) \ \tau_{\psi}(t) ]^{T} \in \Re^{3} $ is the torque vector of the onboard controller in the body frame which produces the multi-rotor motion. These can be formalised by the following code in Isabelle/HOL:

and the torque vector $\boldsymbol{\tau}$ in [Reference Jasim and Veres31],

(2) \begin{equation} \boldsymbol{\tau} =\begin{bmatrix}\tau_{\phi} \\ \tau_{\theta} \\ \tau_{\psi}\end{bmatrix}=\begin{bmatrix}\ell l(\Omega^{2}_{2} - \Omega^{2}_{4}) \\\ell l(-\Omega^{2}_{1} + \Omega^{2}_{3}) \\b(-\Omega^{2}_{1} + \Omega^{2}_{2} - \Omega^{2}_{3} + \Omega^{2}_{4})\end{bmatrix},\end{equation}

is defined by bounding all rotor angular velocities $\Omega_i$ with their maximum value $\Omega_{\max}$ as

where $\ell$ is the length from the centre of mass of the multi-rotor to the rotor, and l and b are the aerodynamic force and torque constants of the rotors. The control input $\textbf{u}$ [Reference Jasim and Veres31],

(3) \begin{equation}\textbf{u} = \dot{\boldsymbol{\omega}}_{r} + K_{\boldsymbol{\omega}}\dot{\boldsymbol{\xi}} + K_{\textbf{q}}\boldsymbol{\xi},\end{equation}

and the control torque [Reference Jasim and Veres31],

(4) \begin{equation}\boldsymbol{\tau} = \hat{I} \textbf{u} + \textbf{u}_d + \hat{\textbf{c}}(\boldsymbol{\omega}),\end{equation}

are defined in the prover through the following code:

where $\hat{I}$ is an estimated matrix of the inertia matrix I of the craft, $\textbf{u}$ represents a new input vector to be designed later on in Equation (3), $\boldsymbol{\omega}(t)$ is the angular velocities vector, $\dot{\boldsymbol{\omega}}(t)$ and $\dot{\boldsymbol{\omega}}_{r}(t)$ are the angular acceleration and reference angular acceleration vectors, respectively. $K_{\omega}=diag[k_{\omega_{1}} k_{\omega_{2}} k_{\omega_{3}}] \in \Re^{3 \times 3}$ is a positive-definite diagonal gain matrix setting the error gains in feedback and $K_{q}=diag[k_{q_{1}} k_{q_{2}} k_{q_{3}}]\in \Re^{3 \times 3}$ is a positive-definite diagonal gain matrix. $[q_{e_{1}} \ q_{e_{2}} \ q_{e_{3}}]^{T}$ is a term defined to reduce the the dimensions of quaternions. $\hat{\textbf{c}}(\boldsymbol{\omega})$ is an estimate of $\textbf{c}(\boldsymbol{\omega})= \Gamma(\boldsymbol{\omega}) I \boldsymbol{\omega}$ as based on $\hat{I}$ and measured $\boldsymbol{\omega}$ . The additional term $\textbf{u}_d$ is added to render the effects of uncertainty and disturbances in order to guarantee robustness against these effects; $\textbf{u}_d$ will be defined later, to counter these, in Equation (14). The derivative of $\dot{\boldsymbol{\omega}}$ [Reference Jasim and Veres31] is

(5) \begin{equation} \begin{split}\dot{\boldsymbol{\omega}} &= I^{-1}\hat{I} \textbf{u} + I^{-1} \textbf{u}_d + I^{-1} [\boldsymbol{\Delta} (\boldsymbol{\omega}) - \boldsymbol{\tau}_d] \\&= \textbf{u} + (I^{-1} \hat{I} - \mathbb{I}) \textbf{u} + I^{-1} \textbf{u}_d + I^{-1}[\boldsymbol{\Delta} (\boldsymbol{\omega}) - \boldsymbol{\tau}_d] \\& = \textbf{u} + I^{-1} \textbf{u}_d - \textbf{y} \\\end{split}\end{equation}

where

(6) \begin{equation} \begin{split}\textbf{y} &=[\mathbb{I} - I^{-1}\hat{I}] \textbf{u} -I^{-1}[\boldsymbol{\Delta} (\boldsymbol{\omega}) - \boldsymbol{\tau}_d] \ \ , \ \ \boldsymbol{\Delta}(\boldsymbol{\omega})=\hat{\textbf{c}}(\boldsymbol{\omega})- \textbf{c}(\boldsymbol{\omega})\end{split}\end{equation}

is formalised in Isabelle/HOL based on the $att\_dyms$ , $cont\_u$ and $cont\_law$ (see the the repository [Reference Jasim and Veres43] for the details). The closed-loop error dynamics [Reference Jasim and Veres31] is

(7) \begin{equation} \begin{split}\dot{\boldsymbol{\eta}} &= A\boldsymbol{\eta} + G[\textbf{y} - I^{-1} \textbf{u}_d]\end{split}\end{equation}

where $\boldsymbol{\eta} = [\boldsymbol{\xi}^T \ \ \dot{\boldsymbol{\xi}}^T]^{T} \in \Re^{6 \times 1}$ and

(8) \begin{equation} A =\begin{bmatrix}0^{3 \times 3}\;\;\; & \mathbb{I}^{3 \times 3} \\-K_{q}^{3 \times 3}\;\;\; & -K_{\omega}^{3 \times 3}\end{bmatrix}, \ G =\begin{bmatrix}0^{3 \times 3} \\ \mathbb{I}^{3 \times 3}\end{bmatrix},\end{equation}

are implemented as

The $set\_of\_definitions$ in the above code are definitions, which rely on many available, pre-defined definitions in Isabelle. The flight envelope assumptions proposed in [Reference Jasim and Veres31] are as follows.

Assumption 1,

(9) \begin{equation}\sup (\| \dot{\boldsymbol{\omega}}_{r}\|) < \alpha.\end{equation}

Assumption 2,

(10) \begin{equation}\lambda_{\min} \le \| I^{-1}\| \le \lambda_{\max}.\end{equation}

(11) \begin{equation} \| \mathbb{I} - I^{-1}\hat{I}\| \le \delta \le 1,\end{equation}

where $\lambda_{\min}>0$ and $\lambda_{\max}>0$ a lower and upper bound, respectively; $\delta>0$ is a constant.

Assumption 3 ,

(12) \begin{equation}\| \boldsymbol{\tau}_d \| \le \gamma,\end{equation}

where $\gamma>0$ is a known upper constant bound. The above assumptions are formalised in Isabelle/HOL as follows:

Stability analysis of the attitude controller as stated in Equations (38)–(47) [Reference Jasim and Veres31] is implemented in Isabelle/HOL using a set of definitions in (definition), several lemmas in (lemma) and short theorems in terms of (theorem). This structure of using several lemmas and theorems during the proof is due to the fact that the reasoning system of the theorem prover cannot handle long proofs with many assumptions, i.e. the Isabelle system is unable to prove statements, which have many equations, if they are formalised in only one lemma or theorem style. However, stability analysis starts by defining the candidate Lyapunov function V [Reference Jasim and Veres31],

(13) \begin{equation} V(\boldsymbol{\eta}) = \boldsymbol{\eta}^{T} Q \boldsymbol{\eta} >0 \ , \ \forall \boldsymbol{\eta} \ne 0\end{equation}

which is formalised as a definition in Isabelle/HOL:

where $Q \in \Re^{6 \times 6}$ is a symmetric positive-definite matrix. Taking the candidate Lyapunov function V, the time derivative of Lyapunov function is derived and the derivations in Equations (39)–(41) in [Reference Jasim and Veres31] are proven symbolically and detailed in the online repository. An outline is shown below:

The term $\textbf{u}_\textbf{d}$ [Reference Jasim and Veres31],

(14) \begin{equation}\textbf{u}_d =\begin{cases}\dfrac{ \zeta(\boldsymbol{\eta},t)}{\|G^{T}Q \boldsymbol{\eta}\|} G^{T}Q \boldsymbol{\eta}, &\ \ \ if \ \ \|G^{T}Q \boldsymbol{\eta}\| \ge \mu\\\dfrac{ \zeta(\boldsymbol{\eta},t)}{\mu} G^{T}Q \boldsymbol{\eta},&\ \ \ if \ \ \|G^{T}Q \boldsymbol{\eta}\| < \mu.\end{cases}\end{equation}

is defined then as the derivative in Equation (43) [Reference Jasim and Veres31]. This is performed using the Cauchy-Schwartz inequality (see “theorem $Eq\_43$ ” in the repository). Based on the $\zeta(\boldsymbol{\eta},t)$ definition, $\varepsilon > 0$ , $\lambda_{\min} >0$ ,

(15) \begin{equation} \zeta(\boldsymbol{\eta},t) \ge \dfrac{\varepsilon}{\lambda_{\min}},\end{equation}

and the upper bound of the norm of $\textbf{y}$ derived in Equation (45) (31) (see “theorem $Eq\_45$ ” in the repository), $\zeta(\boldsymbol{\eta},t)$ in (15) is obtained (see “theorem $Eq\_46$ ” in the repository). The terms $\textbf{u}_\textbf{d}$ and $\zeta(\boldsymbol{\eta},t)$ are implemented in the prover as “ $u_{d}\_def$ ” and “ $zeta\_def$ ” respectively:

Note that the short arrow $\rightarrow$ in the code refers to implies, while the longer $ \longrightarrow$ refers to convergence in HOL.

Finally, based on all the above definitions and assumptions, it has been verified that the proposed control system is asymptotically stable since the time derivative of Lyapunov function [Reference Jasim and Veres31],

(16) \begin{equation}\dot{V}(\boldsymbol{\eta}) = - \boldsymbol{\eta}^{T} P \boldsymbol{\eta} + 2 \boldsymbol{\eta}^{T}QG (\textbf{y}- I^{-1} \dfrac{\zeta(\boldsymbol{\eta},t)}{\|G^{T}Q \boldsymbol{\eta}\|} G^{T}Q \boldsymbol{\eta}) <0,\end{equation}

(17) \begin{equation}\dot{V}(\boldsymbol{\eta}) = - \boldsymbol{\eta}^{T} P \boldsymbol{\eta} + 2 \boldsymbol{\eta}^{T}Q G(\textbf{y}- I^{-1} \dfrac{ \zeta(\boldsymbol{\eta},t)}{\mu} G^{T}Q \boldsymbol{\eta} ) <0.\end{equation}

is strictly negative for $\forall \boldsymbol{\eta} \ne 0$ where $A^{T}Q + QA = -P$ . It has also been proven that the tracking error converges to zero as the time converges to infinity ( $\|\boldsymbol{\eta}\| \longrightarrow 0$ ). The code below illustrates the symbolic proof in the Isabelle theorem prover.

Figure 4. Onboard verification framework of UAVs.

4.2 Onboard verification for a safe flight using MetiTarski

The control system of the UAV can be designed, simulated and verified at the modelling/design stage. The verified controller can then be converted into code and implemented on the autopilot system, which controls the aircraft trajectory. The aircraft can be exposed to gusts of wind, and the payload may also vary. Both may cause, in the extremes, unstable flight when the flight envelope is violated. The autopilot should provide a warning when the aircraft has entered an unstable region of its dynamical states as that may cause a crash or damage. This monitoring of the occurrence of instability can be implemented in an ATP tool, for instance the MetiTarski prover. The MetiTarski ATP has been chosen here to verify the controller stability of the aircraft due to its ability to deal with algebraic inequalities over numerical real numbers quickly. Unlike the previous verification stage using Isabelle, MetiTarski proves inequalities automatically, without the need for any interaction with a pilot or engineer.

MetiTarski can be implemented on the autopilot’s electronic board such as PixHawk [44] or Navio2 [45] Raspberry Pi and many others. Interaction between the autopilot and MetiTarski is easy to create in practice. A possible system architecture is illustrated in Fig. 4. The applicability of this approach is demonstrated here in simulation by interfacing the Simulink/Matlab model with MetiTarski and testing the stability of the control system. The multi-rotor model used was based on the realistic MathWorks Simulink model from [Reference Horton46].

For the stability analysis stated in Equations (16) and (17), the time derivative of the Lyapunov function needs to be tested to check whether it is negative definite or not. If it is not negative definite, then this indicates that the control system is out of its stability region, hence the autopilot can pass a warning message to the pilot to take some suitable action or perform emergency landing.

The verification process starts by formalising the stability Equations (16) and (17) into a FOL syntax. The parameters of these equations are then passed from the autopilot system to the MetiTarski prover via an interface. The test is conducted in the MetiTarski prover to test the condition ( $\dot{V}(\boldsymbol{\eta}) < 0$ ). The stability equations, Equations (16) and (17), were simplified using symbolic computations in Matlab before formalising them into FOL syntax. The parameters included in the stability equations were passed from Simulink/Matlab to MetiTarski to perform the monitoring. The following code illustrates stability checking in the MetiTarski prover as based on Equation (16) (the full code can be found in the online repository [Reference Jasim and Veres43]:

and for Equation (17),

The system entering an unstable region has been demonstrated in Simulink- Matlab by applying an external force to the nonlinear multi-rotor’s dynamics to simulate a strong gust of wind. The control system tried to overcome this dynamical change and it failed to keep the system stable because the external force was high and continuous. Without stability monitoring, this could have lead to the craft crashing.

The approach has been implemented not only to check Lyapunov stability but also to check if the limits of the system’s robustness parameters have been exceeded. This has been achieved by checking whether the maximum value of one or all of the rotor angular velocities $0<\Omega_i<\Omega_{\max}$ have been exceeded. This can be known from checking the torques for $\|\tau\| < \tau_{\max}$ , where $\tau_{\max}$ can be obtained from the known maximum angular velocity of each propeller $\Omega_{\max}$ and Equation (2). In the simulation this meant that the craft was unable to compensate against the external wind and caused the craft unable to follow its prescribed flight path. For more details, see our previous work [Reference Jasim and Veres47]. In such a situation the pilot, or autonomous agent, may decide not to insist on the planned flight path and make the most of the developing situation to save the craft by replanning.

All that has been demonstrated in simulation can also be implemented on a real UAV. This can be achieved as follows. First, implement MetiTarski on the autopilot’s electronic board with the stability and verification framework. Second, implement the interface between the autopilot and MetiTarski prover. Finally, implement the emergency decision framework inside the autopilot code as in Fig. 4.

The steps of a single cycle of stability testing, in a Simulink free onboard environment, can be implemented as follows:

• An interface is created for the autopilot which enables communication with the MetiTarski prover.

• The required parameters are communicated by the autopilot to the MetiTarski.

• The interface formalises the stability condition in FOL for MetiTarski.

• MetiTarski tests stability and communicates the result to the autopilot through the interface.

• The autopilot informs the human pilot or the onboard decision making software.

Using these methods the autopilot can send warning messages to the pilot, or in case of an autonomous mission, to an onboard intelligent software agent, to decide on a suitable counter manoeuvre [Reference Lusk, Glaab, Glaab and Beard48].