An Unscented Kalman Filter Based on the Adams–Bashforth Method with Applications to the State Estimation of Osprey-Type Drones Composed of Tiltable Rotor Mechanisms

TL;DR

This paper introduces a modified Unscented Kalman Filter using the Adams–Bashforth method to improve computational efficiency for state estimation in fast-moving drones.

Contribution

A novel application of the Adams–Bashforth method to the Unscented Kalman Filter for efficient state estimation in nonlinear systems.

Findings

The proposed method reduces computational load while maintaining high estimation accuracy.

Simulation results show improved performance for state estimation in an Osprey-type drone model.

The method is suitable for systems requiring high sampling frequencies.

Abstract

In the state estimation problem for nonlinear systems, the Unscented Kalman Filter (UKF) has gained attention as an algorithm capable of accurate state estimation based on high-fidelity discretization for strongly nonlinear systems. Furthermore, for applying the UKF to continuous-time state–space models, a method employing the Runge–Kutta method in the time-update equation for sigma points has already been proposed to achieve high-precision state estimation. While this method uses high-order numerical approximations, the associated decrease in computational efficiency due to processing time becomes problematic. It is thus unsuitable for the state estimation of relatively fast-moving objects, such as autonomous vehicles and drones, which require high sampling frequencies. In this study, to reduce computational load while achieving relatively high estimation accuracy, we newly apply the Adams–Bashforth method to the UKF algorithm. The effectiveness of the proposed method is demonstrated by first explaining a low-dimensional model’s state estimation problem, followed by a comparison of estimation accuracy and computation time in state estimation simulations for the UAV model of an Osprey-type drone.

Figures

Control allocation process. To generate control inputs for the actual vehicle, the intermediate variable vector $[eqn]$ is obtained using $[eqn]$ , and then the squared rotation speeds $[eqn]$ , $[eqn]$ and tilt angles $[eqn]$ , $[eqn]$ for each rotor are generated via the following relations:

Full text

1. Introduction

The navigation problem for moving objects is a fundamental challenge with important engineering applications [1]. The problem of position estimation is inherently formulated as a nonlinear filtering problem, and the theory of optimal filtering has been established by researchers such as [2]. However, in practice, only the linear Kalman filter (KF) proposed by Kalman [3] and its nonlinear extensions are feasible for implementation. This is because optimal nonlinear filters are computationally too expensive to implement.

The motion models of many real-world systems, particularly autonomous mobile robots, are described in continuous time. However, since observation information from sensors is obtained discretely through sampling, the observation model becomes discrete time. Such a system is called a continuous-discrete (CD) system [4]. While the approach that directly applies the Kalman filter to this CD system is referred to as CD-KF, in practice, it is common to discretize or numerically integrate the continuous-time motion model and construct a discrete-time Kalman filter. A Kalman filter in which both the motion model and the observation model are discrete-time is called a discrete–discrete Kalman filter (DD-KF). For nonlinear dynamics or nonlinear observation systems, the CD-EKF (continuous-discrete extended Kalman filter) or DD-EKF (discrete–discrete extended Kalman filter), which apply the Kalman filter to their linearly approximated systems, are widely used [5,6,7,8]. Furthermore, as a stable algorithm that does not require linear approximation via Jacobians, the Unscented Kalman Filter (UKF) has been proposed for DD nonlinear systems [9,10,11]. This will be referred to as DD-UKF hereafter. The UKF employs a nonlinear transformation called the Unscented Transform (UT). For a system of order n, it generates $[eqn]$ sample points called sigma points, applies a nonlinear transformation to each, and then calculates the conditional expectation of the state vector and the covariance matrix by taking the sample average of the transformed points. By using this UT, the UKF enables state estimation based on high-fidelity discretization.

Since the proposal of the UKF, many researchers have advanced its development and improvement. Särkkä [12] derived continuous-time and continuous-discrete versions of the UKF and applied them to nonlinear continuous-time filtering and re-entry vehicle tracking problems. van der Merwe [13] proposed the square-root UKF to improve numerical stability, and Ito and Xiong [14] analyzed the convergence of Gaussian filters, including the UKF. Julier [15] proposed the spherical simplex UT to reduce the number of sigma points. Meanwhile, Haykin [16] advanced research on the integration of neural networks and Kalman filtering, demonstrating new possibilities for nonlinear adaptive filtering. Elliott [17] contributed to the foundations of filtering theory through the theory of hidden Markov models and stochastic control. Recently, Singh [18] comprehensively reviewed the development of Gaussian filtering, including the UKF. In particular, recent studies have combined the UKF with outlier handling and event-triggered mechanisms. For example, Ref. [19] proposed a dual-event-triggering ANFIS-based UKF for cluster cooperative navigation with measurement anomalies, demonstrating potential applications to navigation problems in complex environments.

From the perspective of numerical integration methods, Takeno and Katayama [20] applied the Runge–Kutta (RK) method to the time-update equations for sigma points and demonstrated that high-precision state estimation is possible for several models. Furthermore, Takeno and Katayama [21] applied Heun’s method, an improved version of the Euler method, to the prediction step. Kulikov and Kulikova [22,23] have advanced research on high-precision numerical implementation of the continuous-discrete UKF. However, discussions on these UKFs have focused primarily on estimation accuracy, with little consideration given to computational efficiency. Due to its characteristic of performing nonlinear transformations on each of the $[eqn]$ sigma points in the UT algorithm, the UKF has the drawback of increased computation time as the system order n grows. For state estimation of large-scale real-world models, both accuracy and real-time performance are required, making the improvement of computational efficiency in UKF algorithms critically important [24]. In recent years, research addressing this computational efficiency challenge has begun to emerge. Wang et al. [25] proposed an adaptive step-size UKF based on the degree of nonlinearity, achieving a balance between accuracy and computational cost by dynamically adjusting the step size in highly nonlinear regions. Researchers in the aerospace domain [26] introduced a high-efficiency UKF utilizing parallel computation for sigma point propagation, demonstrating significant speedup for multi-target trajectory estimation. The integration of machine learning techniques with UKF has also been explored; a GAN-enhanced UKF framework [27] dynamically predicts and updates filter parameters in real-time, improving estimation accuracy without sacrificing computational performance. In the realm of numerical implementation, Kulikova and Kulikov [28] developed square-root information-type methods for continuous-discrete extended Kalman filtering, enhancing numerical stability which is crucial for efficient computation. Furthermore, theoretical advances in sensor scheduling for continuous-discrete systems [29] provide insights into optimizing the trade-off between resource allocation and estimation accuracy. A comparative evaluation of nonlinear filters [30] further contextualizes the performance of the UKF against other methods in practical tracking applications. These recent developments underscore the growing recognition of computational efficiency as a critical factor in the practical deployment of nonlinear Kalman filters.

