The purpose of this first interval-based INS is to estimate reliably the heading ψ of the boat as an interval ψ , using the dual GPS of the boat, which gives an estimation ψ _GPS, the gyrometers of the AHRS ω (rotation around the vertical axis), and the magnetometers of the AHRS, which gives an estimation ψ _mag (we will assume here that the magnetic declination and any physical bias between the sensors have been already taken into account, or their uncertainties are within the width of the intervals). It is possible to describe our problem as follows: x ̇ = f x e v o l u t i o n    e q u a t i o n y = g x o b s e r v a t i o n    e q u a t i o n 0 = h x c o n s i s t e n c y    e q u a t i o n (1) where x ( t ) = cos ψ ( t ) sin ψ ( t ) is the state vector, y ( t ) = cos ψ GPS ( t ) sin ψ GPS ( t ) cos ψ mag ( t ) sin ψ mag ( t ) is the output vector, f ( x ) = 0 − ω ω 0 x , g ( x ) = x 1 x 2 x 1 x 2 , and h ( x ) = x 1 2 + x 2 2 − 1 . We will assume that:

• For all t ∈ t 0 , t f , we have intervals enclosing ω(t):

In 2D, the vertical gyrometer is measuring the heading velocity, so we have: ψ ̇ = ω . (5)

If we apply at each time step the contraction procedure described in the Data fusion with interval analysis section, the observation equation in combination with Eq. (3) will provide precise intervals for ψ at regular time instants ⁵ (from ψ GPS if the GPS is available and the distance between the two antennas is considered consistent with the user-provided antenna distance, or from ψ mag otherwise), while Eqs. (5) and (2) should propagate an estimation for other time intervals. For that, we could discretize the differential Eq. 5 using the Euler method, i.e., ψ t + d t = ψ t + d t ⋅ ω t ; however, in practice manipulations such as additions and comparisons directly on angles are not recommended due to the modulo 2π (or 360°) equivalence of angles. Instead, we can use the cos and sin of the angles, so if x ( t ) = cos ψ ( t ) sin ψ ( t ) , due to Eq. 5, we have x ̇ ( t ) = − ω ( t ) sin ψ ( t ) ω ( t ) cos ψ ( t ) (evolution equation). This can be written as x ̇ = A x with A = 0 − ω ω 0 , and additionally, it is possible to find an exact solution for this type of differential equation, as x t = k ⁡ exp A ⋅ t if ω is assumed to be a piecewise constant (i.e., A is not varying much during the time step dt). By writing x t + d t with respect to x t , it can be discretized as x t + d t = exp A d t x t . If A is antisymmetric, exp A d t is a rotation matrix, of angle ω ⋅ dt in our case. Although not strictly necessary, it is also possible to use the equation x 1 2 + x 2 2 = 1 (consistency equation) to try to limit potential overestimation of the uncertainty when evaluating cos ψ ( t ) and sin ψ ( t ) . Then the polar contractor from Desrochers and Jaulin (2016b) ⁶ can be used to contract ψ from cos ψ ( t ) and sin ψ ( t ) . If an empty set appears at any time during the contraction procedure at a time step, ψ is set to its last known value.

To get a first validation of the interval-based INS designed, comparisons with other systems have been made ⁷ . A first experiment has been set up with a high-end iXBlue QUADRANS Fiber-Optic Gyrocompass (FOG) INS (see Figure 2) as a kind of ground truth to get a first idea on the precision of our system. Both systems have been mounted together and moved manually with various rotations and oscillations at different speeds (only planar movements, see Figures 3 and 4). The mean error and the standard deviation between the systems are both around 1.5°, which means that there was probably a physical bias of 1.5° between the two systems, and assuming the QUADRANS is the ground truth (specifications suggest an error of less than 0.2°), the interval-based INS has a typical error of 1.5° ⁸ , which seems small enough for the autonomous boat used as a target system to be able to follow a track within the typical precision of a standard GPS, as it will be demonstrated in the Online control results in real experiments section (also, note that the typical error of an SBG Ellipse2-A-G4A3-B1, whose hardware—especially the gyrometers—is used in the interval-based INS, is around 1°).

