Experimental Validation of a Guaranteed Nonlinear Model Predictive Control

Abstract

This paper combines the interval analysis tools with the nonlinear model predictive control (NMPC). The NMPC strategy is formulated based on an uncertain dynamic model expressed as nonlinear ordinary differential equations (ODEs). All the dynamic parameters are identified in a guaranteed way considering the various uncertainties on the embedded sensors and the system’s design. The NMPC problem is solved at each time step using validated simulation and interval analysis methods to compute the optimal and safe control inputs over a finite prediction horizon. This approach considers several constraints which are crucial for the system’s safety and stability, namely the state and the control limits. The proposed controller consists of two steps: filtering and branching procedures enabling to find the input intervals that fulfill the state constraints and ensure the convergence to the reference set. Then, the optimization procedure allows for computing the optimal and punctual control input that must be sent to the system’s actuators for the pendulum stabilization. The validated NMPC capabilities are illustrated through several simulations under the DynIbex library and experiments using an inverted pendulum.

1. Introduction

Despite their maturing technology, robotic solutions still face many challenges related mainly to their design and control techniques. In practical applications, control of automated devices is often formulated depending on the nature of the operational task for which they are designed. Several research works have been realized in the past few years to develop optimal and reliable controllers, ensuring various reliability goals. For instance, authors in [1] investigated enhancing the reliability of control systems by using multiple controllers for a given plant, which is a delicate way to improve the reliability in control structures subject to both the controllers and plant failures. In [2], a design approach of reliable centralized and decentralized controllers is proposed, such that the resulting control systems remain stable and accurate despite sensor or actuator disturbances. Moreover, reliable H-infinity controllers for linear systems were designed in [3], in such a way that they ensure the asymptotic stability of the control systems, mainly when perturbations occur in some control components (sensors and actuators). These classical control strategies are generally designed based on dynamic and/or kinematic models and simplified control schemes. However, the drawback of these deterministic approaches is the lack of techniques to prove the robustness analysis of such controllers to different uncertain parameters in the plant.

Significant advances in stabilizing uncertain systems with parametric uncertainties have been recently treated in a very different way. Among the techniques used to achieve this end, interval analysis methods (IA) can provide guaranteed solutions to deal with complex nonlinear problems and to manage various uncertainties involved in the hardware and software parts of the automated device [4]. Indeed, IA tools turn out to be more promising and reliable to handle all the noises and uncertainties in the experimental and theoretical data, contrary to these classical approaches that are generally based on simplified models with deterministic values. IA methods are extremely popular in several research fields due to their efficiency and reliability in order to achieve a robust control in the presence of model uncertainties [5,6] or for state estimation [7]. For example, a robust controller based on interval arithmetic for constrained discrete-time nonlinear systems with additive uncertainties is presented in [5]. In addition, authors in [6] developed a robust feedback controller synthesis such that all the unknown parameters are assumed to vary between known bounds. This controller guarantees at each time step that the closed-loop system behavior is greater than the lower bound of the set-point interval and is lower than its upper bound.

Model-based Predictive Control (MPC) methodology (also known as the receding horizon control or the moving horizon control) has become a powerful control technique widely used in a large number of areas (e.g., the automotive industry and chemical processes [8]). There are several reasons for this: (i) its capability to anticipate the future on the reference trajectory (i.e., future changes in set-point), and (ii) its ability to handle nonlinearities of multi-variable systems as well as all the system’s constraints, i.e., input and state limits. The primary purpose of the linear and nonlinear MPC is to apply a plant model to predict its behavior along a receding horizon. At each sampling time, the MPC technique computes the optimal control input that minimizes some cost functions and fulfills all the safety constraints (e.g., actuators position, speed, and acceleration bounds). Overall, it has been proven that the MPC technique is much more practical and sophisticated to solve linear and nonlinear control problems [9,10,11,12,13]. Amongst the existing works, a real-time constrained MPC with safety and stability constraints is proposed in [9,10], where all the constraints are expressed as inequalities with respect to the optimization variable. Furthermore, a linear time-varying MPC method is designed in [12,13]. It is based on the online exact linearization of the nonlinear model so as to formulate the optimization problem completely as a Quadratic Program (QP). It can be easily handled by linear solvers to decrease the computational complexity. Nevertheless, they assume that uncertainties, related to the system’s modeling and the degree of correctness of measurements, are omitted. Consequently, a control strategy in a guaranteed way is required to ensure robustness toward all the uncertainties occurring in dynamic parameters’ estimation.

Currently, many researchers apply IA methods to synthesize a guaranteed linear or nonlinear MPC controllers [14,15]. These guaranteed methods allow the calculation of smooth and safe input intervals over the prediction horizon, i.e., the control vector is expressed here as an interval-vector (or a box that is the Cartesian product of intervals). Appendix A gives a brief presentation outlining different IA tools. For instance, authors in [14] use IA methods to build a new nonlinear MPC law considering the discrete system. It is based on forward-backward contraction using a dynamic model to compute feasible and admissible inputs. This method allows authors to calculate a validated control for an inverse pendulum system. Other techniques include the adjoint method for sensitivity analysis enabling the handling of the rate of change of the system dynamics considering the control variable [6,16,17,18]. Otherwise, this last approach is used to calculate the derivative of dynamics with respect to the input variable that must be limited between some bounds. In addition, a robust and guaranteed MPC of constrained discrete-time nonlinear systems with additive uncertainties is designed in [6,17] to quantify different uncertainties. Moreover, a branch-and-bound algorithm was applied in [15] to synthesize a NMPC controller. This combinatorial optimization algorithm is promising but needs high computation time to obtain optimal solutions. In addition, a reliable nonlinear MPC via interval arithmetic is developed in [19]. This algorithm uses validated simulation and branching methods to create input interval-vectors (or input boxes) and respect the system’s constraints by consecutive bisection of the initial feasible boxes of actuator’s inputs.