Regarding Unmanned Aerial Vehicle (UAV) navigation, comprehensive reviews exist on topics such as autonomous navigation in GPS-denied environments [31], UAV positioning using GNSS [32], and sensor-based autonomous localization [33]. These reviews highlight the importance of filtering techniques in UAV navigation. Therefore, in this study, to improve the computational efficiency of the UKF, we propose a UKF that applies the Adams–Bashforth (AB) method, instead of the RK method, to the time-update equations for sigma points. The RK method is a single-step numerical integration method; it calculates the next numerical point $[eqn]$ based solely on the current point $[eqn]$ . Furthermore, when the order of accuracy is $[eqn]$ , s function evaluations are required to generate the next numerical point. On the other hand, the AB method is a multi-step numerical integration method. For an accuracy order of 4, the next numerical point $[eqn]$ is calculated based on the past four points $[eqn]$ , and only one function evaluation is needed to generate the next numerical point. Thus, regardless of the accuracy order s of the AB method, the number of function evaluations is always one. Therefore, compared to the fixed-step RK method, computation with the AB method is more efficient.

It should be noted, as a related study focusing on the degradation of computational efficiency of the RK method when applying numerical integration to nonlinear Kalman filters, that He et al. [34] have already conducted research applying a multi-step method based on the Adams–Bashforth–Moulton (ABM) method to the EKF. The objective of this study is to demonstrate that the UKF applying the AB method to the time-update equations for sigma points in the UT can achieve estimation with computational efficiency superior to that of the RK method, while maintaining comparable estimation accuracy. In addition to comparison with the RK method, we also compare the differences in estimation accuracy and computation time resulting from varying orders (2nd to 6th) of the AB method itself.

The effectiveness of the proposed method is demonstrated through state estimation in a MATLAB 2022a simulation environment using two nonlinear models. Specifically, as a preliminary experiment, we first apply the method to a low-dimensional falling object model used in previous research and compare estimation performance under identical conditions. Subsequently, as an application to a more complex model, we perform a similar comparison of estimation performance using a UAV model of an Osprey-type drone.

The structure of this paper is as follows. Section 2 formulates the continuous-discrete system and the state estimation problem. Section 3 discusses the integration of the UT with the Runge–Kutta method and the Adams–Bashforth method, presenting the time-update equations for sigma points based on each method. Section 4 details the modeling of the Osprey-type drone. After defining the coordinate systems, the rotational and translational motions are explained in detail. Section 5 describes the control law and control allocation. Specifically, a controller based on the computed torque method and the control input allocation problem are detailed. Section 6, as a preliminary experiment, applies the sigma-point time-update equations presented in Section 3 to the falling object model used in previous research, comparing estimation accuracy and computation time. Section 7 first derives the state equations and sigma-point time-update equations for the Osprey-type UAV model, then performs comparisons of estimation accuracy and computation time similar to those in Section 6. Furthermore, we provide discussions based on a broader range of estimation results, including not only comparisons with the RK method but also differences in estimation performance due to the order of the AB method and comparisons of estimation performance under different sampling periods. Section 6 summarizes the paper and presents concluding remarks.

2. Continuous-Discrete System and Problem Formulation

2.1. Continuous-Discrete System

The motion of many real-world systems, particularly moving objects such as UAVs, is described in continuous time, while observations from sensors are obtained discretely through sampling. Such systems are called CD systems, and the foundations of their filtering theory are detailed in Jazwinski [4].

In this study, we consider a system described by the following nonlinear continuous-time state equation and discrete-time observation equation:

$$[eqn]$$ $$[eqn]$$

Here, t denotes continuous time, $[eqn]$ ( $[eqn]$ ) denotes sampling instants, and $[eqn]$ is the sampling interval. $[eqn]$ is the state vector, $[eqn]$ is the input vector, and $[eqn]$ is the output vector. $[eqn]$ and $[eqn]$ are nonlinear functions. Furthermore, $[eqn]$ is the system noise, and $[eqn]$ is the observation noise, both assumed to be Gaussian white noises with zero mean and covariance matrices $[eqn]$ and R, respectively.

2.2. Discrete–Discrete Unscented Kalman Filter

Discretizing the continuous-time state Equation (1) using some numerical integration method yields the following discrete-time state equation:

$$[eqn]$$

where $[eqn]$ is an approximation function obtained through the numerical integration method, and $[eqn]$ is the system noise arising from discretization. The UKF can be applied to the DD system composed of Equations (2) and (3).

The UKF is an algorithm that performs state estimation for nonlinear systems using the Unscented Transform (UT) [9,10,11]. Hereafter, the filtered estimate is denoted as $[eqn]$ , its covariance matrix as $[eqn]$ , the one-step predicted estimate as $[eqn]$ , and its covariance matrix as $[eqn]$ . The initial estimate $[eqn]$ and initial covariance matrix $[eqn]$ are assumed to be known.

The UKF algorithm consists of a prediction step and an update step. The detailed derivation and specific equations are provided in Appendix A; here, we present only an overview.

Prediction Step: From the estimate $[eqn]$ and covariance matrix $[eqn]$ at time $[eqn]$ , $[eqn]$ sigma points $[eqn]$ and corresponding weights $[eqn]$ are generated. These sigma points are propagated through the nonlinear function $[eqn]$ , and the one-step predicted estimate $[eqn]$ and its covariance matrix $[eqn]$ are computed by taking the weighted average.

Update Step: The predicted sigma points $[eqn]$ are substituted into the observation function $[eqn]$ to compute the predicted output $[eqn]$ , its covariance $[eqn]$ , and the cross-covariance between state and output $[eqn]$ . The innovation $[eqn]$ is computed using the actual observation $[eqn]$ , and the estimate $[eqn]$ and covariance matrix $[eqn]$ are updated via the Kalman gain $[eqn]$ .

For the specific computational procedure and complete set of equations for the UKF, please refer to Appendix A.

3. Unscented Transform and Discretization Methods

3.1. Discrete-Time State–Space Model

Consider the following continuous-time nonlinear system:

$$[eqn]$$

where $[eqn]$ is the state vector and $[eqn]$ is the initial value. Discretizing this continuous-time system using some numerical integration method yields the following time-update equation:

$$[eqn]$$

The observation equation is given similarly to Equation (2) as

$$[eqn]$$

where v is Gaussian white noise with zero mean and covariance matrix R. Equations (5) and (6) constitute the discrete-time state–space model.

The objective of this study is to propose a nonlinear filtering algorithm that sequentially estimates the state vector $[eqn]$ based on the observation data $[eqn]$ . Takeno and Katayama [20] demonstrated that combining the UKF with the RK method enables high-precision state estimation. In this study, we show that by combining the AB method with the UT instead of the RK method, relatively high-precision state estimation can be achieved while reducing computational load.

For the general theory of numerical solutions for stochastic differential equations, please refer to Kloeden and Platen [35].

3.2. Integration of UT and the Runge–Kutta Method

When applying the RK method to the time update of sigma points in the UKF, numerical integration is performed for each sigma point. The sigma point matrix $[eqn]$ at time k with its element $[eqn]$ is defined as follows:

$$[eqn]$$

Here, $[eqn]$ denotes the element of the state vector, and $[eqn]$ denotes the index of the sigma point. Let the step size be $[eqn]$ and the time instant be $[eqn]$ .

