A Symbolic Regression Method for Dynamic Modeling
and Control of Quadrotor UAVs

A QUAV model represents a nonlinear system of four inputs and six degrees of freedom. It is composed of four individual electric motors and an X-shape airframe, as shown in Fig. 1. The torque and thrust can be changed by adjusting the rotation speed of each propeller, thereby further accomplishing different motions.

The dynamics of a quadrotor are explained using two reference frames. The inertial frame $X_{e}$-$Y_{e}$-$Z_{e}$ is defined relative to the ground, with the gravity pointing along the positive direction of the $Z_{e}$ inertial axis. Let the vector $[x,y,z]^{\mbox{\tiny\sf T}}\in\mathbb{R}^{3}$ denote the position of the center of mass in the inertial frame. The velocity and Euler angles of the quadrotor with respect to the inertial frame are represented by the vector $[\dot{x},\dot{y},\dot{z}]^{\mbox{\tiny\sf T}}\in\mathbb{R}^{3}$ and $[\phi,\theta,\psi]^{\mbox{\tiny\sf T}}\in\mathbb{R}^{3}$, respectively. The body frame $X_{b}$-$Y_{b}$-$Z_{b}$ is associated with the orientation of the quadrotor, with the thrust direction always consistent with the negative direction of the $Z_{b}$ body axis. And the vector $[w_{x},w_{y},w_{z}]^{\mbox{\tiny\sf T}}\in\mathbb{R}^{3}$ represents the angular velocity in the body frame. It is worthwhile pointing out that the roll and pitch angles are limited to ($-\pi/2,\pi/2$) and the yaw angle is limited to ($-\pi,\pi$), which is physically meaningful. Then, the complete 12 dimensional state vector is denoted as

$\displaystyle\xi=$

$\displaystyle\left[\begin{array}[]{cccccccccccc}x&y&z&\dot{x}&\dot{y}&\dot{z}&\phi&\theta&\psi&w_{x}&w_{y}&w_{z}\end{array}\right]^{\mbox{\tiny\sf T}}.$

Let $\xi_{i}$, $i=1,\cdots,12$, be the $i$-th element of the vector $\xi$. The control input $u_{i}$ represents the rotation speed of each motor, for $i=1,\cdots,4$. Denote $u=\left[\begin{array}[]{cccc}u_{1}&u_{2}&u_{3}&u_{4}\end{array}\right]^{\mbox{\tiny\sf T}}$. The dynamic model is of the following form

	$\displaystyle\dot{\xi}_{1}$	$\displaystyle=\xi_{4}$
	$\displaystyle\dot{\xi}_{2}$	$\displaystyle=\xi_{5}$
	$\displaystyle\dot{\xi}_{3}$	$\displaystyle=\xi_{6}$
	$\displaystyle\dot{\xi}_{i}$	$\displaystyle=\zeta_{i}(\xi,u),\;i=4,\cdots,12$		(1)

with the some functions $\zeta_{i}(\xi,u),\;i=4,\cdots,12$, to be determined.

Symbolic regression (SR) is an established method that can effectively generate theories of causation in the form of mathematical equations from a collected data set. Unlike the traditional regression methods that require a fixed-form model derived from priori knowledge, SR automatically searches both the parameters and the form of equations. The expression of equations is composed of common functions ($\sin$, $\cos$, $\sqrt{\cdot}$), algebraic operators ($+$, $-$, $\times$, $\div$), constants and relative variables. The selection of equations depends on its complexity and fitness. And the fitness function $e(f^{*}(x),y)=\sum_{i=1}^{n}\lvert f^{*}(x^{i})-y^{i}\rvert$ is widely used in SR, where $x$ denotes the input variable and $y$ the output variable, $n$ is the amount of pairs of input-output data $(x^{i},y^{i})$, and $f^{*}$ is the equation learned by SR. Since the computation complexity caused by a large search space, genetic programming (GP) [23] is typically implemented in SR. GP is a meta-heuristic algorithm that adopts the principle of biological evolution to update candidate solutions.