Contributions: This paper focuses on developing a new validated and reliable nonlinear model predictive control (NMPC), which is more specifically designed using an uncertain mathematical model. It encompasses two main research dimensions. (i) This paper is an extension of our previous work introduced in [20] that aims to identify in a guaranteed way the inertial and frictional parameters of our new inverted pendulum (manufactured in our lab and shown in Figure 1). The IA and set-inversion tools have been deployed to express all the needed parameters as intervals instead of constants’ values so that the coverage ratio between the physical reality and the mathematical model is significant. (ii) The second contribution lies in the synthesis of a validated NMPC strategy relying on the previous guaranteed identification to compute smooth and safe input intervals. These approaches have been investigated through various simulations with uncertain high-order ordinary differential equations (ODEs) describing the system’s behavior. The validated simulations are performed using our DynIbex (DynIbex Library: https://perso.ensta-paris.fr/chapoutot/dynibex/ (accessed on 19 August 2021)) library (a sophisticated solver of ODEs using Runge–Kutta schemes [21]) and several experiments on the actual nonlinear inverted pendulum.

The remainder of this paper is organized as follows: firstly, an overview of the direct dynamic model of our inverted pendulum is presented in Section 2. Secondly, Section 3 introduces the guaranteed NMPC approach based on IA. Thirdly, simulation and experimental results and discussions are reported in Section 4, on which the guaranteed identification results using the SIVIA algorithm (Set-Inversion Via Interval Analysis [7]) are recapped, and the capabilities of the proposed controller are illustrated via the nonlinear inverted pendulum. Finally, Section 5 closes the paper by summarizing the contribution and providing an overview of the future work.

2. Background of Dynamic Modeling of the Inverted Pendulum

2.1. Robotic Device

The inverted pendulum is a simple system with two serial rotational joints (active and passive). We have designed and manufactured one in our laboratory (see Figure 1). Our experimental device is equipped with a brushless DC motor to actuate the first joint and two incremental encoders to measure the rotational angles of both the motor shaft and the passive joint. Its kinematics is sketched in Figure 1. This device must be identified and stabilized in a validated and guaranteed way, considering all the errors and uncertainties related to the system modeling and data measurements.

2.2. Dynamic Modeling of the Inverted Pendulum

The Euler–Lagrange equation is a systematic method used here to express the system’s equations of motion, describing various nonlinear derivatives to get the ODEs (for more details, see [22]). After going through this process, the nonlinear equations of a direct dynamic model for rotary nonlinear inverted pendulum are given by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{y}}_1=y_2,\hfill \\ {\dot{y}}_2={\displaystyle \frac{1}{\Delta }}\{-{\mu }_{3}cos\left(y_3\right)\left[{\mu }_{1}sin\left(y_3\right)cos\left(y_3\right)y_2^2+{\Gamma }_2-{\mu }_{g}sin\left(y_3\right)\right]+\hfill \\ {\mu }_4\left[{\mu }_{3}sin\left(y_3\right)y_4^2-{\mu }_{1}sin\left(2y_3\right)y_2{y}_4+{\Gamma }_1\right]\},\hfill \\ {\dot{y}}_3=y_4,\hfill \\ \\ {\dot{y}}_4={\displaystyle \frac{1}{\Delta }}\{\left[{\mu }_{1}sin{\left(y_3\right)}^2+{\mu }_2\right]\left[{\mu }_{1}cos\left(y_3\right)sin\left(y_3\right)y_2^2+{\Gamma }_2-{\mu }_{g}sin\left(y_3\right)\right]-\hfill \\ {\mu }_{3}cos\left(y_3\right)\left[{\mu }_{3}sin\left(y_3\right)y_4^2-2{\mu }_{1}cos\left(y_3\right)sin\left(y_3\right)y_2{y}_4+{\Gamma }_1\right]\},\hfill \end{array}\right.\end{array}$$

(1)

where $y_1$ and $y_2$ indicate, respectively, the angular position and velocity of the rotative arm. Likewise, $y_3$ and $y_4$ are the angular and speed variables of the pendulum arm. We draw all these agreements as $y=[y_1,y_2,y_3,y_4]\in {\mathbb{R}}^4$ referring to the state vector. ${\Gamma }_i$ is the torque applied at the joint i, which depends on the viscous friction variables, and

$$\Delta =-{\mu }_3^{2}cos{\left(y_3\right)}^2+{\mu }_1{\mu }_{4}sin{\left(y_3\right)}^2+{\mu }_2{\mu }_4,$$

$${\mu }_1=l_p^2\left[{\displaystyle \frac{m_p}{4}}+M\right],{\mu }_2=l_a^2[m_p+M]+{\displaystyle \frac{J_m}{N_g^2}}+J_a,$$

$${\mu }_3=l_p{l}_a\left[{\displaystyle \frac{m_p}{2}}+M,\right],{\mu }_4=l_p^2\left[{\displaystyle \frac{m_p}{4}}+M\right]+J_p,{\mu }_g=\left[{\displaystyle \frac{m_p}{2}}+M\right]l_{p}g,$$

where $m_a$ and $J_a$ denote the horizontal arm’s mass and its inertia, respectively—similarly for the pendulum arm with the mass $m_p$ and inertia $J_p$. M is the mass of the load attached to the pendulum arm, $J_m$ is the DC motor inertia, and $N_g$ its gear ratio. $l_a$, $l_p$, and $r_0$ are the device’s dimensional parameters (as defined in Figure 1), and g is the gravitational acceleration.