When applying the 4th-order Runge–Kutta method, the time-update equation for the sigma points is expressed in the following general form:

$$[eqn]$$

Here, $[eqn]$ are the function evaluations at each stage of the RK method. The specific definitions and detailed derivations are provided in Appendix B.

By taking the weighted average of these updated sigma points, the predicted value of the state vector (Equation (A10)) and the predicted covariance matrix (Equation (A11)) are obtained.

3.3. Integration of UT and the Adams–Bashforth Method

This section describes the integration of the AB method with the UT, which is the core of this study. The ordinary AB method is a multi-step method that calculates the next point using information from multiple past points. In this study, we extend and apply this method to the sigma point matrix of the UKF.

First, the result of substituting each column (each sigma point) of the sigma point matrix $[eqn]$ at time k into the nonlinear function $[eqn]$ is defined as the following $[eqn]$ matrix:

$$[eqn]$$

This $[eqn]$ is a matrix that stores the results of evaluating the nonlinear function for all sigma points at time k. While the ordinary AB method uses function evaluation values for past state vectors, the feature of the proposed method is that it stores function evaluation matrices for past sigma point matrices $[eqn]$ and uses them as “past information.”

Below, we present the time-update equations for sigma points using the 2nd- to 6th-order AB methods [36,37].

2nd-order Adams–Bashforth method:

$$[eqn]$$

3rd-order Adams–Bashforth method:

$$[eqn]$$

4th-order Adams–Bashforth method:

$$[eqn]$$

5th-order Adams–Bashforth method:

$$[eqn]$$

6th-order Adams–Bashforth method:

$$[eqn]$$

Here, $[eqn]$ and $[eqn]$ are $[eqn]$ matrices. These equations are an extension of the ordinary AB method formulas to the sigma point matrix, representing an original formulation for integrating the numerical integration method into the specific filtering method of UKF. In particular, the method of storing and utilizing the results of function evaluations for past sigma point matrices is a novel aspect of this approach not found in existing research.

3.4. Comparison of Computational Complexity: Runge–Kutta Method vs. Adams–Bashforth Method

We compare the computational complexity of the RK and AB methods in the UT from the following two perspectives.

3.4.1. Computational Complexity Comparison Based on “Number of Stages”

The difference width $[eqn]$ in discretization is generally called the step size in numerical computation. The RK method computes the derivative $[eqn]$ multiple times while advancing one step (this is referred to as the “number of stages”). For example, the 4th-order RK method requires four stages of computation, and its global truncation error is $[eqn]$ .

On the other hand, since the AB method utilizes previously computed sigma point matrices $[eqn]$ , only one computation of $[eqn]$ is sufficient to derive the time-update equation for the current sigma points. That is, regardless of the accuracy order of the AB method, the number of stages is always $[eqn]$ that of the RK method.

3.4.2. Computational Complexity Comparison Based on “Number of Elements in the Sigma Point Set”

Generally, when simulating large-scale models in real space, the number of state variables n can become very large (millions to tens of millions). In the UKF, since the size of the sigma point set is $[eqn]$ , this effect becomes pronounced. The RK method performs four stages of derivative computation for each sigma point, for each element of the n-dimensional state vector and for each element of the sigma point set.

In contrast, the AB method requires only one stage of function computation for each element of the sigma point set. Furthermore, for past sigma point sets, it only involves multiplying the entire matrix by coefficients according to the accuracy order. Therefore, the impact of an increase in the number of state variables is significantly smaller compared to the RK method.

This theoretical comparison of computational complexity will be verified through numerical experiments presented later in Section 6 and Section 7.

4. Modeling of the Osprey-Type Drone

4.1. Definition of Coordinate Systems

Figure 1 shows the coordinate systems used in this study for the Osprey-type drone. The world coordinate system $[eqn]$ is a right-handed coordinate system with origin $[eqn]$ and axes $[eqn]$ , $[eqn]$ , $[eqn]$ , where $[eqn]$ is positive vertically downward. The body coordinate system $[eqn]$ has its origin $[eqn]$ at the center of gravity of the vehicle. It is also a right-handed coordinate system with axes $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ positive vertically downward. The positive direction of the $[eqn]$ axis is designated as the forward direction of the vehicle.

The coordinate system for the first coaxial rotor, $[eqn]$ , has origin $[eqn]$ and axes $[eqn]$ , $[eqn]$ , $[eqn]$ . Similarly, the coordinate system for the second coaxial rotor, $[eqn]$ , has origin $[eqn]$ and axes $[eqn]$ , $[eqn]$ , $[eqn]$ . The coordinate system for the i-th ( $[eqn]$ ) coaxial rotor is integrated into $[eqn]$ . In the body coordinate system $[eqn]$ , the angular velocities about the $[eqn]$ , $[eqn]$ , $[eqn]$ axes are $[eqn]$ , respectively. Furthermore, in the rotor coordinate system $[eqn]$ , the tilt angles about the $[eqn]$ , $[eqn]$ axes are $[eqn]$ . The initial tilt angles $[eqn]$ and $[eqn]$ are 0. Due to servo motor characteristics, the range of $[eqn]$ is $[eqn]$ , while $[eqn]$ is unlimited.

Let the Euler angles in the world coordinate system be $[eqn]$ . The rotation matrices about the x, y, z axes, $[eqn]$ , $[eqn]$ , $[eqn]$ , can be expressed as follows (S denotes sin, and C denotes cos):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

From these $[eqn]$ , $[eqn]$ , $[eqn]$ , the transformation matrix from the body coordinate system to the world coordinate system, $[eqn]$ , can be expressed as

$$[eqn]$$

Also, let the transformation matrix for angular velocity from the world coordinate system to the body coordinate system be $[eqn]$ . Then, $[eqn]$ is

$$[eqn]$$

and its inverse is given by

$$[eqn]$$

Here, it is found in [38,39] that

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ . Note that $[eqn]$ is invertible unless $[eqn]$ . That is, it is invertible as long as it does not take the specific $[eqn]$ where gimbal lock occurs.

From the tilt angles $[eqn]$ and $[eqn]$ , where $[eqn]$ in each rotor coordinate system, the transformation matrix from each rotor coordinate system to the body coordinate system, $[eqn]$ , can be expressed as

$$[eqn]$$ $$[eqn]$$

4.2. Rotational Motion

First, the angular velocity $[eqn]$ and angular acceleration $[eqn]$ of the i-th ( $[eqn]$ ) coaxial rotor can be expressed as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ are the angular velocities of each brushless motor in the first coaxial rotor, and $[eqn]$ , $[eqn]$ are those of the second coaxial rotor.

Next, the reaction torque $[eqn]$ of the i-th ( $[eqn]$ ) coaxial rotor can be expressed as

$$[eqn]$$

where $[eqn]$ represents the torque coefficient.

Also, the thrust $[eqn]$ of the i-th ( $[eqn]$ ) coaxial rotor can be expressed as

$$[eqn]$$

where $[eqn]$ represents the thrust coefficient of the coaxial rotor.

Using the Newton–Euler method, the torque $[eqn]$ generated by the coaxial rotor can be expressed as