In the initial phase of modeling, we generate multiple groups of random propeller rotation speeds within different ranges as the input data set. Then, by applying the input data set on the UAV simulator, we collect the output data set which includes velocity, acceleration, angular velocity, and angular acceleration. To perform symbolic regression modeling, we import the input and output data into the Eureqa software, select the formula building-blocks, e.g., $\left\{{\rm constant},+,-,\times,\div,\sin,\cos\right\}$ and the absolute error metric for training. As a result, the analytic equations through efficient training and sifting are obtained as follows,

	$\displaystyle\zeta_{4}(\xi,u)=$	$\displaystyle-0.00768u_{1}\cdot\sin(\xi_{7})\cdot\sin(\xi_{9})-1.63e^{-5}u_{2}\cdot u_{4}$
		$\displaystyle\cdot\sin(\xi_{7})\cdot\sin(\xi_{9})-0.00665\xi_{8}\cdot u_{1}\cdot\cos(\xi_{7})$
		$\displaystyle\cdot\cos(\xi_{9})-1.63e^{-5}\xi_{8}\cdot u_{2}\cdot u_{4}\cdot\cos(\xi_{7})\cdot\cos(\xi_{9})$
	$\displaystyle\zeta_{5}(\xi,u)=$	$\displaystyle 7.87e^{-6}u_{1}^{2}\cdot\sin(\xi_{7}+\xi_{9})+7.78e^{-6}u_{2}^{2}$
		$\displaystyle\cdot\sin(\xi_{7}+\xi_{9})+7.78e^{-6}u_{3}^{2}\cdot\sin(\xi_{7}+\xi_{9})$
		$\displaystyle+7.78e^{-6}u_{4}^{2}\cdot\sin(\xi_{7}+\xi_{9})-0.0441$
	$\displaystyle\zeta_{6}(\xi,u)=$	$\displaystyle 9.44-6.90e^{-6}\cos(\xi_{7})\cdot(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+u_{4}^{2})$
		$\displaystyle-0.000289\xi_{6}^{2}$
	$\displaystyle\zeta_{7}(\xi,u)=$	$\displaystyle\sin(\xi_{7})\cdot\frac{\sin(\xi_{8})}{\cos(\xi_{8})}\cdot\xi_{11}+\cos(\xi_{7})\cdot\frac{\sin(\xi_{8})}{\cos(\xi_{8})}\cdot\xi_{12}$
		$\displaystyle+\xi_{10}$
	$\displaystyle\zeta_{8}(\xi,u)=$	$\displaystyle\cos(\xi_{7})\cdot\xi_{11}-\sin(\xi_{7})\cdot\xi_{12}$
	$\displaystyle\zeta_{9}(\xi,u)=$	$\displaystyle\frac{\sin(\xi_{7})}{\cos(\xi_{8})}\cdot\xi_{11}+\frac{\cos(\xi_{7})}{\cos(\xi_{8})}\cdot\xi_{12}$
	$\displaystyle\zeta_{10}(\xi,u)=$	$\displaystyle 0.697\xi_{11}\cdot\xi_{12}+8.33e^{-5}(u_{2}^{2}+u_{3}^{2}-u_{1}^{2}-u_{4}^{2})$
	$\displaystyle\zeta_{11}(\xi,u)=$	$\displaystyle-0.708\xi_{10}\cdot\xi_{12}+8.03e^{-5}(u_{1}^{2}+u_{3}^{2}-u_{2}^{2}-u_{4}^{2})$
	$\displaystyle\zeta_{12}(\xi,u)=$	$\displaystyle 0.0219\xi_{10}\cdot\xi_{11}+4.86e^{-6}(u_{1}^{2}+u_{2}^{2}-u_{3}^{2}-u_{4}^{2}).$

In particular, the functions $\zeta_{i}(\xi,u)$, $i=7,8,9$, represent the conversion relationship between $\{\dot{\xi}_{7},\dot{\xi}_{8},\dot{\xi}_{9}\}$ and $\{\xi_{10},\xi_{11},\xi_{12}\}$, which are considered as known transform equations. The remaining functions are obtained from different training data via symbolic regression. For instance, in order to acquire the function $\zeta_{4}(\xi,u)$, we generate four groups of training data (2,000 samples) using the inputs in different ranges of 0-300rad/s, 0-700rad/s, 0-800rad/s and 0-1000rad/s. Both model validation and controller design are implemented on the equations learned from different training data, and the equation that performs the best is selected in the final model. The other functions in (1) are obtained in a similar way. It is worth mentioning that the propeller speeds for $u_{i},i=1,\cdots,4,$ are constrained between 0 and 1000rad/s.

To validate the learned dynamic model, a test data set is generated using the following excitation inputs for 10 seconds,

	$\displaystyle u_{1}=$	$\displaystyle 300+300\sin(10t)$
	$\displaystyle u_{2}=$	$\displaystyle 300+300\sin(10t+\pi/4)$
	$\displaystyle u_{3}=$	$\displaystyle 300+300\sin(10t+\pi/3)$
	$\displaystyle u_{4}=$	$\displaystyle 300+300\sin(10t+\pi/6).$