Since the second joint is passive, ${\Gamma }_2$ can only be written as a function of friction conditions. The relationship between friction and velocity is well modeled in the literature [20]. In this paper, we use a simple friction model, in which it is related to the viscous terms, which are proportional to the joint speed:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\Gamma }_1\\ \\ {\Gamma }_2\end{array}\right]=\left[\begin{array}{c}\tau -f_{v1}{y}_2\\ \\ -f_{v2}{y}_4\end{array}\right],\end{array}$$

(2)

where $\tau$ is the motor’s torque and $f_{v_i}$ is the viscous Coulomb friction coefficient of the joint i. These coefficients are often considered as constant values in the majority of literature works. Here, the main goal is to identify them in a guaranteed way by reducing the error between the theoretical model output and measured data output, i.e., they must belong to specific subsets consistent with the model and sensors’ outputs. The guaranteed estimated parameters are reported in [20] and outlined in

Table 1. Identified inertial and frictional parameters of the inverted pendulum.

Symbol	Measurement Unit	Interval Enclosure
${\mu }_1$	kg·m${}^2$	$[9.63318e^{-4};1.22604e^{-3}]$
${\mu }_2$	kg·m${}^2$	$[2.19585e^{-3};2.79471e^{-3}]$
${\mu }_3$	kg·m${}^2$	$[7.227e^{-4};8.833e^{-4}]$
${\mu }_4$	kg·m${}^2$	$\left[1.08271e^{-3};1.37799e^{-3}\right]$
${\mu }_g$	kg·m${}^2$·s${}^{-2}$	$[6.24299e^{-2};8.44641e^{-2}]$
$f_{v_1}$	N·m·s·rad${}^{-1}$	$[0.043012;0.13002]$
$f_{v_2}$	N·m·s·rad${}^{-1}$	$[0.000454;0.001174]$

.

2.3. Torque, Speed, and Acceleration Computation

Because our experimental datasets have noises, the recorded signals are consequently processed by applying an eighth-order low-pass Butterworth filter. The input of the system (torque signal) and the joints’ velocity and acceleration are required to stabilize the nonlinear inverted pendulum. The speed and acceleration are computed using a non-causal derivative filter. The function firpm of Matlab allows us to synthesize this kind of filter to determine the suitable filter’s parameters (impulse response coefficients). The second stage is the actual application of the designed filter to the joints’ positions ($y_1$ and $y_3$), using filtfilt command, enabling to process a data in both the forward and reverse directions (zero-phase distortion).

On the other hand, the torque generated at the base of the rotary arm (i.e., at the load gear) can be expressed as a function of the angular velocity (input variable) denoted u as below:

$$\begin{array}{c}\hfill \tau ={\displaystyle \frac{N_m{N}_g{\eta }_m{\eta }_g(u_0-u)}{R_v}}\end{array}$$

(3)

where $N_m$ and $N_g$ are respectively the motor and gear ratios, ${\eta }_m$ and ${\eta }_g$ are the motor and gearbox efficiencies, u and $u_0$ are the instantaneous and no-load angular velocities, respectively, and $R_v$ is the motor speed-torque constant.

3. Guaranteed Controller Design via Validated Simulation

In this section, a validated Nonlinear Model Predictive Control (NMPC) is developed to stabilize the inverted pendulum in the upright position considering the system’s constraints and uncertainties. Firstly, we give a general overview of the validated numerical integration methods in Section 3.1. Secondly, we present the validated NMPC in Section 3.2.

3.1. Validated Numerical Integration

The dynamic model given by Equation (1) can simply be written by this nonlinear ODEs:

$$\begin{array}{c}{\dot{y}}_t=\mathbf{f}(t,y_t,u_t)\hfill \end{array}$$

(4)

At each validated numerical integration, we assume that $u_t$ is constant to compute corresponding tight enclosures of the state variable. We are now considering an Initial Value Problem for ODEs (IVP-ODEs) over the time interval $t\in [0,T]$,