$$[eqn]$$

Here, since the angular velocity and angular acceleration caused by the servo motor tilt are instantaneous and minute, they can be expressed as

$$[eqn]$$

Therefore, $[eqn]$ and $[eqn]$ become as follows:

$$[eqn]$$ $$[eqn]$$

4.2.1. Rotational Torque τB Generated by Rotor Thrust

The rotational torque $[eqn]$ generated on the vehicle by rotor thrust can be expressed using the position vector $[eqn]$ of each rotor in the body coordinate system as follows:

$$[eqn]$$ $$[eqn]$$

The calculation details for $[eqn]$ are shown below:

$$[eqn]$$ $$[eqn]$$

Therefore, $[eqn]$ can be expressed as

$$[eqn]$$

4.2.2. Reaction Torque τc Generated by Rotor Rotation

The reaction torque $[eqn]$ generated by rotor rotation can be expressed using the reaction torque $[eqn]$ of each rotor as follows:

$$[eqn]$$

The calculation details for $[eqn]$ are shown below:

$$[eqn]$$ $$[eqn]$$

Therefore, $[eqn]$ can be expressed as

$$[eqn]$$

4.2.3. Vehicle Coriolis Force τcori

The vehicle Coriolis force $[eqn]$ can be expressed using the vehicle’s angular velocity $[eqn]$ and its inertia matrix $[eqn]$ as follows:

$$[eqn]$$

Here, since

$$[eqn]$$

we have

$$[eqn]$$

Therefore, letting $[eqn]$ , we obtain

$$[eqn]$$

4.2.4. Vehicle Angular Acceleration ω˙B

Using the equations above, the vehicle’s angular acceleration $[eqn]$ can be expressed as

$$[eqn]$$

Furthermore, neglecting external disturbance torques $[eqn]$ acting on the vehicle and gyroscopic effects (i.e., $[eqn]$ ), we have

$$[eqn]$$

4.3. Translational Motion

The thrust $[eqn]$ in the body coordinate system can be expressed using the thrust $[eqn]$ of each rotor as

$$[eqn]$$

Therefore, letting $[eqn]$ , we obtain

$$[eqn]$$

Using the vehicle mass m and gravitational acceleration $[eqn]$ , the vehicle’s translational acceleration $[eqn]$ can be expressed as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the external disturbance force acting on the vehicle. Neglecting effects such as friction, we set $[eqn]$ .

5. Controller Design and Control Allocation

5.1. Computed Torque Method

The computed torque method is a technique to determine the generalized force (comprising thrust and torque in the body coordinate system) from the desired translational and angular accelerations in the world coordinate system, given that the vehicle’s physical information is known. Utilizing the relation $[eqn]$ in the translational acceleration of Equation (44) and the rotational angular acceleration of Equation (40), and defining the generalized coordinates as $[eqn]$ , the system can be represented as

$$[eqn]$$

This is because, from Equation (20),

$$[eqn]$$

and thus

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Here, since the rotation matrix $[eqn]$ is orthogonal, $[eqn]$ . Also, as each principal moment of inertia $[eqn]$ , $[eqn]$ . Then, the inverse system of Equation (46) is

$$[eqn]$$

Let the target value for the generalized coordinate vector X be $[eqn]$ . Constructing its augmented acceleration vector $[eqn]$ using a PD servo system yields

$$[eqn]$$

where $[eqn]$ , $[eqn]$ is the derivative gain, and $[eqn]$ is the proportional gain. Substituting $[eqn]$ from Equation (49) for $[eqn]$ in Equation (48) gives

$$[eqn]$$

Substituting this $[eqn]$ into the original plant Equation (46) yields

$$[eqn]$$

Here, we set $[eqn]$ and $[eqn]$ . For application to actual vehicles, considering modeling errors, the condition $[eqn]$ , which ideally achieves a damping ratio of 1 (i.e., critical damping), is used as a base, but the damping condition for this error may need slight modification. In the simulations, the proportional gains for translational position and attitude angles are defined as $[eqn]$ , and the derivative gains as $[eqn]$ , with the following settings:

$$[eqn]$$

5.2. Control Input Allocation Problem

Appropriate rotor speeds and tilt angles are allocated to each coaxial rotor to achieve the generalized force in the body coordinate system, determined earlier by the computed torque method. A single coaxial rotor constructed in this study utilizes two brushless motors. The thrust of coaxial rotor 1 is $[eqn]$ , and that of coaxial rotor 2 is $[eqn]$ . The reaction torque for the first rotor is $[eqn]$ , and for the second rotor is $[eqn]$ . However, directly solving for the rotation speeds $[eqn]$ and tilt angles $[eqn]$ results in an underdetermined problem: six generalized forces versus eight unknown decision variables. This would yield an effectiveness matrix of size $[eqn]$ , requiring solution via pseudo-inverse or numerical optimization such as constrained QP. Below, we present a simpler method based on experimental results for coaxial rotor thrust.

Referring to the experimental results on coaxial rotors by Itakura et al. [40], for a total thrust of $[eqn]$ N of the coaxial rotor, the thrust of the upstream propeller alone was $[eqn]$ N, and that of the downstream propeller alone was $[eqn]$ N. Therefore, the thrust efficiency is

$$[eqn]$$

This serves as an indicator of deviation from theoretical ideal conditions. On the other hand, the interference efficiency, which is the increase rate of total thrust relative to the upstream propeller’s thrust, is

$$[eqn]$$

Introducing this interference efficiency allows for the substitution of $[eqn]$ with $[eqn]$ and $[eqn]$ with $[eqn]$ , thereby eliminating $[eqn]$ and $[eqn]$ from the unknown variables.

With this interference efficiency, the reaction torque $[eqn]$ from Equation (37) becomes

$$[eqn]$$

Meanwhile, the rotational torque due to thrust $[eqn]$ from Equation (35) becomes

$$[eqn]$$

Similarly, the thrust $[eqn]$ from Equation (43) is rewritten as

$$[eqn]$$

Therefore, by introducing an intermediate 6-dimensional vector $[eqn]$ , defining $[eqn]$ , $[eqn]$ , and for $[eqn]$ ,

$$[eqn]$$

the relation $[eqn]$ finally reduces to a linear equation:

$$[eqn]$$

Here, the effectiveness matrix A is

$$[eqn]$$

Figure 2. Control allocation process. To generate control inputs for the actual vehicle, the intermediate variable vector $[eqn]$ is obtained using $[eqn]$ , and then the squared rotation speeds $[eqn]$ , $[eqn]$ and tilt angles $[eqn]$ , $[eqn]$ for each rotor are generated via the following relations:

$$[eqn]$$

Figure 2 illustrates the sequential calculation flow for this control allocation problem.

6. Preliminary Simulation: Application to the Falling Body Model

In this section, we describe state estimation for the falling body model [41] as a preliminary experiment. The falling body model is often used as a nonlinear model and was also employed as a simulation model by Takeno and Katayama [20]. The purpose here is to verify the basic effectiveness of the proposed method by performing state estimation under conditions similar to previous studies. Since the RK method used by Takeno et al. is a 4th-order discretization method, this section focuses on the 4th-order AB method for comparison.