The reliability of the model is evaluated by comparing the six target variables ($\ddot{x},\ddot{y},\ddot{z},\dot{w}_{x},\dot{w}_{y},\dot{w}_{z})$ obtained from the actual system and the SR model, respectively, and the result is shown in Fig. 2.

Generally speaking, symbolic regression can discover the accurate model equations beneath data, especially for angular acceleration $\dot{w_{x}},\dot{w_{y}},\dot{w_{z}}$, which indicates that symbolic regression is an effective modeling approach. Moreover, the RMSE, MAE and $R^{2}$ indexes are applied to quantitively reflect the approximation ability of SR; the results are summarized in Table I. The accuracy of the analytic model is negatively correlated with the values of RMSE, MAE and $|1-R^{2}|$.

The objective of this part is to design a proportional-integral (PI) controller for the outer loop subsystem to ensure the roll and pitch angles to be automatically tuned to drive the quadrotor to the desired horizontal position $x_{d},y_{d}$. Based on (1), the horizontal position subsystem is described as follows:

	$\displaystyle\dot{\xi}_{1}$	$\displaystyle=\xi_{4}$
	$\displaystyle\dot{\xi}_{2}$	$\displaystyle=\xi_{5}$
	$\displaystyle\dot{\xi}_{4}$	$\displaystyle=\zeta_{4}(\xi_{7},\xi_{8},\xi_{9},u)$
	$\displaystyle\dot{\xi}_{5}$	$\displaystyle=\zeta_{5}(\xi_{7},\xi_{9},u)$		(4)

where $\xi_{7}=\phi$ and $\xi_{8}=\theta$ are regarded as the virtual controls to this subsystem.

More specifically, we aim to solve the desired angle trajectories for $\xi_{7}$ and $\xi_{8}$, i.e., $\xi_{7d}=\phi_{d}$ and $\xi_{8d}=\theta_{d}$, from the following equations

	$\displaystyle\zeta_{4}(\xi_{7d},\xi_{8d},\xi_{9},u)=k_{px}\dot{e}_{x}+k_{ix}e_{x}+\ddot{x}_{d}$
	$\displaystyle\zeta_{5}(\xi_{7d},\xi_{9},u)=k_{py}\dot{e}_{y}+k_{iy}e_{y}+\ddot{y}_{d}.$		(5)

The solutions to (6) can be explicitly expressed as follows

	$\displaystyle\xi_{7d}=$	$\displaystyle\sin^{-1}((k_{py}\dot{e}_{y}+k_{iy}e_{y}+\ddot{y}_{d}+0.0441)$
		$\displaystyle/(7.87e^{-6}u_{1}^{2}+7.78e^{-6}u_{2}^{2}+7.78e^{-6}u_{3}^{2}$
		$\displaystyle+7.78e^{-6}u_{4}^{2}))-\xi_{9}$
	$\displaystyle\xi_{8d}=$	$\displaystyle(k_{px}\dot{e}_{x}+k_{ix}e_{x}+\ddot{x}_{d}+0.00768u_{1}$
		$\displaystyle\cdot\sin(\xi_{7d})\cdot\sin(\xi_{9})+1.63e^{-5}u_{2}\cdot u_{4}$
		$\displaystyle\cdot\sin(\xi_{7d})\cdot\sin(\xi_{9}))/(-0.00665u_{1}\cdot\cos(\xi_{7d})$
		$\displaystyle\cdot\cos(\xi_{9})-1.63e^{-5}u_{2}\cdot u_{4}\cdot\cos(\xi_{7d})\cdot\cos(\xi_{9})).$		(6)

When $\xi_{7}$ and $\xi_{8}$ achieve $\xi_{7d}$ and $\xi_{8d}$, respectively, (to be elaborated in the next subsection), the closed-loop system (4) with $\xi_{7}=\xi_{7d}$ and $\xi_{8}=\xi_{8d}$ becomes

	$\displaystyle 0$	$\displaystyle=\ddot{e}_{x}+k_{px}\dot{e}_{x}+k_{ix}e_{x}$
	$\displaystyle 0$	$\displaystyle=\ddot{e}_{y}+k_{py}\dot{e}_{y}+k_{iy}e_{y}$		(7)

With proper selection of the parameters $k_{px}$, $k_{ix}$ $k_{py}$, and $k_{iy}$, the control objective (2) can be achieved.

The design of $\xi_{7d}$ and $\xi_{8d}$, as a virtual controller to the system (4), is of a PI structure containing the proportional terms $k_{px}\dot{e}_{x}$ and $k_{py}\dot{e}_{y}$ and the integral terms $k_{ix}e_{x}=k_{ix}\int\dot{e}_{x}$ and $k_{iy}e_{y}=k_{iy}\int\dot{e}_{y}$, as well as feedforward compensation.