$$\left\{\begin{array}{c}{\dot{y}}_t=\mathbf{f}(t,y_t,u_t)\hfill \\ y_0\in \left[{\mathbf{y}}_0\right]\subseteq {\mathbb{R}}^n\hfill \\ u_t\in \left[{\mathbf{u}}_{}\right]\subseteq {\mathbb{R}}^m\hfill \end{array}\right.$$

(5)

where the sets $\left[{\mathbf{y}}_0\right]$ and $\left[{\mathbf{u}}_{}\right]$, both expressed as boxes, are respectively the initial condition related to the initial state vector and the considering input interval. This IVP-ODE has a unique solution $y_t(t;y_0,u_t)$ at $t>0$ since $\mathbf{f}:\mathbb{R}\times {\mathbb{R}}^n\times {\mathbb{R}}^m\to {\mathbb{R}}^n$ is continuous in t and Lipschitz in $y_t$, but, for our purpose, we assume $\mathbf{f}$ to be smooth enough, i.e., of class $C^k$.

Compared to the classical numerical integration for the IVP-ODE problem, validated approaches consist of solving this problem in a complete and validated way (i.e., each actual guaranteed solution is rigorously returned and enclosed in a tight interval). These methods are generally based on Taylor series [23] or Runge–Kutta methods [21]. Basically, the main principle of the validated methods is to obtain the tight enclosure of the IVP-ODEs problem. As presented in [23], the purpose of a validated numerical algorithm is to solve Equation (5) so as to get a sequence of boxes $\left[{\mathbf{y}}_0\right],\dots ,\left[{\mathbf{y}}_n\right]$ at time iterations $t_0=0<\dots <t_n=T$ and at given input boxes $\left[{\mathbf{u}}_0\right],\dots ,\left[{\mathbf{u}}_n\right]$. It is achieved in such a way that the inclusion function verifies $\left[{\mathbf{y}}_{j+1}\right]\supseteq \left[\mathbf{f}\right]\left(t_j,\left[{\mathbf{y}}_j\right],\left[{\mathbf{u}}_j\right]\right),\forall j\in [0,n]$. Such approaches work in two stages at each integration step to compute the guaranteed solutions. They are summarized as:

i

Computation of a prior enclosure of the solution $\left[{\tilde{\mathbf{y}}}_{j+1}\right]$, such that $\mathbf{f}\left(t_j,\left[{\mathbf{y}}_j\right],\left[{\mathbf{u}}_j\right]\right)\in \left[{\tilde{\mathbf{y}}}_{j+1}\right]$ for all t in the time interval $t\in [t_j,t_{j+1}]$. This stage enables proving the existence and the uniqueness of the solution.

ii

Computation of a tight enclosure of state variable $\left[{\mathbf{y}}_{j+1}\right]$ at time instance $t_{j+1}$, such that $\mathbf{f}\left(t_{j+1},\left[{\mathbf{y}}_j\right],\left[{\mathbf{u}}_j\right]\right)\in \left[{\mathbf{y}}_{j+1}\right]$. It uses the solution $\left[{\tilde{\mathbf{y}}}_{j+1}\right]$ to bound the truncation error, i.e., the distance between the exact and the numerical solutions.

The tight and prior enclosures calculated along one integration step between $t_j$ and $t_{j+1}$ (of size $h_j=t_{j+1}-t_j>0$) can be visualized in Figure 2.

3.2. Validated Nonlinear Model Predictive Control

The purpose of this section is to apply the validated simulation algorithm presented in the previous section to compute a solution satisfying the NMPC constraints and minimizing an interval cost function.

The validated NMPC aims to stabilize the nonlinear inverted pendulum in the upright position. It is based on interval analysis tools to compute guaranteed control variables over the receding horizon using the validated numerical integration method of the IVP-ODEs problem (5). The control intervals are calculated such that the convergence to the set-point interval is ensured (i.e., $y_i\to \left[{\mathbf{y}}_r\right],\forall i$), and all the state and input constraints are satisfied (i.e., $y_i\in \left[{\mathbf{y}}_i\right]$ and $u_j\in \left[{\mathbf{u}}_j\right],\forall i,j$). The validated NMPC is decomposed into two steps:

• The IVP-ODEs problem (5) is simulated in a complete and validated way starting from the initial domain of inputs referring to the actuators’ saturation. Then, for each tested input set, we check if the constraints on state variables are respected using validated simulation and branching methods (Branch-and-Prune algorithm).
• Since safe inputs are computed over a finite time horizon, the optimal controller box is calculated by minimizing the interval objective function, optimizing the error between the predicted outputs and the reference set, and the norm of the inputs.

3.2.1. Filtering and Branching Procedures

The validated NMPC strategy should provide an optimal input box at each time-step k. This latter can be denoted as ${\left[\mathbf{U}\right]}_k^{☆}={\left[\mathbf{u}\right]}_k^{☆}\times {\left[\mathbf{u}\right]}_{k+1}^{☆}\times \dots \times {\left[\mathbf{u}\right]}_{k+N_p-1}^{☆}$, where $N_p$ is the number of input samples generated over the prediction horizon. Branching procedure allows us to generate safe input intervals ${\left[\mathbf{U}\right]}_k^{☆}$ over the prediction horizon that satisfies all the state constraints and convergence to the reference interval $\left[{\mathbf{y}}_r\right]$. Therefore, the system can be described using a simple interval notation by the following constraint differential equation:

$$\begin{array}{c}{\dot{y}}_t=\mathbf{f}(t,y_t,u_t),\hfill \\ \mathrm{s}.\mathrm{t}.:u_t\in \mathbb{U},\forall t\ge 0\hfill \\ y_t\in \mathbb{Y},\forall t\ge 0\hfill \end{array}$$

(6)

where $\mathbb{U}$ and $\mathbb{Y}$ are respectively the input and state variables saturation. They can be written as

$$\mathbb{U}=\left\{u_t\in {\mathbb{R}}^m|u_{min}\le u_t\le u_{max}\right\}\mathrm{and}\mathbb{Y}=\left\{y_t\in {\mathbb{R}}^n|y_{min}\le y_t\le y_{max}\right\}$$

(7)

where $u_{min}$ and $u_{max}$ are the admissible domain bounds for the control variable. Similarly, $y_{min}$ and $y_{max}$ are those of the state components.

Starting from the current state boxes $\left[{\mathbf{y}}_t\right]$ computed at time t and input feasible domain $\left[{\mathbf{u}}_k\right]=[u_{min},u_{max}]$, the computation of a new state domain $\left[{\mathbf{y}}_{t+T_c}\right]$ is performed by the validated simulation algorithm, where $T_c$ denotes the sampling time of the control signal, and then the horizon time $T_p=T_c\times N_p$. This strategy bisects the initial domain $\left[{\mathbf{u}}_k\right]$ on the left and right intervals that reduces the width of the initial interval $\left[{\mathbf{u}}_k\right]$. For each bisected input interval, the validated simulation algorithm provides a state domain of which the safety and convergence criteria are verified, as can be seen in this equation:

$$\left[\mathbf{f}\right](t,\left[{\mathbf{u}}_k\right],\left[{\mathbf{y}}_t\right])=\left[{\mathbf{y}}_{t+T_c}\right]\mathrm{and}\left[{\mathbf{y}}_{t+T_c}\right]\subseteq [y_{min},y_{max}]$$

(8)

The filtering and branching algorithm is given in Algorithm 1. The function bisect-check is a pruning function that performs a bisection procedure of the initial input interval if its width exceeds some tolerance $tol$ value. One side of bisected intervals is kept by considering these criteria: (i) A branch leading to unsafe state is avoided (i.e., $\left[{\mathbf{y}}_{t+T_c}\right]⊈[{\mathbf{y}}_{min},{\mathbf{y}}_{max}]$); and (ii) a branch leading to a state far from the reference interval $\left[{\mathbf{y}}_r\right]$ or partially a cost greater than the other branch is eliminated (see interval optimization in Section 3.2.2). Moreover, the function Test sensitivity() enables us to simulate the adjoint function $x\left(u\right)={\displaystyle \frac{\partial y}{\partial u}}$ expressing the rate of variation of the state concerning the control variable. Its derivative can be derived as below:

$$\dot{x}={\displaystyle \frac{\partial \mathbf{f}}{\partial y}}\times x+{\displaystyle \frac{\partial \mathbf{f}}{\partial u}}$$

(9)

Using the Test sensitivity() function, it returns the Boolean formula:

$$\left\{\exists ϵ>0,0\in x\left(\left[{\mathbf{u}}_{\mathbf{k}}\right]\right)\mathrm{AND}x\left(\left[{\mathbf{u}}_k\right]\right)⊈[-ϵ,ϵ]\right\}\mathrm{OR}\left\{\exists \tau \in [t,t+T_c],y\left(\tau \right)⊈\mathbb{Y}\right\}$$

(10)

Sensitivity analysis allows for checking if the computed input $\left[{\mathbf{u}}_{\mathbf{k}}\right]$ is leading the system to move in forward and backward directions at the same time, as can be seen by the first part of the Boolean expression. The second part of this formula enables us to verify if the state constraints are respected. Since the safe inputs’ intervals are computed using the validated simulation, the interval cost function is optimized to calculate its sub-optimal solution. The optimization part is detailed in the next section.

Algorithm 1 Filtering and branching

Require:$\left[{\mathbf{y}}_0\right],\left[{\mathbf{y}}_r\right],T_c,N_p,tol$

while not success do

while $k\le N_p$ do

if $w\left(\left[{\mathbf{u}}_k\right]\right)>tol$ then

$\left[{\mathbf{u}}_k\right]=\mathrm{bisect}-\mathrm{check}({\left[{\mathbf{u}}_k\right]}_{\mathrm{left}},{\left[{\mathbf{u}}_k\right]}_{\mathrm{right}})$

end if

$\left[{\mathbf{y}}_{t+T_c}\right]=$ Simulation $(\left[{\mathbf{y}}_t\right],\left[{\mathbf{u}}_k\right],T_c)$

if $\left[{\mathbf{y}}_{t+T_c}\right]\subseteq \left[{\mathbf{y}}_r\right]$ then

$k=k+1,\left[{\mathbf{y}}_t\right]=\left[{\mathbf{y}}_{t+T_c}\right],t=t+T_c$

$success=true$

else

if Test sensitivity() then

if $w\left(\left[{\mathbf{u}}_k\right]\right)>tol$ then

$\left[{\mathbf{u}}_k\right]=\mathrm{bisect}-\mathrm{check}({\left[{\mathbf{u}}_k\right]}_{\mathrm{left}},{\left[{\mathbf{u}}_k\right]}_{\mathrm{right}})$

else

$success=false$

end if

else

$\left[{\mathbf{y}}_t\right]=\left[{\mathbf{y}}_{t+T_c}\right]$

$k=k+1$, $t=t+T_c$

$success=true$

end if

end if

end while

end while

3.2.2. Interval Cost Function Formulation

Since the feasible intervals of inputs are computed fulfilling different constraints, the optimization procedure is required to find the optimal control box and then the punctual input that we inject into the actuator’s control unit $u_1$. The continuous cost function can be expressed over the prediction horizon $T_p=T_c\times N_p$ as

$$J(y_t,u_t)={\int }_t^{t+T_p}F(y\left(\tau \right),u\left(\tau \right))d\tau$$

(11)

where F is in general a quadratic function that minimizes the norm of the inputs and the error between the predicted outputs y and the reference $y_r$. F can be expressed in the following quadratic form:

$$F(y\left(\tau \right),u\left(\tau \right))={\left(y\left(\tau \right)-y_r\right)}^T\mathbf{Q}\left(y-y{\left(\tau \right)}_r\right)+u{\left(\tau \right)}^T\mathbf{R}u\left(\tau \right)$$

(12)

where $\mathbf{Q}$ and $\mathbf{R}$ are both positive semi-definite weighting matrices.

The continuous objective function (11) then takes the form

$$\begin{array}{cc}\hfill J(y_t,u_t)& ={\int }_t^{t+T_p}\left[{\left(y\left(\tau \right)-y_r\right)}^T\mathbf{Q}\left(y\left(\tau \right)-y_r\right)+u{\left(\tau \right)}^T\mathbf{R}u\left(\tau \right)\right]d\tau \hfill \\ \hfill & =\sum _{k=1}^{N_p}{\int }_{t+(k-1)T_c}^{t+k{T}_c}\left[{\left(y\left(\tau \right)-y_r\right)}^T\mathbf{Q}\left(y\left(\tau \right)-y_r\right)\right]d\tau +\sum _{k=1}^{N_p}{\int }_{t+(k-1)T_c}^{t+k{T}_c}u{\left(\tau \right)}^T\mathbf{R}u\left(\tau \right)d\tau \hfill \end{array}$$

(13)

The purpose is to express the cost function (13) in a discrete form using interval analysis methods. Since the control signal is constant along the input sampling period $[t+(k-1)T_c,t+k{T}_c]$, the second integral can be simplified as

$$\begin{array}{ccc}\hfill J(y_t,u_t)& \hfill & \hfill =\sum _{k=1}^{N_p}{\int }_{t+(k-1)T_c}^{t+k{T}_c}\left[{\left(y\left(\tau \right)-y_r\right)}^T\mathbf{Q}\left(y\left(\tau \right)-y_r\right)\right]d\tau +T_c\sum _{k=1}^{N_p}\left[u_k^T\mathbf{R}{u}_k\right]\end{array}$$

(14)

A validated simulation method generates a list of tight enclosures’ $\left[{\mathbf{y}}_0\right],\left[{\mathbf{y}}_1\right],\dots ,\left[{\mathbf{y}}_n\right]$ solutions of the IVP-ODEs (5) from the initial conditions $\left[{\mathbf{y}}_0\right]$ and $\left[{\mathbf{u}}_0\right]$. Therefore, each solution y can be enclosed into a box thanks to the validated methods (i.e., $\forall t\in [t+(k-1)T_c,t+k{T}_c],y\left(t\right)\in \left[{\mathbf{y}}_t\right]$). Hence, one useful consequence of the Rectangle Property is that it allows us to bound the first integral of the cost (14) between two bounds as can be defined below:

$$\begin{array}{c}\hfill \sum _{k=1}^{N_p}{\int }_{t+(k-1)T_c}^{t+k{T}_c}\left[{\left(y\left(\tau \right)-y_r\right)}^T\mathbf{Q}\left(y\left(\tau \right)-y_r\right)\right]d\tau \in T_c\sum _{k=1}^{N_p}\left[{\left(\left[{\mathbf{y}}_t\right]-y_r\right)}^T\mathbf{Q}\left(\left[{\mathbf{y}}_t\right]-y_r\right)\right]\end{array}$$

(15)

By exploiting (15), the interval cost function can be approximated as

$$\begin{array}{cc}\hfill J(y_t,u_t)& \in T_c\sum _{k=1}^{N_p}\left[{\left(\left[{\mathbf{y}}_t\right]-y_r\right)}^T\mathbf{Q}\left(\left[{\mathbf{y}}_t\right]-y_r\right)+{\left[{\mathbf{u}}_k\right]}^T\mathbf{R}\left[{\mathbf{u}}_k\right]\right]\hfill \\ \hfill & \le ub\left\{T_c\sum _{k=1}^{N_p}\left[{\left(\left[{\mathbf{y}}_t\right]-y_r\right)}^T\mathbf{Q}\left(\left[{\mathbf{y}}_t\right]-y_r\right)+{\left[{\mathbf{u}}_k\right]}^T\mathbf{R}\left[{\mathbf{u}}_k\right]\right]\right\}\hfill \end{array}$$

(16)

where $ub$ signifies the upper bound of the interval ($lb$ denotes its lower bound).

Algorithm 2 recaps the optimization procedure to determine the inputs boxes allowing for minimizing the cost function (16). It returns the sub-optimal solutions that minimize the formulated cost function. It is based on consecutive bisections of the first input interval $\left[{\mathbf{u}}_{\mathbf{1}}\right]$ until achieving a small width value $tol$. The punctual and guaranteed input value $u_1$ that we send to the DC motor represents the smaller value of the upper and lower bound of this interval.

Algorithm 2 Optimization

Require:$\left[\mathbf{y}\right],\left[{\mathbf{u}}_1\right],\left[{\mathbf{u}}_2\right],\dots ,\left[{\mathbf{u}}_{N_p}\right],tol$

while$w\left(\left[{\mathbf{u}}_1\right]\right)\ge tol$do

${\left[\mathbf{U}\right]}_{\mathrm{left}}={\left[{\mathbf{u}}_1\right]}_{\mathrm{left}}\times \left[{\mathbf{u}}_2\right]\times \dots \times \left[{\mathbf{u}}_{N_p}\right]$

${\left[\mathbf{U}\right]}_{\mathrm{right}}={\left[{\mathbf{u}}_1\right]}_{\mathrm{right}}\times \left[{\mathbf{u}}_2\right]\times \dots \times \left[{\mathbf{u}}_{N_p}\right]$

if $J(\left[\mathbf{y}\right],{\left[\mathbf{U}\right]}_{\mathrm{left}})\ge J(\left[\mathbf{y}\right],{\left[\mathbf{U}\right]}_{\mathrm{right}})$ then

${\left[\mathbf{U}\right]}^{☆}={\left[\mathbf{U}\right]}_{\mathrm{right}}$

else

${\left[\mathbf{U}\right]}^{☆}={\left[\mathbf{U}\right]}_{\mathrm{left}}$

end if

end while

$u_1=min\left\{lb\left(\left[{\mathbf{u}}_1\right]\right),ub\left(\left[{\mathbf{u}}_1\right]\right)\right\}$

$send\left(u_1\right)$

To summarize, the proposed approach is focused on two main stages:

Algorithm 1:

• returns the safe input tubes ($N_p$ intervals) that allow for respecting the state constraints. This algorithm uses the validated simulation methods to compute the guaranteed inputs. As already mentioned, we assume that the input interval $\left[{\mathbf{u}}_t\right]$ is constant for each validated simulation to calculate the state tubes equivalent to each input interval. If the state limits are not fulfilled, the bisection procedure of the control interval is performed, and validated simulations are re-started.

Algorithm 2:

• Because of the guaranteed input tubes ($N_p$ intervals), the optimization algorithm is started to compute the sub-optimal input interval that minimizes as much as we can the formulated interval cost-function. It is based on the consecutive bisections of the guaranteed interval $\left[{\mathbf{u}}_1\right]$ computed by the first algorithm.

4. Experiment Results

The control framework introduced in this paper is applied for the stabilization of the nonlinear inverted pendulum shown in Figure 1. The pendulum is actuated by a DC motor whose angular speed is the input variable. It is controlled using a Phidget Motor Control. This control card is connected to the laptop via a USB cable to control the motor’s velocity, using C language.

To evaluate the validated NMPC, we will initially summarize in Section 4.1 the guaranteed identification results and its validation using the identified subsets of parameters. In the second place, we will validate the proposed controller in Section 4.2 through simulation tests and then experiments using the inverted pendulum.

4.1. Guaranteed Dynamic Model Identification

The guaranteed identification of dynamic parameters $p=[{\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4,{\mu }_g,f_{v_1},f_{v_2}]\in {\mathbb{R}}^{n_p=7}$ as can be seen by the dynamic model (1) was proposed in [20]. These parameters are estimated as intervals instead of real numbers to account for all the uncertainties and errors related to the sensors measurements and the system modeling. The identification approach is based on the IA and set-inversion techniques (SIVIA algorithm) to determine the feasible sets of coefficients [7]. A brief resume of the SIVIA algorithm principle is given in Appendix A. The feasible sets of all the dynamic parameters are recapped in Table 1.

The DC motor and gearbox are of type MDP-Maxon, and their characteristics are provided by the manufacturer. Here is a recap of their properties: the motor inertia $J_m$ = $9.06e^{-7}$ kg·m${}^2$, its no-load speed $u_0$ = 11,700 tr·min${}^{-1}$, the slope speed-torque $Rv$ = 34,000 tr·min${}^{-1}$·N${}^{-1}$·m${}^{-1}$, the motor and gearbox efficiencies ${\eta }_m=90.6\%$ and ${\eta }_g=90\%$, and gear ratios $N_m=5.2$ and $N_g=5$.

To demonstrate the usefulness of our identification, we can solve in a complete and validated way the IVP-ODEs (5) using the identified parameters. The set of possible tight enclosures $\left[{\mathbf{y}}_1\right],\dots ,\left[{\mathbf{y}}_n\right]$ at a sequence of time-instants $t_1,\dots ,t_n$ are calculated. To do so, our library DynIbex is employed to solve this IVP-ODEs [21]. Furthermore, the guaranteed identification is compared to the traditional least square method identification (LSMI) (where the dynamic parameters are estimated as constant values) to show its efficiency. Figure 3 shows the measured pendulum angle $y_3$ (black lines), as well as the ones calculated with the identified parameters with LSMI (red lines) and IA (blue lines) approaches. All the model validation tests are conducted under the same settings and conditions (e.g., same initial conditions $\left[{\mathbf{y}}_0\right]$ and $\left[{\mathbf{u}}_0\right]$, same inputs, etc.). The simulation duration is fixed to be 4 s and the precision of ${10}^{-7}$. We can notice that the simulated tubes of the pendulum angle $y_2$ given by the sophisticated solver DynIbex, using IA and LSMI parameters, are quite similar to the real measured pendulum angle. Nonetheless, the simulated $y_2$ with IA parameters comply as closely as possible with the measured one at different initial conditions. All this is due to the fact that all the system variables are taken here as intervals characterizing different uncertainties.

In order to emphasize the potential capabilities of the predicted model, the associated coverage rates can be computed. We evaluate the number of cases where the measurements data are included in the simulated tubes, and then the percentage of inclusion is calculated for each tested scenario. In other words, it is calculated by checking at each time step whether the time variable related to the measurement is included in the simulated time interval $t_i\in \left[{\mathbf{t}}_{\mathbf{i}}\right]$. If so, we inspect if the corresponding measure $y_i$ is included in the simulated tube $\left[{\mathbf{y}}_i\right]$ given by the DynIbex library. Then, the coverage rate can be approximated to be $\frac{N_q}{N_t}\times 100$, where $N_t$ is the number of times that each sampled time is incorporated in the simulated time interval, and $N_q$ is the one related to the second condition regarding $y_i$. The computed values are recapped in

Table 2. Degree of coverage between the model and physical reality.

Scenario	(a)	(b)	(c)	(d)
With IA parameters	46%	61%	38%	25%
With LSMI parameters	41%	34%	20%	21%

. By the way, this criterion is too strict to test, which explains the small calculated percentages. This is due to the timing lag between the measurement and the estimated model. Even with this issue, the coverage ratios obtained using IA parameters indicate a good match with the measurements process. It also confirms that the model is estimated with high precision when the uncertainties are considered.

Overall, the experimental results provided by the guaranteed identification of inertial and frictional parameters are more relevant than those obtained by the classical LSMI. Indeed, the estimation approach based on IA and set inversion techniques allows us to successfully identify these coefficients as intervals considering all the errors and uncertainties related to the sensors and the system modeling. Despite its proven effectiveness to more closely approximate the model output to the actual system’s behavior, it still has some drawbacks. The main ones are related globally to the computation time, which depends mainly on the number of parameters and the initial domain. This issue can be enhanced with a contraction approach [24].

4.2. NMPC Results

4.2.1. Simulation Results

In the simulations, the parameters of this proposed validated NPMC are tuned as: prediction horizon $N_p=10$, control sampling time $T_c=16$ ms and the final time of the simulation $T_f=0.3$ s. The safety saturation on state variables are: $y_1\in [-2\pi ;2\pi ]$, $y_2\in [-52.4$ rad·s${}^{-1};52.4$ rad·s${}^{-1}]$, $y_3\in [-2\pi ;2\pi ]$, $y_4\in [-100$ rad·s${}^{-1};100$ rad·s${}^{-1}]$, and the DC motor’s torque constraints are $\tau \in [-8.05$ N·m; 8.05 N·m]. The reference interval of the desired pendulum arm is $y_r\in [\pi -0.1;\pi +0.1]$. The tolerance parameter applied for the bisection procedure is adjusted as $tol=0.25$. Since the filtering and branching algorithm begins from the admissible input domain $[-8.05$ N·m; 8.05 N·m], it can lead to approximately ${64}^{10}$ branches with $N_p=10$. The interval cost function is computed only in the optimization process using these weighting matrices $\mathbf{Q}=diag\left[1000,1000,1000,1000\right]$ and $\mathbf{R}=1$.

Several software tools have been developed to simplify the resolution of the validated NMPC strategy via interval analysis. These tools are globally in the form of libraries to be integrated into a program (e.g., DynIbex library [21], Numerica [25] or RealPaver [26]). However, it has been shown that these solvers are not much more practical due to their high computation time. In other words, the fast convergence of this kind of solvers is not guaranteed over a finite time (i.e., extreme long response time to get the solution of the problem). That’s why the real-time validation of this guaranteed approach is not easy to realize. To do so, it is divided into two steps: a numerically offline calculation of the input variables via validated simulations and then the injection of these inputs into the physical pendulum.

Two scenarios have been executed with different initial conditions on the angular pendulum position $\left[{\mathbf{y}}_{3_{ini}}\right]$ to discuss the capabilities of the proposed controller. Figure 4a–c display the results of the first scenario starting from $\left[{\mathbf{y}}_{3_{ini}}\right]=[0^{\circ };0^{\circ }]$ and Figure 5a–c show those of the second one starting from $\left[{\mathbf{y}}_{3_{ini}}\right]=[{149}^{\circ };{150}^{\circ }]$. As can be seen from Figure 4a, the pendulum arm starts from the downward position, and it is stabilized via the validated NMPC in its vertical position interval $\left[{\mathbf{y}}_r\right]$ with a small settling-time (around $t_{r_{5\%}}\approx 0.18$ s). Similarly, the simulation has also been executed when the initial position is close to the desired position $\pi$. The generated enclosures of the inverted pendulum position are plotted in Figure 5a. As the previous results, the validated NMPC succeeds in stabilizing the pendulum feasibly in its terminal position.

On the other hand, the calculated tubes of the pendulum’s velocity are depicted in Figure 4b and Figure 5b. As can be noticed in both of these scenarios, the validated velocity respects the imposed constraints and stagnate around zero when the pendulum reaches its final desired position. In addition, the torques applied by the DC motor for each tested scenario are shown in Figure 4c and Figure 5c. The relationship between the motor’s velocity and torque is given by Equation (3). In the first case, a strong input variable with the bang-bang process is applied to enable the pendulum to catch the equilibrium and stabilize in its terminal reference interval (see Figure 4c). Since the pendulum is stabilized, the torque applied by the system’s actuator changes the sign to maintain the reference interval. In the second case, the actuator applies maximum inputs for a short time to try to reach the desired interval (see Figure 5c). All the calculated torques fulfill the system’s intrinsic constraints.

4.2.2. Experimental Validation Using the Inverted Pendulum

To further show the usefulness of the proposed controller, it is tested on the experimental platform (Figure 1). As mentioned before, the real-time implementation of the validated strategy cannot be accomplished due to the high computation time of interval methods. In addition, we are currently working to improve this point by adopting parallel computation and relaxation methods. Consequently, we inject the computed input components into the control part of the pendulum that must somehow behave as the simulated model. Indeed, the dynamic model is proven to be as close as possible to the system’s actual behavior (see Section 4.1).

Figure 6 shows the measured pendulum angle for the tested scenarios, starting from $0^{\circ }$ and ${150}^{\circ }$. The stabilization is performed in an open-loop by incorporating the simulated input data into the DC motor controller (the computed inputs are drawn in Figure 4c and Figure 5c). Compared to the simulation results, the calculated control inputs stabilize the inverted pendulum in the steady-state regime. The second scenario gives coherent results enabling the pendulum to be stabilized around the target position $\pi$ thanks to the validated NMPC (see Figure 6 (Right)). However, for the first scenario (starting from $0^{\circ }$), we have some synchronization problems between the DC motor sampling period (around 0.028 s) and the sampling horizon of the NMPC method $T_c$. This is due to the lack of good deeds that the system’s actuator executes at the right moment and its sensibility to weak signals that destabilize it. Despite these issues, the proposed NMPC enables the system to reach another equilibrium position around $\pi /2$ (see Figure 6 (Left)), so we will try to handle this issue separately.

To sum up, the validated NMPC seems more promising for enabling the system to successfully converge to the reference target with high accuracy and small settling time. The imposed state and control bounds are respected, critical, and necessary for the pendulum’s stability and safety. Nevertheless, the computation time of the domain $\left[{\mathbf{u}}_k\right]$ is remarkable due to a large number of bisections of the initial input domain and the high computation time required by validated simulations.

5. Conclusions and Future Works

The research presented in this paper focuses on developing a reliable and validated nonlinear model predictive control (NMPC), which is more specifically designed using an uncertain mathematical model. The proposed controller is based on validated simulation mixed with interval analysis techniques to handle the system’s nonlinearities and uncertainties. Two different procedures are employed: branching procedure to calculate the valid input intervals from an uncertain nonlinear model, allowing the convergence to the set-point interval and the system constraints’ insurance. From the generated branches, the optimization procedure is applied to determine the sub-optimal control variable. The efficiency and robustness of the proposed method are validated through numerical and experimental tests using a new nonlinear inverted pendulum.

As regards ongoing works, we will focus on studying two issues: firstly, the computational complexity improvement of the validated simulation methods by taking into consideration contraction methods and parallel computing strategies; secondly, the projection of the validated NMPC on a complex mechatronic system (ROV—Remotely Operated underwater Vehicle). It allows us to perform some tasks (trajectory tracking) in the aquatic environment while respecting all the system’s constraints. In addition, another relevant aspect enabling completing a system’s autonomy is online estimation of the model’s parameters via nonlinear observers and interval methods with respect to the environmental conditions, e.g., aerodynamic coefficients. This issue allows for reaching high accuracy and stability of the trajectory tracking task.