6.1. Model Overview (Falling Body)

The differential equations for the falling body model are as follows [41]:

$$[eqn]$$ $$[eqn]$$

Here, $[eqn]$ is the altitude of the falling object, g is the gravitational acceleration, $[eqn]$ is the atmospheric density at sea level (0 m), $[eqn]$ is the decay coefficient, and $[eqn]$ is the ballistic coefficient.

The extended state vector is defined as follows:

$$[eqn]$$

Letting $[eqn]$ , the state–space model becomes

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The observation equation is nonlinear and is assumed to be as follows:

$$[eqn]$$

Here, $[eqn]$ is the horizontal distance between the falling object and the radar, and $[eqn]$ is the radar altitude. Furthermore, v is the observation noise, assumed to be Gaussian white noise with zero mean and variance R.

6.2. Estimation Algorithms (Falling Body)

In this study, we compare three types of numerical integration methods integrated with the UKF: the Euler method, the RK method, and the AB method. Hereafter, the UKF using the Euler method is called “Euler-UKF”, that using the RK method “RK-UKF”, and that using the AB method “AB-UKF”.

The detailed derivation and specific equations for the sigma point time update of each method are presented in Appendix C. This section provides only an overview of each method.

6.2.1. Euler-UKF (Falling Body)

In Euler-UKF, the Euler method is applied to Equations (64)–(66). The Euler method has first-order accuracy and requires one function evaluation per update step. The specific sigma point time-update equations are given in Appendix C.1.

6.2.2. RK-UKF (Falling Body)

In RK-UKF, the 4th-order RK method (hereafter denoted as RK4) is applied to Equations (64)–(66) [20]. The RK4 method has fourth-order accuracy and requires four function evaluations per update step. The specific sigma point time-update equations are given in Appendix C.2.

6.2.3. AB-UKF (Falling Body)

In AB-UKF, the 4th-order AB method (hereafter denoted as AB4) is applied to Equations (64)–(66). The AB4 method achieves fourth-order accuracy while requiring only one function evaluation per update step. Instead, it stores and utilizes the function evaluation results for past sigma point matrices, $[eqn]$ , $[eqn]$ , and $[eqn]$ . The specific sigma point time-update equations are given in Appendix C.3.

6.3. Estimation Conditions (Falling Body)

State estimation simulations are performed using the three methods described above. The physical parameters are as follows: $[eqn]$ [m/s^2^], $[eqn]$ [kg·s^2^/m^4^], $[eqn]$ [m]. The true ballistic coefficient is $[eqn]$ [kg/m^2^]. The horizontal distance between the falling object and the radar is $[eqn]$ = 10,000 [m], and the radar altitude is $[eqn]$ [m]. The step size is $[eqn]$ [s], the total number of sampling points is $[eqn]$ , and the UT parameter is $[eqn]$ . The initial state vector, system noise covariance, and observation noise covariance are as follows [20]:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The initial state estimate and its covariance matrix are as follows:

$$[eqn]$$ $$[eqn]$$

6.4. Simulation Results

6.4.1. State Estimation Accuracy

Figure 3a shows the estimation results for the ballistic coefficient $[eqn]$ , and Figure 3b shows the time evolution of the state estimation error. The state estimation error was calculated using the following RMSE (Root Mean Squared Error) formula:

$$[eqn]$$

The results are averages from 100 simulation runs for each method.

Figure 3a shows that the parameter $[eqn]$ converges well for all three methods: Euler-UKF, RK-UKF, and AB-UKF. On the other hand, Figure 3b indicates that after the convergence of $[eqn]$ , the state estimation error is suppressed more effectively in the order of Euler-UKF, RK-UKF, and AB-UKF.

Table 1 presents the average RMSE values of the state estimation error obtained from 100 simulation runs for each of the three methods. The computations were performed on a laptop with 16 GB RAM, an AMD Ryzen7 4800H with Radeon Graphics CPU, running Windows 10 and MATLAB 2022a.

As shown in Table 1, the Euler method exhibits the largest RMSE, while the RK and AB methods show smaller RMSE values than the Euler method. This result is attributed to the difference in discretization order: the Euler method has first-order accuracy, whereas the RK and AB methods have fourth-order accuracy.

6.4.2. Computation Time

To more clearly demonstrate the differences between the methods, the estimation conditions from Section 6.3 were modified: $[eqn]$ and the total number of sampling points $[eqn]$ . Table 2 shows the comparison results for the total computation time of the entire UKF algorithm (from Equations (A3) to (A18)) and the computation time for the prediction step only (Equation (A9)).

Table 2 shows that the computation time is shortest for the Euler method, followed by the AB method, and then the RK method. The particularly short prediction time of the Euler method stems from the difference in discretization order.

6.5. Discussion

The results presented in Section 6.4 can be attributed to the following two factors.

Small number of state variables: The falling body model has three state variables, with $[eqn]$ representing a parameter estimation problem. Consequently, as is evident from the sigma point time-update equations, the discretization formulas are actually applied only to the state variables $[eqn]$ and $[eqn]$ , making this a model where differences between the RK and AB methods are less likely to manifest.Small proportion of UT within the overall UKF algorithm: The falling body model has a relatively small number of state variables and simple discretization equations, making it difficult for differences in computational load to become pronounced. As a result, the time spent on UT constitutes a small proportion of the overall UKF estimation algorithm, and the differences in computation time due to the discretization method are less noticeable.

Considering the above two points, the next section, Section 7, addresses a state estimation problem for a UAV model, which is higher dimensional than the falling body model.

7. Application to the Osprey-Type Drone

7.1. Model Overview (Drone)

The UAV vehicle conditions are based on the Osprey-type drone with 2-DOF tiltable coaxial rotors by Itakura et al. [40]. The vehicle’s general appearance and the derivation of its dynamic model have already been detailed in Section 4.

The vehicle’s position in the world coordinate system, $[eqn]$ , follows the translational motion equation expressed by Equations (44) or (45).

On the other hand, the vehicle’s angular acceleration $[eqn]$ is given by Equations (40) or (41). Assuming the Euler angle variations are relatively small, i.e., $[eqn]$ , then $[eqn]$ . Therefore, Equation (41) can be rewritten as

$$[eqn]$$

Now, defining a 12-dimensional state vector as $[eqn]$ , the continuous-time state equation for the UAV model can be derived from Equations (45) and (74) as follows:

$$[eqn]$$

Here, $[eqn]$ in Equation (75) are defined as follows:

$$[eqn]$$

7.2. Estimation Algorithms (Drone)

Similar to the falling body model, we compare three types of numerical integration methods integrated with the UKF for the UAV model: the Euler method, the RK method, and the AB method. The detailed derivation and specific equations for the sigma point time update of each method are presented in Appendix D. This section provides only an overview of each method.

7.2.1. Euler-UKF (Drone)

In Euler-UKF, the Euler method is applied to Equation (75). The Euler method has first-order accuracy and requires one function evaluation per update step. The specific sigma point time-update equations are given in Appendix D.1.