In this part, a backstepping controller is designed to ensure $\xi_{7}=\phi$ and $\xi_{8}=\theta$ achieve $\xi_{7d}=\phi_{d}$ and $\xi_{8d}=\theta_{d}$, respectively, for the objective (2) as explained in the previous subsection, that is,

$\displaystyle\lim_{t\rightarrow\infty}e_{\phi}(t)=0,\;\lim_{t\rightarrow\infty}e_{\theta}(t)=0$

(8)

for

$\displaystyle e_{\phi}=\phi_{d}-\phi,\;e_{\theta}=\theta_{d}-\theta.$

Also, it is to ensure the objective (3). In other words, this subsection is concerned about the subsystem governing

		$\displaystyle\left[\begin{array}[]{cccccccc}\xi_{3}&\xi_{6}&\xi_{7}&\xi_{8}&\xi_{9}&\xi_{10}&\xi_{11}&\xi_{12}\end{array}\right]^{\mbox{\tiny\sf T}}$		(10)
	$\displaystyle=$	$\displaystyle\left[\begin{array}[]{cccccccc}z&\dot{z}&\phi&\theta&\psi&w_{x}&w_{y}&w_{z}\end{array}\right]^{\mbox{\tiny\sf T}}$		(12)

and the control objectives are (3) and (8) where $z_{d},\psi_{d}$ are the desired trajectories and $\phi_{d}$ and $\theta_{d}$ designed in (6).

We introduce the following coordinate transformation

$$\begin{bmatrix}\dot{\phi}\\ \dot{\theta}\\ \dot{\psi}\\ \end{bmatrix}=R\begin{bmatrix}w_{x}\\ w_{y}\\ w_{z}\\ \end{bmatrix}$$

(13)

with

$$R=\begin{bmatrix}1&\frac{\sin(\phi)\sin(\theta)}{\cos(\theta)}&\frac{\cos(\phi)\sin(\theta)}{\cos(\theta)}\\ 0&\cos(\phi)&-\sin(\phi)\\ 0&\frac{\sin(\phi)}{\cos(\theta)}&\frac{\cos(\phi)}{\cos(\theta)}\end{bmatrix},$$

then the variables $[\ddot{\phi},\ddot{\theta},\ddot{\psi}]^{T}$ can be obtained as follows:

$$\begin{bmatrix}\ddot{\phi}\\ \ddot{\theta}\\ \ddot{\psi}\\ \end{bmatrix}=R\begin{bmatrix}\dot{w}_{x}\\ \dot{w}_{y}\\ \dot{w}_{z}\\ \end{bmatrix}+\dot{R}\begin{bmatrix}w_{x}\\ w_{y}\\ w_{z}\\ \end{bmatrix}.$$

(14)

Now, the subsystem governing the states (12) can be equivalently transformed to a systems governing the following new state vector

$\displaystyle s=$

$\displaystyle\left[\begin{array}[]{cccccccc}z&\phi&\theta&\psi&\dot{z}&\dot{\phi}&\dot{\theta}&\dot{\psi}\end{array}\right]^{\mbox{\tiny\sf T}}.$

(16)

Let $s_{i}$ be the $i$-th element of $s$, for $i=1,\cdots,8$. Now, the dynamics for $s$ can be reorganized as follows

	$\displaystyle\dot{s}_{1}$	$\displaystyle=s_{5}$
	$\displaystyle\dot{s}_{2}$	$\displaystyle=s_{6}$
	$\displaystyle\dot{s}_{3}$	$\displaystyle=s_{7}$
	$\displaystyle\dot{s}_{4}$	$\displaystyle=s_{8}$
	$\displaystyle\dot{s}_{i}$	$\displaystyle=\gamma_{i}(s,f),\;i=5,\cdots,8$		(17)

where

$\displaystyle\gamma_{5}(s,f)=$

$\displaystyle 9.44-6.9e^{-6}\cos(s_{2})\cdot f_{1}-0.000289s_{5}^{2}$

and $\gamma_{i}(s,f),\;i=6,7,8$ can be derived in a similar manner (with the details omitted). The vector $f=\left[\begin{array}[]{cccc}f_{1}&f_{2}&f_{3}&f_{4}\end{array}\right]^{\mbox{\tiny\sf T}}$ represents the new control variables defined as follows

	$\displaystyle f_{1}=u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+u_{4}^{2}$
	$\displaystyle f_{2}=u_{2}^{2}+u_{3}^{2}-u_{1}^{2}-u_{4}^{2}$
	$\displaystyle f_{3}=u_{1}^{2}+u_{3}^{2}-u_{2}^{2}-u_{4}^{2}$
	$\displaystyle f_{4}=u_{1}^{2}+u_{2}^{2}-u_{3}^{2}-u_{4}^{2}.$		(18)