The purpose of the second batch of experiments is to demonstrate that the interval-based INS could be competitive compared with existing solutions in some scenarios where faults or inconsistencies between sensors occur. Contrary to the preceding experiment, we are not interested in precision but rather in the behavior when handling obvious outliers. For that, two other systems were tested: an SBG Ellipse3-D-G4A2-B2 (with gyrometers full scale of 450°/s and accelerometers of 8 g) with SDK 6.2 and a Pixhawk 4 mini [firmware ArduRover V4.1.0-dev (12 645 674)] with simpleRTK2B+heading–Basic Starter Kit [firmware UBX_F9_100_HPG_113_ZED_F9P configured with instructions from ArduSimple (2021)]. Both systems are made of three magnetometers, three gyrometers, three accelerometers, and a dual GPS, and they are supposed to provide heading information from a fusion of all these sensors, similar to the interval-based INS of the boat. For all systems, default settings were used with magnetometers disabled (as long as the GPS is available) and dual antenna settings (approximately 1 m between the antennas, to fit the size of the boat). The antenna type is the same (antenna 1 at 46 cm behind the Ellipse-D, antenna 2 at 58.5 cm in front, set to “Rough lever arm” setting in sbgCenter), and initial date and lever arms are configured where applicable.

The test procedure has been the following:

• The INS is first left at a constant angle oriented toward the South (180°) during some time to ensure it is correctly initialized.

Figure 5 shows the results of the experiment for the Ellipse-D. The fused heading and the GPS heading coincide (except that the update rate of the GPS is 5 Hz, while the fused heading is determined at a frequency of 50 Hz) until the gyrometers are saturated. At this point, the fused heading outputs a wrong value of 200°, possibly evaluated from the gyrometer data; however, it does not correspond with the heading of the IMU, which should have been 270° or the heading of the GPS, which should have been 180° at this point of the experiment. After the IMU is put back in the same direction of the GPS, the fused heading indicates 100°, instead of 180°. It should be noted, however, that in other tests, where the gyrometer saturation was probably shorter, the Ellipse-D was sometimes able to recover from that disturbance. Also, although not tested, it might be possible to make it trust more the GPS instead of the gyrometers with other settings. Indeed, the commercial purpose of this type of INS is probably to be able to be more robust to GPS outliers due to, e.g., multipaths as well as increasing the heading output rate compared with a dual GPS without IMU aid. However, it should be possible to detect easily the gyrometer saturations and to choose in that case automatically to trust more the GPS data as long as they seem to be correct (e.g., distance between antennas still consistent with what was evaluated previously and with what was specified by the user as the initial value), as it was the case in this experiment.

Figure 6 shows the results of the experiment for the Pixhawk. The fused heading and the GPS heading coincide (except the update rate of the GPS is 5 Hz, while the fused heading is determined at a frequency of 50 Hz) until the gyrometers are saturated. At this point, the fused heading seems off but converges slowly toward the GPS value, until the IMU is moved back in the same direction of the GPS, when the gyrometers detect the rotation even though the fused heading seems again to converge finally toward the GPS value. This is probably due to a continuous weighted fusion between GPS and gyrometers data, which seems to be consistent with the fact that an EKF (Extended Kalman Filter) is known to be used internally (as this is an open-source design, this could be checked in more detail) ⁹ .

Figure 7 shows the results of the experiment for the interval-based INS. The fused heading and the GPS heading coincide (except that the update rate of the GPS is 1 Hz, while the fused heading is 20 Hz) except for a short time (1–2 s) when the gyrometers are saturated or when the IMU is moved back to its original position after the saturation. It could be interesting to check in further experiments whether a GPS rate of 5 Hz (like in the Ellipse-D) would give better results especially in the experiment with the QUADRANS, since it is not expected that it would significantly change the behavior in the experiment of Figure 7.

Although many improvements can already be foreseen, those results show that designing an interval-based INS naturally leads to consider which solution should be returned in case of inconsistencies in the data fusion process (i.e., empty intervals at some point in the computations). This can help to determine better when there are sensors that contradict each other and, whenever possible, which. Even if at the moment, it is not demonstrated that the designed interval-based INS would surpass other existing devices or methods, especially in terms of precision and accuracy, the foreseen theoretically increased reliability (with respect to inconsistencies, with a better propagation of all sensor and model errors, etc.) is alone worth to be reported, and this fact is already experimentally proven by the previously described experiments, which enable us to conclude that it appears to be better than COTS devices at least in some specific situations where inconsistencies could occur ¹⁰ .

Since GPS only provides a new position estimate every second, the boat model used in the Jet-boat simulator section could be used to enhance the state estimation with an improved temporal resolution, e.g., by following the method in Le Bars et al. (2018) (and the problem is somehow similar to the fusion between gyrometer and GPS heading described in the Principle of operation section). In this context, we can describe our problem as: x ̇ = f x , u evolution equation y = g x observation equation , (7) where x(t) is the state vector of the boat, u(t) is its input vector, y(t) is its output vector, and:

• For some time instant t ₀, we have a box x t 0 containing the state vector:

The evolution equation is given by Eq. (6), and the observation equation is: y = x y ψ (11) because (x, y) can be measured regularly (i.e., not at all times) by the GPS, and ψ is measured by the INS at all times.

To solve this problem, we will use the contractor approach described previously:

• First, we need to use the observation equation to contract the known GPS positions of the robot.

Remember also that any interval for which we do not have specific initial data first needs to be initialized as wide as possible to stay consistent with any further data, but − ∞ , ∞ should only be used as a last resort to minimize the risk that no contraction can occur. In practice, we can almost always assume a bounded position and speed due to the physical limitations of the system in its actual operating conditions.

The approach described in the previous section works also if we simulate the loss of the GPS position during a large portion of time (note that although the heading logs are here from the dual GPS, the system would have fell back to magnetometer in case of real GPS loss, which has no specific influence on the problem described in this section). During the loss of the GPS, the uncertainties will quickly increase over time (typical problem of dead reckoning), which translates to intervals (or boxes in 2D) with an increasing width, as explained in Le Bars et al. (2018). Figure 11 illustrates the trajectory estimation of the boat assuming the GPS was lost after 150 s, and details about the estimation of x and y are depicted in Figures 12 and 13 (see the Interval analysis section for more details about tubes). The figures have been generated in a computing time of around 30 s using the Codac library [see Rohou (2021)] and VIBes [see Drevelle and Nicola (2014)]. The GPS accuracy when it is in use is assumed to be 2.5 m, and its data are available every second. Additionally, a water current of 0.025 ; − 0.022 m/s (assumed to be uniform, although it is not always the case especially in a river) has been taken into account in the state equation model, and the coefficients α _f  = 1.3, α = 11.5, β = 23, and L = 1 have been estimated manually. Since most of the model uncertainties are likely to be related to the damping, it is assumed that α f ∈ 1.3 + − 0.02 ; 0.02 to ensure consistency of the model with the real positions. Therefore, the overall envelope of the trajectory estimation stays around a width of 5 m as long as the GPS is in use, but grows over time as soon as no position information is available any more. Despite the large uncertainty even at the last position, it is close to the GPS estimation if the center of the x and y intervals are taken as first approximation (which is often the most likely to be close to reality, similar to what a Kalman filter would do). During the period when the GPS is taken into account, it is assumed that the position it gives is the one at the time of reception of the data, i.e., this approximation is assumed to be within the width of the corresponding interval, which is 5 m. Due to this, it is allowed that the GPS trajectory (which is a step function) exits the position envelope, as long as it is at less than 2.5 m.

As mentioned in the Description of the autonomous boat section, a DVL sensor has been installed on the boat for some experiments. Therefore, a state estimation method inspired from, e.g., Le Bars et al. (2010) can be used (as a comparison, for a Kalman-based method, see, e.g., Reitbauer and Schmied (2021). There, a similar problem is simulated for a ground robot. Odometers in combination with a Kalman filter are used to get the vehicle speed). Indeed, if p = x y is the position of the robot, v _r is the speed vector measured directly by the DVL, and R is the Euler rotation matrix measured by the INS ¹¹ , we have: p ̇ = R ⋅ v r . (12)

We can assume that:

• For some time instant t ₀, we have a box p t 0 containing the state vector:

A propagation with respect to time on the differential Eq. 12 discretized using the explicit Euler method can be used to get an evaluation of the trajectory. Figure 14 illustrates the trajectory estimation of the boat assuming the GPS was lost after 20 s, and details about x and y estimations are visualized in Figures 15 and 16. Contrary to the previous situation without DVL, the position uncertainty grows slower when the GPS is not available anymore. There are also moments where the model is close to be inconsistent with the GPS position in Figure 15, although it is not a problem due to the GPS uncertainties as explained before ¹² .

Whatever the sensors or model used, these results demonstrate that interval methods can naturally provide a reliable estimation of the position error during the survey. The consistency with the GPS position considered as ground truth validates also that the interval-based INS used inside the boat outputs correct data. As long as the different hypotheses are true (e.g., assumptions on the precision of the discretized model, sensors data, etc.), the system can then theoretically provide a guaranteed explored, potentially explored, and surely not explored area as proposed in Desrochers and Jaulin (2016a) (here the position of the boat is in some way considered as the explored area, but in practice, the payload sensor characteristics should be used, and e.g., it would be the explored underwater area that would need to be guaranteed for the search of a target underwater).

Note that there were boat dynamic problems due to the DVL shape, the use of a smaller DVL such as Water Linked (2021), or improving its mechanical integration in the hull could easily limit the disturbances.