7.2.2. RK-UKF (Drone)

In RK-UKF, the RK4 method is applied to Equation (75). The RK4 method has fourth-order accuracy and requires four function evaluations per update step. The specific sigma point time-update equations are given in Appendix D.2.

7.2.3. AB-UKF (Drone)

In AB-UKF, the AB method is applied to Equation (75). Regardless of its order, the AB method requires only one function evaluation per update step and stores and utilizes the function evaluation results for past sigma point matrices. The specific forms of the sigma point time-update equations corresponding to the 2nd to 6th order AB methods are given in Appendix D.3.

7.3. Estimation Conditions (Drone)

This section describes the UAV physical parameters and the covariance of system and observation noise used in the estimation simulation.

7.3.1. Parameter Settings for the UAV Model

Similar to the falling body model, state estimation simulations are performed using the three methods described above. Table 3 summarizes the definitions and values of each parameter.

All initial states are set to $[eqn]$ , and the initial estimation information is set to $[eqn]$ and $[eqn]$ .

7.3.2. Determination of System Noise

In this simulation experiment, the system noise is set based on the global truncation error arising from discretization by the Euler, RK, and AB methods. When the step size is denoted by h, the global truncation error for the Euler method is proportional to $[eqn]$ , and for the RK4 and AB4 methods, it is proportional to $[eqn]$ . Therefore, letting the variance values of the respective system noises be $[eqn]$ and $[eqn]$ , they can be expressed as

$$[eqn]$$ $$[eqn]$$

Thus, for comparison purposes in this simulation experiment, the value $[eqn]$ is set as follows:

$$[eqn]$$

Here, the exponent $[eqn]$ signifies taking an intermediate order of magnitude between $[eqn]$ and $[eqn]$ .

7.3.3. Determination of Observation Noise

Typically, UAVs are equipped with GPS sensors for position observation and gyro sensors for attitude observation. Determining a precise noise value for GPS sensor error is difficult due to various influencing factors; thus, for simplicity, the noise value is set to 2 m here. For the gyro sensor, we assume the use of the MPU-9250, a 9-axis sensor module from InvenSense.

The MPU-9250 datasheet states a gyro sensor noise value of 0.01 [deg/ $[eqn]$ ]. Using this noise value and the step size h [s], the observation noise variance $[eqn]$ [rad^2^] for attitude angles $[eqn]$ and $[eqn]$ [rad^2^/s^2^] for attitude angular velocities $[eqn]$ are determined.

First, convert the dimension 0.01 [deg/ $[eqn]$ ] to [rad/ $[eqn]$ ] to obtain

$$[eqn]$$

Squaring the value obtained in Equation (80) yields

$$[eqn]$$

Here, the dimension [Hz] can also be expressed as [1/s]; thus, the dimension [rad^2^/Hz] can also be expressed as [rad^2^·s]. Therefore, the observation noise variance $[eqn]$ [rad^2^] for attitude angles $[eqn]$ is obtained by dividing Equation (81) by h [s]:

$$[eqn]$$

Then, the observation noise variance $[eqn]$ [rad^2^/s^2^] for attitude angular velocities $[eqn]$ is obtained by dividing Equation (82) by $[eqn]$ [s]:

$$[eqn]$$

Consequently, the observation equation is

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the observation matrix for position $[eqn]$ :

$$[eqn]$$

The covariance $[eqn]$ of the observation noise $[eqn]$ is represented by the following diagonal matrix:

$$[eqn]$$

Remark

The observation setting used in this paper is linear (direct observation of position and attitude angles), and the advantages of nonlinear filtering are not fully utilized. In actual UAV navigation, nonlinear observations such as radar observations in GPS-denied environments are common. While the main focus of this paper is the analysis of the impact of state equation complexity on computational efficiency, and the linearity of observations is not an essential problem, extension to nonlinear observations is an important future task.

7.4. Target Trajectory

The simulation starts from the initial position $[eqn]$ and initial attitude angles $[eqn]$ . First, over 10 s from the start time, the vehicle moves to the target position $[eqn]$ , while all target attitude angles remain zero. Subsequently, from 10 s to 70 s, the vehicle is commanded to follow the target position, velocity, and acceleration given below [42].

**Target position $[eqn]$ **:

$$[eqn]$$

**Target velocity $[eqn]$ **:

$$[eqn]$$

**Target acceleration $[eqn]$ **:

$$[eqn]$$

This target trajectory is a circular path with a radius of 1 m and center $[eqn]$ in the $[eqn]$ -plane, rotated 45 deg about the x-axis, and translated −1 m along the z-axis. Since this vehicle is left-right symmetric, independent control in three degrees of freedom (x, y, z directions) can be verified by following this circular path. The vehicle follows this trajectory with a period of 20 s per revolution. Simultaneously, the target roll angle changes every revolution: from 0 deg to 45 deg, 45 deg to 0 deg, and 0 deg to −45 deg. The target pitch angle changes every revolution: from 0 deg to −45 deg, −45 deg to 0 deg, and 0 deg to 45 deg. The target yaw angle is always 0 deg. Following this trajectory verifies the possibility of independent control in all six degrees of freedom.

7.5. Description of Comparative Experiments

This comparative experiment focuses on two aspects of the UKF based on each discretization method: “estimation accuracy” and “computational efficiency.” Details of the comparison methods are described below. Comparisons are performed in two patterns: “comparison of three methods: Euler, RK4, and AB4” and “comparison of 2nd to 6th-order AB methods.” The results of this comparative experiment are averages from 50 Monte Carlo simulations.

7.5.1. Estimation Accuracy Comparison Experiment

The estimation accuracy comparison is performed by calculating the RMSE values between the estimates from each discretization-based UKF and the true values.

The RMSE value at estimation time k ( $[eqn]$ ) is calculated using the following equations:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the true value and the estimate of position x at time i, respectively. Here, N is the total number of sampling points, determined from the simulation time of 70 s and the discretization period h [s] as $[eqn]$ . For example, if $[eqn]$ [s], then $[eqn]$ . Equations (91) and (92) are similarly applied to other variables: positions $[eqn]$ , velocities $[eqn]$ , attitude angles $[eqn]$ , and attitude angular velocities $[eqn]$ .

Furthermore, when comparing estimation accuracy under different discretization periods ( $[eqn]$ ), the sum of the individual RMSE values $[eqn]$ $[eqn]$ for the 12 state variables at the final discrete time $[eqn]$ is denoted as $[eqn]$ , and this value is used for comparison.

Here, we discuss the theoretical meaning of $[eqn]$ . In the UAV model, when the a priori error covariance matrix $[eqn]$ at time $[eqn]$ is expressed as

$$[eqn]$$

the trace of this error covariance matrix is denoted as $[eqn]$ . The aforementioned $[eqn]$ is the sum of RMSE values based on actual sampled estimates, averaged over time up to $[eqn]$ and averaged over the ensemble of 50 simulation results. Therefore, while transient values differ from the actual $[eqn]$ , they can be treated as equivalent values in the steady state.