Let

	$\displaystyle s_{1d}=z_{d},\;s_{2d}=\phi_{d},\;s_{3d}=\theta_{d},\;s_{4d}=\psi_{d}$
	$\displaystyle e_{i}=s_{id}-s_{i},\;i=1,\cdots,4.$

Now, the objectives (3) and (8) are converted to design of $f_{i}\;i=1,\cdots,4$ for (17) such that

$\displaystyle\lim_{t\rightarrow\infty}e_{i}(t)=0,\;i=1,\cdots,4.$

(19)

The design follows a backstepping procedure, see, e.g., [24]. First, we consider the Lyapunov function candidates

$$V_{i}(e_{i})=\frac{1}{2}e_{i}^{2},\;i=1,\cdots,4.$$

(20)

The derivative of $V_{i}(e_{i})$ along the trajectory of (17) satisfies

$$\dot{V}_{i}(e_{i})=e_{i}\dot{e}_{i}=e_{i}(\dot{s}_{id}-s_{i+4}).$$

(21)

Here, $s_{i+4}$ is the virtual control whose desired value is

$$s_{(i+4)d}=\dot{s}_{id}+c_{i}e_{i}.$$

(22)

Let

$$\delta_{i}=s_{(i+4)d}-s_{i+4}.$$

(23)

Now, we consider the Lyapunov function candidates involving both $e_{i}$ and $\delta_{d}$ as follows

$$W_{i}(e_{i},\delta_{i})=\frac{1}{2}e_{i}^{2}+\frac{1}{2}\delta_{i}^{2},\;i=1,\cdots,4,$$

(24)

whose time derivative satisfies

		$\displaystyle\dot{W}_{i}(e_{i},\delta_{i})$
	$\displaystyle=$	$\displaystyle e_{i}\dot{e}_{i}+\delta_{i}\dot{\delta}_{i}$
	$\displaystyle=$	$\displaystyle e_{i}(\dot{s}_{id}-s_{i+4})+\delta_{i}(\dot{s}_{(i+4)d}-\dot{s}_{i+4})$
	$\displaystyle=$	$\displaystyle e_{i}(s_{(i+4)d}-c_{i}e_{i}-s_{i+4})+\delta_{i}(\dot{s}_{(i+4)d}-\gamma_{i+4}(s,f))$
	$\displaystyle=$	$\displaystyle e_{i}(\delta_{i}-c_{i}e_{i})+\delta_{i}(\dot{s}_{(i+4)d}-\gamma_{i+4}(s,f))$
	$\displaystyle=$	$\displaystyle-c_{i}e_{i}^{2}+\delta_{i}(e_{i}+\dot{s}_{(i+4)d}-\gamma_{i+4}(s,f)).$

Letting

	$\displaystyle\gamma_{i+4}(s,f)=$	$\displaystyle e_{i}+\dot{s}_{(i+4)d}+c_{i+4}\delta_{i}$
	$\displaystyle=$	$\displaystyle s_{id}-s_{i}+\ddot{s}_{id}+c_{i}(\dot{s}_{id}-s_{i+4})$
		$\displaystyle+c_{i+4}(\dot{s}_{id}+c_{i}(s_{id}-s_{i})-s_{i+4})$		(25)

gives

$\displaystyle\dot{W}_{i}(e_{i},\delta_{i})=-c_{i}e_{i}^{2}-c_{i+4}\delta_{i}^{2}.$

If the two parameters $c_{i}$ and $c_{i+4}$ are selected as positive constants, one has the target (19) achieved.

The control input $f$ can be solved from the equation (25) for $i=1,\cdots,4$. For example, the input $f_{1}$ which stabilizes the altitude subsystem is thus presented as follows

	$\displaystyle f_{1}=$	$\displaystyle(s_{1d}-s_{1}+\ddot{s}_{1d}+c_{1}(\dot{s}_{1d}-s_{5})+c_{5}(\dot{s}_{1d}$
		$\displaystyle+c_{1}(s_{1d}-s_{1})-s_{5})-9.44+0.000289s_{5}^{2})$
		$\displaystyle/(-6.9e^{-6}\cos(s_{2})).$

The inputs $f_{2},f_{3},f_{4}$ can be solved in a similar manner. And the control laws $u_{1},u_{2},u_{3},u_{4}$ can be solved from (18).