7.5.2. Computational Efficiency Comparison Experiment

The computational efficiency comparison experiment measures and compares computation times at the following three points:

• Measurement Time I: Computation time for sigma point time update (computation of Equation (A9)).
• Measurement Time II: Computation time for the entire UKF algorithm (computation from Equations (A3) to (A18)).
• Measurement Time III: Total computation time of the state estimation program for the UAV model.

Figure 4 shows the overall diagram of the state estimation program for the UAV model and the locations of “Measurement Time I”, “Measurement Time II”, and “Measurement Time III”.

The controller design for generating generalized forces follows the computed torque method described in Section 5. The target trajectory involves a circular path from 10 s to 70 s, with yaw angle always at 0 deg, and roll and pitch angles changing to specified angles each revolution. However, the control inputs are first transformed from the generalized forces (thrust and torque for the motion system) provided by the computed torque method into actuator commands (rotor speeds and tilt angles) via the control allocation law. These are then multiplied by the propeller thrust coefficient $[eqn]$ and torque coefficient $[eqn]$ to reproduce the thrust and torque (including reaction torque) from the actuators, which are applied as inputs to the plant.

Below, we detail the data collection methods for Measurement Times I, II, and III.

Measurement Time I: Since one simulation yields N measurement data points due to the total number of sampling points N, the average of all these is taken to calculate the average computation time per loop. However, due to the nature of the AB method, the average is taken over $[eqn]$ data points for 2nd order, $[eqn]$ for 3rd order, $[eqn]$ for 4th order, $[eqn]$ for 5th order, and $[eqn]$ for 6th order.Measurement Time II: Similar to Measurement Time I, the average computation time per loop is calculated, and the total computation time for N loops is also calculated.Measurement Time III: This is the total program computation time, so one simulation yields a single measurement data point.

Using these measurement methods, 50 simulations are repeatedly executed, and their average computation times are calculated for each case.

Remark

The purpose of this study’s computational efficiency evaluation is to compare algorithm-specific computational loads using the MATLAB environment. The validity of this approach is based on the following points:

• Elimination of Hardware Dependence: MATLAB simulations eliminate hardware-dependent factors such as CPU architecture, memory bandwidth, and cache size, enabling a pure evaluation of the relative computational load differences between the proposed numerical integration method (the AB method) and the conventional method (the RK method).
• Indicator of Computational Amount: In an onboard environment (flight control computer), the sampling frequency is high and computing resources are limited, making it important to reduce the number of floating-point operations and execution time to achieve the same accuracy. The reduction in computation time measured in MATLAB serves as a basic indicator of the effectiveness of processing speed improvement on an actual aircraft.

However, final verification requires implementation on an actual autopilot (e.g., Pixhawk) and benchmarking, which is considered a future challenge. This simulation demonstrates the theoretical and algorithmic advantages prior to actual implementation.

7.6. Results and Discussion of Estimation Accuracy

7.6.1. Comparison of Three Methods in Estimation Accuracy: Euler, RK4, and AB4

First, to show the estimates for the UAV model and to verify whether the UKF based on the multi-step AB method can perform estimation as stably as the single-step Euler and RK methods, Figure 5 shows the time evolution of state estimates by UKF based on the three methods for $[eqn]$ .

As seen from Figure 5, the UKF estimates based on the AB method converge stably and achieve high-precision estimation. Also, no significant deviation from the RK method estimates is observed.

Next, Figure 6 shows the time evolution of the RMSE values for the state estimates by UKF based on the three methods for $[eqn]$ , over the 70 s simulation period. From Figure 6, under the condition $[eqn]$ , although the RMSE differences among the methods are very small for all state variables, relatively larger differences are observed in the RMSE values for positions $[eqn]$ , attitude angles $[eqn]$ , and attitude angular velocities $[eqn]$ . Also, for all 12 state variables, the time variation of RMSE values becomes smaller after around 50 s of simulation time.

In addition, estimation accuracy is compared under different discretization periods: $[eqn]$ . Figure 7 shows the change in $[eqn]$ with varying discretization period h. Figure 7 shows that as the discretization period h increases, the differences among the methods become larger, and ultimately the $[eqn]$ value for the RK method is the smallest. Furthermore, for the AB4 method, estimation values could not be obtained (diverged) for h values greater than or equal 0.08. Details regarding this are discussed later.

The specific values for each point in Figure 7 are shown in Table 4. The “N/A” in Table 4 indicates that estimation values could not be obtained due to divergence.

7.6.2. RMSE Comparison of AB Methods with Different Orders of Accuracy

The estimation values are omitted, as their differences from Figure 6 were very small. Figure 8 shows the time evolution of the RMSE values for UKF estimates based on the second- to sixth-order AB methods for $[eqn]$ . From Figure 8, under the condition $[eqn]$ , all AB methods from the second to sixth order converge to values comparable to the RMSE values for UKF estimates based on the Euler and RK methods shown in Figure 6, without significant deviation.

Next, estimation accuracy is compared under different discretization periods: $[eqn]$ . Figure 9 shows the change in $[eqn]$ with varying discretization period h. Similar to Figure 7, Figure 9 shows that as the discretization period h increases, the differences among the methods become larger. Also, the sixth-order AB method diverged for h greater than or equal 0.03; the fifth-order, for h greater than or equal 0.05; and the fourth-order, for h greater than or equal 0.08, preventing estimation value acquisition.

The specific values for each point in Figure 9 are shown in Table 5. Additionally, for the 2nd and 3rd-order AB methods, further verification confirmed that they similarly diverge for values greater than or equal $[eqn]$ and $[eqn]$ , respectively. Based on these results and additional experiments, the maximum step sizes at which each order can perform stable estimation without divergence are summarized as follows:

• 2nd-order AB method: Stable up to $[eqn]$ (verified by additional experiments up to $[eqn]$ ). Divergence occurs at $[eqn]$ .
• 3rd-order AB method: Stable up to $[eqn]$ . Divergence occurs at $[eqn]$ .
• 4th-order AB method: Stable up to $[eqn]$ . Divergence occurs at $[eqn]$ .
• 5th-order AB method: Stable up to $[eqn]$ . Divergence occurs at $[eqn]$ .
• 6th-order AB method: Stable up to $[eqn]$ . Divergence occurs at $[eqn]$ .

These experimentally obtained stability limits are consistent with the theoretical property of linear multi-step methods that the absolute stability region shrinks as the order increases [36,37].

In practical applications, the choice of step size involves a trade-off between estimation accuracy and computational efficiency. A larger step size reduces computational load but may degrade accuracy or cause divergence, while a smaller step size improves stability and accuracy at the cost of increased computation time. Users should select the step size based on the specific requirements of their application, referring to the experimentally verified stability limits above as a guideline for the maximum allowable step size.

7.6.3. Comparative Analysis: RK Method vs. AB Method

Comparing the results in Table 4 (re-entry vehicle problem) and Table 5 (UAV problem) reveals the following:

• Accuracy Difference Between the RK and AB Methods: At the same step size (e.g., $[eqn]$ ), the RMSE of the RK4-UKF (UAV: 0.47646) is nearly identical to the RMSE of the AB4-UKF (UAV: 0.47304), indicating no significant difference in accuracy between the two. This indicates that the RK method achieves high-order accuracy through multiple evaluations within a single step, while the AB method achieves similar accuracy with fewer function evaluations by utilizing past information.
• Relationship between Order and Stability: For both problems, the higher the order of the AB algorithm, the smaller the upper limit of the step size for stable estimation. This is due to the theoretical property of numerical analysis that the absolute stability region shrinks as the order of the linear multi-step algorithm increases.
• Problem Complexity and Computational Efficiency: For the simple 3D reentry vehicle problem (Table 2), the computational time of the AB4-based UKF increased by approximately 2.9% compared to the Euler-based UKF, while for the complex 12D UAV model (Table 6), it was reduced by approximately 5.1%. This reversal phenomenon is thought to be due to the characteristics of matrix operations in high-dimensional models (regularity of memory access patterns, cache efficiency).

These results indicate a trade-off between computational efficiency and numerical stability. When applying this method to real problems, it is necessary to select an appropriate order based on the required estimation accuracy and acceptable step size.

7.7. Results and Discussion of Computational Efficiency

7.7.1. Comparison of Three Methods in Computational Efficiency: Euler, RK4, and AB4

Table 6 shows the computation times for Time I, Time II, and Time III for the three methods (Euler, RK4, and AB4) with step size $[eqn]$ . From Table 6, the reduction rate in computation time for the 4th-order AB method is 12.3% compared to Euler and 41.1% compared to RK for Time I; 7.1% compared to Euler and 29.9% compared to RK4 for Time II; and 5.1% compared to Euler and 20.9% compared to RK4 for Time III.

Next, Figure 10 shows the change in computation time with varying discretization period h. Figure 10 shows that as h becomes smaller, the difference between the two methods increases. This is considered to be due to the influence of the increasing number of program loops as the total number of sampling points N becomes larger for smaller h.

7.7.2. Computational Time Comparison of AB Methods with Different Orders of Accuracy

Next, computation times for Time I, Time II, and Time III are compared under the same conditions for AB methods from 2nd to 6th order for $[eqn]$ , not just the AB4 method. The result is given in Table 7. The computational difference among AB methods from 2nd to 6th order lies only in the number of past sigma point matrices used as “previous information”. Therefore, the differences among methods are smaller than the results shown in Table 6.

Next, Figure 11 shows the change in computation time with varying discretization period h. From Figure 11, unlike the comparison between RK and AB methods, almost no difference is seen among the methods. The reason for this result is considered to be that the difference in computational load due to the order of the AB method stems only from the difference in the number of sigma point matrices handled, and the number of derivative calculations per step is the same for any order. Also, from this result, it can be understood that while estimation with the AB method diverges if the step size is too large, it does not diverge under conditions with h values smaller than 0.01, allowing estimation with higher-order accuracy and shorter computation time compared to the RK method under those conditions.

8. Conclusions

This paper proposed a UKF that newly integrates the Adams–Bashforth method into the sigma point time-update equation of the UKF (AB-UKF) as an alternative to the Runge–Kutta method. The proposed method maintains estimation accuracy comparable to the RK-UKF proposed by Takeno and Katayama [20] while enabling estimation with more efficient computation time.

The main findings obtained from the numerical simulation results are summarized below.

Maintenance of estimation accuracy: For both models, AB-UKF achieved estimation accuracy comparable to RK-UKF. For the UAV model ( $[eqn]$ ), no significant difference was observed between the RMSE of AB4-UKF (0.47304) and that of RK4-UKF (0.47646).Comparison of computational efficiency: Comparison of total computation time relative to the Euler-based UKF confirmed the effectiveness of the proposed method.

• –For the simple three-dimensional reentry vehicle problem (Table 2), the computation time of AB4-based UKF was approximately 2.9% longer than Euler-based UKF (slower than Euler), whereas RK4-based UKF was approximately 9.4% longer.
• –For the complex 12-dimensional UAV model (Table 6), the computation time of AB4-based UKF was approximately 5.1% shorter than Euler-based UKF (faster than Euler). In contrast, RK4-based UKF exhibited a more pronounced slowdown, with an approximately 20.0% increase. These results demonstrate that the computational time reduction effect achieved by the proposed AB method integration becomes more pronounced for models with more complex state equations and higher dimensions.Relationship between order and stability: Verification using the UAV model quantitatively showed that higher-order AB methods have a smaller upper limit on the step size that allows stable estimation. The maximum stable step sizes for each order are: AB2: $[eqn]$ , AB3: $[eqn]$ , AB4: $[eqn]$ , AB5: $[eqn]$ , and AB6: $[eqn]$ . Divergence occurs at step sizes exceeding these limits.

These results indicate that high-precision numerical integration methods are not necessarily superior in terms of computational efficiency; the RK4-based UKF was consistently slower than the Euler-based UKF in both problems (reentry problem: +9.4%, UAV problem: +20.0%). The proposed AB-UKF is an effective option for improving computational efficiency while maintaining estimation accuracy, particularly for high-dimensional models.

8.1. Limitations of This Study

This study has the following limitations:

• Linearity of observations: The observations in the UAV model are linear (direct observation of position and attitude angles), and the advantages of nonlinear filtering are not fully utilized. In actual UAV navigation, nonlinear observations such as radar observations in GPS-denied environments are common.
• Simulation environment: The evaluation of computational efficiency is based on MATLAB simulations and has not been verified for real-time operation on actual hardware (flight control computers).
• Application to specific models: The effectiveness of this method has been verified on two models (the falling body model and the UAV model), but additional verification is required for generalization to other nonlinear systems.

8.2. Future Work

Based on the above limitations, future work includes the following:

• Extension to nonlinear observations: Verify the effectiveness of the proposed method under nonlinear observations using only range and azimuth. Specifically, we will tackle the problem of estimating a 12-dimensional state vector using only range and azimuth information from one or two radars, and verify whether AB-UKF can maintain its advantage in terms of computational efficiency even under nonlinear observations.
• Reduction of the number of sigma points: While the standard UKF uses $[eqn]$ sigma points, the number can be reduced to $[eqn]$ or $[eqn]$ by using the Spherical Simplex Unscented Transformation [15,43,44,45] or the Simplex Unscented Transformation [46]. Combining these methods with the AB method is expected to further improve computational efficiency.
• Implementation on actual hardware: Implement the proposed algorithm on a microcontroller (e.g., Pixhawk) for the UAV model used in this simulation, and verify the degree of improvement in estimation accuracy and computational efficiency in actual flight environments.
• Adaptive step size control: Integrate the adaptive step size control based on the degree of nonlinearity proposed by Wang et al. [25] into AB-UKF to further optimize estimation accuracy and computational efficiency.

The proposed method provides a new option for the trade-off between computational efficiency and estimation accuracy in state estimation for high-dimensional nonlinear systems, and its future development is expected to lead to applications in practical UAV navigation systems.