\section{INTRODUCTION}
Accurate trajectory tracking for autonomous Unmanned Aerial Vehicles (UAVs), such as quadrotors, is necessary for maintaining autonomy. Although industrial applications of autonomous UAVs, such as commercial deliveries, search-and-rescue \cite{mishra2020drone} and wireless power transfer \cite{qian2022practical}, have attracted much attention in recent years, precisely tracking high-speed and high-acceleration UAV trajectories is an extremely challenging control problem, particularly in unknown and dynamic environments with unpredictable aerodynamic forces.

To achieve safe, precise and reliable quadrotor trajectory tracking, there are two main problems that need to be solved: How can we achieve robust and feasible estimation (or modelling) of the aerodynamic effects on quadrotors in complex dynamic environments? And; How can the whole control framework be integrated with aerodynamic effect estimation to solve the uncertainties and disturbances while tracking trajectory references precisely and reliably?

Previous work has shown that the primary source of uncertainties are aerodynamic effects deriving from drag forces and moment variations caused by the rotors and the fuselage \cite{torrente2021data}. Prominent aerodynamic effects appear at flight speeds of 5 $ms^{-1}$ in wind tunnel experiments by \cite{faessler2017differential}. These effects acting on quadrotors are chaotic and hard to model directly, as they are generated from a combination of the individual propellers and airframe \cite{hoffmann2007quadrotor}, turbulent effects caused by rotor–rotor and airframe–rotor interactions \cite{russell2016wind}, and the propagation of other turbulence \cite{kaya2014aerodynamic}.

Most current approaches to quadrotor trajectory tracking treat aerodynamic effects as simple external disturbances, and do not account for higher-order effects or attempt to deviate from a determined plan \cite{tal2020accurate,ru2017nonlinear,bicego2020nonlinear}. While these solutions are efficient and feasible for lightweight on-board computers, aggressive maneuvers at high speed, e.g., greater than 5 $ms^{-1}$, introduce large positional and attitude tracking errors. Recent data-driven approaches, such as Gaussian Processes (GP) \cite{torrente2021data,wang2022kinojgm} and neural networks \cite{spielberg2021neural} combined with Model Predictive Control (MPC), show accurate modelling of aerodynamic effects. However, due to the nonparametric nature of GP, the GP-based approaches perform poorly in complex environments - where large datasets contain drastic changes in wind speed and heading. In these instances, learning-based (neural networks) approaches perform better than those GP-based approaches \cite{loquercio2021learning}. Achieving adaptability and robustness in complex environments is still challenging, however, primarily because training datasets are collected from simulated platforms and real-world historical records that do not fully describe the complex environments.

In comparison to existing data-driven approaches, Reinforcement Learning (RL), an interactive learning process, is able to learn complex and changeable disturbances - i.e., the errors between the true and estimated values - using much less model information \cite{zhang2021model}. The key challenge of most existing RL approaches \cite{lillicrap2015continuous} is that policy optimization biases toward actions with high variance value estimates, since some of these values will be overestimated by random chance \cite{ma2021conservative}. In risk-sensitive or safety-critical applications such as autonomous quadrotor navigation these actions should be avoided. Recent work on distributional RL \cite{bellemare2017distributional} was proposed to approximate and parameterize the entire distribution of future rewards, instead of the expected value. Distributional RL algorithms have been operated to achieve advanced results on continuous control domains \cite{hessel2018rainbow}. In princple, they provide more complete and richer value-distribution information to enable a more stable learning process \cite{bellemare2017distributional}. Previous distributional RL algorithms parameterize the policy value distribution in different ways, including canonical return atoms \cite{bellemare2017distributional}, the expectiles \cite{rowland2019statistics}, the moments \cite{nguyen2021distributional}, and the quantiles \cite{dabney2018distributional,zhang2019quota}. The quantile approach is especially suitable for autonomous UAV trajectory tracking due to its risk-sensitive policy optimization.

Robust MPC for tracking control of uncertain systems like quadrotors is rapidly developing thanks to advances in hardware and algorithmic efficiency \cite{nguyen2021model}. The robust control approach enables a `worst case' formulation to analyze the stability and performance of a system under bounded uncertainties and disturbances \cite{zhou1998essentials}. However, with consideration of the uncertainties in real-world scenarios, such worst-case design renders the optimal control actions inherently inadequate and overly conservative in practise \cite{mayne2016robust,mesbah2016stochastic}. To avoid the conservatism of the worst case design, Stochastic MPC (SMPC) \cite{schwarm1999chance} uses the probabilistic descriptions, such as stochastic constraints (also called chance constraints), to predict probability distributions of system states within acceptable levels of risk in the receding-horizon optimization \cite{primbs2009stochastic}.

The core challenges for SMPC include: 1) optimizing the feedback control laws over arbitrary nonlinear functions \cite{munoz2020convergence}; 2) the chance constraints are non-convex and intractable \cite{mesbah2016stochastic,zhang2021stochastic}; and 3) the computational complexity will grow dramatically as more uncertainties are added. To address the first challenge, one solution is to use affine parameterization of the control policy over finite horizons. However, this approach cannot guarantee convexity, i.e., the second challenge, where the policy set may still be convex \cite{lofberg2003approximations}. Another solution is an Affine Disturbance Feedback (ADF) control parameterization, proposed in \cite{goulart2008input}. This ADF control parameterization can address the first two challenges, which are optimizing the dynamic function and guaranteeing the decision variables to be convex, respectively. However, the main weakness is that the computational complexity grows quadratically with the prediction horizon, i.e., the third challenge. To overcome this difficulty, a Simplified Affine Disturbance Feedback (SADF) proposed in \cite{zhang2021stochastic}, where the SADF is equivalent to ADF but a finite-horizon optimization can be computed more efficiently using CasADi \cite{andersson2019casadi}, a nonlinear MPC solver. \cite{zhang2021stochastic} achieve good results by implementing SADF with zero-mean disturbance, however it is unclear how the SADF would perform on systems with non-zero-mean disturbances, such as a quadrotor.

To address the two stated issues, we propose \textbf{Qua}ntile-approximation-based \textbf{D}istributional \textbf{Re}inforced-\textbf{D}isturbance-estimation for Stochastic MPC (QuaDRED-SMPC), a systematic, safe and feasible quadrotor trajectory tracking framework for use with high variance aerodynamic effects. The details are as follows:
\begin{enumerate}
    \item Aerodynamic Disturbance Estimator: a Quantile-approximation-based Distributional Reinforced-disturbance-estimator (QuaDRED), described by Algorithm~\autoref{DRL_QuaDRED}, is proposed for aerodynamic disturbance estimation.
    QuaDRED builds upon prior QR-DQN \cite{dabney2018distributional} and QUOTA \cite{zhang2019quota} insofar as QuaDRED is a quantile-approximated distributional RL which uses a set of quantiles to approximate the full value distribution. In Section \uppercase\expandafter{\romannumeral4}-A, theoretical guarantees on convergence of the QuaDRED are provided based on policy evaluation (\textit{Proposition 3}) and policy improvement (\textit{Proposition 4}), respectively.

    \item Trajectory Tracker: Similar to \cite{zhang2021stochastic}, a Simplified Affine Disturbance Feedback (SADF) is used for control parameterization in SMPC (Algorithm~\autoref{SADF_SMPC}), where the convexity can be guaranteed in this process \cite{lofberg2003approximations} and computational complexity can be reduced. Different from prior work assuming zero mean disturbance, we consider the control performance and stability  under non-zero-mean disturbance. We use an Input-to-State Stability (ISS) \cite{jiang2001input} property to find conditions that imply stability and convergence of the tracker.
    \item The QuaDRED-SMPC framework is proposed to track quadrotor trajectory accurately under high variance aerodynamic effects. The overall control framework is shown in Fig.~\ref{RL_Control_framework}. In Section \uppercase\expandafter{\romannumeral4}-B, the closed-loop stability of the QuaDRED-SMPC framework is demonstrated under Lipschitz Lyapunov function \cite{rifford2000existence}.
\end{enumerate}

Our contributions can be summarized as follows:
\begin{itemize}
\item[1)]
QuaDRED, a distributional RL with quantile approximation that can sufficiently estimate variable aerodynamic disturbances. In the cases tested, we show that QuaDRED outperforms traditional RL, such as Deep Deterministic Policy Gradient (DDPG) \cite{lillicrap2015continuous}, and prior Distributional RL approaches, such as C51 \cite{bellemare2017distributional}.
\end{itemize}
\begin{itemize}
\item[2)]
The integration of a trajectory tracker with an aerodynamic disturbance estimator, a quadrotor trajectory tracking framework that integrates QuaDRED into a stochastic optimal control problem.
\end{itemize}
\begin{itemize}
\item[3)]
Convergence and stability guarantees: mathematical proofs are provided for the convergence of distributional-RL-based estimator, and the closed-loop stability of stochastic-MPC-based tracker with consideration of non-zero-mean and bounded disturbances.
\end{itemize}

\section{Problem Formulation}
A quadrotor dynamic model has six Degrees of Freedom (DoF), i.e., three linear motions and three angular motions \cite{torrente2021data}. We consider a nonlinear discrete system of quadrotor dynamics with state $\bm{x}\in\mathbb{X}\subseteq\mathbb{R}^{n}$, an additive disturbance $\bm{w}\in\mathbb{W}\subseteq\mathbb{R}^{n_w}$, and control input $\bm{u}\in\mathbb{R}^{n_u}$, defined for all time steps $k\in\mathbb{N}$ by:
\begin{equation}
\bm{x}_{k+1} = \bm{f}(\bm{x}_{k},\bm{u}_{k}, {\bm{e_{f}}}_k) + \bm{w}_k
\label{quadrotor_nonlinear_discrete_system}
\end{equation}
where $\bm{x}_{k}$ and $\bm{u}_{k}$ are the discrete-time state and input vectors of the quadrotor dynamic model. $\bm{w}_{k}$ is a disturbance caused by aerodynamic effects. $\bm{f}(\bm{x},\bm{u}, \bm{e}_{f})$ in \autoref{quadrotor_dynamics} is the continuous-time nominal model of the quadrotor integrating the aerodynamic effect $\bm{e}_{f}$. The state and input vectors of the nominal model are $\bm{x} = [\bm{P}_{WB},\bm{V}_{WB},\bm{q}_{WB},\bm{\omega} _{B}]^T$ and $\bm{u} = T_i, \forall i\in(0,3)$. $\bm {P}_{WB}$, $\bm{V}_{WB}$ and $\bm{q}_{WB}$ are the position, linear velocity and orientation of the quadrotor, and $\bm{\omega}_{B}$ is the angular velocity \cite{wang2022kinojgm}.
\begin{equation}
\begin{aligned}
&\dot{\bm{P}}_{WB} = \bm{V}_{WB}\\
&\dot{\bm{V}}_{WB} = \bm{g}_W + \frac{1}{m}(\bm{q}_{WB}\odot \bm{c} + \bm{e}_f)\\
&\dot{\bm{q}}_{WB} = \frac{1}{2}\Lambda(\bm{\omega} _{B})\bm{q}_{WB}\\
&\dot{\bm{\omega}} _{B} = \bm{J}^{-1}(\bm{\tau}_B-\bm{\omega}\times J\bm{\omega} _{B})
\end{aligned}
\label{quadrotor_dynamics}
\end{equation}
where $\bm{c}$ is the collective thrust $\bm{c} = [0,0,\sum T_i]^{\rm{T}}$ and $\bm{\tau}_B$ is the body torque; $\bm{g}_W = [0,0,-g]^{\rm{T}}$. The operator $\odot$ denotes a rotation of the vector by the quaternion. The skewsymmetric matrix $\Lambda(\bm{\omega})$ is defined in \cite{wang2022kinojgm}.

We linearize and reformulate \cite{goulart2008input} \autoref{quadrotor_nonlinear_discrete_system} for MPC over a finite horizon $N$:
\begin{equation}
\bm{x}_{t} = \bm{A}x_{0|t} + \bm{B}\bm{u}_{t} + \bm{G}\bm{w}_{t}
\label{linear_discrete_system}
\end{equation}
where $\bm{x}_{t}=[x_{0|t}^{\rm{T}}, x_{1|t}^{\rm{T}}, ..., x_{N|t}^{\rm{T}}]^{\rm{T}}$ and $\bm{u}_{t}=[u_{0|t}^{\rm{T}}, u_{1|t}^{\rm{T}}, ..., u_{N|t}^{\rm{T}}]^{\rm{T}}$ are the sequential states and inputs, and $\bm{w}_{t}=[w_{0|t}^{\rm{T}}, w_{1|t}^{\rm{T}}, ..., w_{N|t}^{\rm{T}}]^{\rm{T}}$ denotes a sequential stochastic disturbance over a horizon of $N$. $\bm{A}$, $\bm{B}$ and $\bm{G}$ are matrices defined in \cite{zhang2021stochastic}.

\begin{comment}
\begin{equation}
\begin{aligned}
&\bm{A} := \begin{bmatrix} I \\
A \\
A^2 \\
\vdots \\
A^N \end{bmatrix},
\bm{B} := \begin{bmatrix} 0 & 0 & \hdots & 0 \\
B & 0 & \hdots & 0 \\
AB & B & \hdots & 0 \\
\vdots & \ddots & \vdots & \vdots \\
A^{N-1}B & \hdots & AB & B \end{bmatrix} \\
&\bm{G} := \begin{bmatrix} 0 & 0 & \hdots & 0 \\
G & 0 & \hdots & 0 \\
AG & G & \hdots & 0 \\
\vdots & \ddots & \vdots & \vdots \\
A^{N-1}G & \hdots & AG & G \end{bmatrix}
\end{aligned}
\end{equation}
\end{comment}

The following assumptions are made:
\begin{comment}
\textit{Assumption 1}: The sequential disturbance $\bm{w}_{t}=[w_{0|t}^{\rm{T}}, w_{1|t}^{\rm{T}}, ..., w_{N|t}^{\rm{T}}]^{\rm{T}}$ is independent, and bounded -i.e., the probability density function (PDF) of $w_{n|t}, \forall n\in[0,N]$ is supported in a compact (closed and bounded) set $\mathbb{W}$ which contains the origin in its interior \cite{munoz2020convergence}.

Note that stochastic disturbance sequence $[w_{0|t}^{\rm{T}}, w_{1|t}^{\rm{T}}, ..., w_{N|t}^{\rm{T}}]^{\rm{T}}$ is only bounded but not a zero-mean distribution -i.e., $\mathbb{E}{(w_k)}\neq0$.
\end{comment}

\textit{Assumption 1}: Matrix G is column full rank.
\begin{comment}
\textit{Theorem 1(see \cite{jiang2001input})}: The dynamic system (\ref{quadrotor_nonlinear_discrete_system}) is Input-to-output Stable (ISS) if, and only if, there exist $\mathcal{H}_{\infty}$-functions $\alpha_1(\cdot)$, $\alpha_2(\cdot)$ and $\alpha_3(\cdot)$, a $\mathcal{H}$-function $\sigma (\cdot)$ and a continuous function $V$: $\mathbb{X}\rightarrow\mathbb{R_+}$ such that
\begin{equation}
\begin{aligned}
&\alpha_1(\left \| x \right\|) \leq V(x) \leq \alpha_2(\left \| x \right\|), \forall x\in\mathbb{R}^n\\
&V(f(x)) - V(x) \leq -\alpha_3(\left \| x \right\|) + \sigma (\left \| x \right\|)\\
\end{aligned}
\label{ISS}
\end{equation}
\end{comment}

\textit{Assumption 2}: The aerodynamic effect $\bm{e_{fk}}$ is available with no delay at each sampling timestamp.

\textit{Proposition 1}: There exists a control law $\bm{u}_b$ that ensures the nominal model $\bm{f}(\bm{x}_{k},\bm{u}_{k}, \bm{e_{fk}})$ is ISS if the stochastic disturbance $\bm{w}_t$ is independent and identically distributed (i.i.d.) a zero-mean distribution, i.e., $\mathbb{E}{(w_k)}=0$.

\textit{Proof}: Based on \textit{Assumption 2}, the nominal model $\bm{f}(\bm{x}_{k},\bm{u}_{k}, \bm{e_{fk}})$ is seen as $\bm{f}(\bm{x}_{k},\bm{u}_{k}^{'})$, in which the aerodynamic force $\bm{e_{fk}}$ is a constant term. Then we prove $\bm{f}(\bm{x}_{k},\bm{u}_{k}^{'})$ is ISS \cite{jiang2001input}. According to \cite{munoz2020convergence} and \cite{goulart2008input}, there exists a continuous function $V_b$: $\mathbb{X}\rightarrow\mathbb{R_+}$ that is an ISS-Lyapunov function. Therefore, the nominal model $\bm{f}(\bm{x}_{k},\bm{u}_{k}, \bm{e_{fk}})$ with $\mathbb{E}{(w_k)}=0$ is ISS. \hfill $\blacksquare$

\begin{comment}
The dynamic system (\ref{quadrotor_nonlinear_discrete_system}) has probabilistic constraints for state and input:
\begin{equation}
\begin{aligned}
&(\bm{x}_k,\bm{u}_k)\in \mathbb{Z}\\
&\mathbb{P}\left\{ g_j^{\rm{T}}\bm{x}^{k+1}\leq h_j\right\}\geq p_j, j\in\mathbb{N}_{[1,n_c]}\\
&\mathbb{P}\left\{ s_j^{\rm{T}}\bm{u}^{k+1}\leq i_j\right\}\geq q_j, j\in\mathbb{N}_{[1,n_u]}\\
\end{aligned}
\label{probabilistic_constraints}
\end{equation}
where $\mathbb{Z}\subseteq\mathbb{R}^{n}\times\mathbb{R}^{n_u}$ is a compact and convex set which contains the origin in its interior. $g_j\in \mathbb{R}^n$, $s_j\in \mathbb{R}^{n_u}$, $h_j\in \mathbb{R}$, $i_j\in \mathbb{R}$ and $p_j, q_j \in [0,1]$. $n_c$ is the number of probabilistic constraints.
\end{comment}

\section{QuaDRED-SMPC}
In this section, we present the proposed QuaDRED-SMPC control framework. Traditional non-interacted methods, e.g., Gaussian Process \cite{torrente2021data} and RDRv \cite{faessler2017differential}, are insufficient for quadrotor dynamic disturbance estimation. This work addresses the limitation, and proposes a novel and feasible disturbance estimation with continuous environmental interactions for variable winds.

\subsection{Quantile-approximation-based Distributional Reinforced-disturbance-estimation \label{D_RL}}
\begin{comment}
Maximize the expectation of value-action function $Q$ - a.k.a., Q-function - is the aim of traditional RL \cite{sutton2018reinforcement}, the fundamental of which is the iterative use of the Bellman equation. In this paper, we define the \textit{Bellman operator} $\mathcal{T}^\pi$ as:
\begin{equation}
\mathcal{T}^{\pi} Q(\bm{s},\bm{a}):=\mathbb{E}R(\bm{s},\bm{a})+\gamma\underset{\bm{p},\pi}{\mathbb{E}}Q(\bm{s'},\bm{a'})
\label{Bellman_traditional_RL}
\end{equation}
where $\mathbb{E}$ is an \textit{expectation operator} over the distribution. As usual, $\bm{s}\in S$, $\bm{a}\in A$, $\bm{p}\in P: S\times A\times S$, $R$ and $\bm{\gamma}\in[0,1]$ are the state vector, action vector, transition probability, immediate reward function and discount rate, respectively, in a tuple Markov Decision Process \cite{puterman2014markov}: $MDP:=\left \langle S,A,P,R,\gamma \right \rangle$. $\pi$ is a stationary policy mapping one state $s \in S$ to one action $a \in A$.
\end{comment}

We consider a distributional Bellman equation \cite{bellemare2017distributional}, the aim of which is different from traditional RL, i.e., maximizing the expectation of value-action function $Q$. In the policy evaluation setting, given a deterministic policy $\pi$, the state-action distribution $Z^\pi$ and the \textit{Bellman operator} $\mathcal{T}^\pi$ are defined as \cite{bellemare2017distributional,dabney2018distributional}:
\begin{equation}
\mathcal{T}^{\pi} Z(\bm{s},\bm{a})\overset{D}{:=} R(\bm{s},\bm{a})+\gamma Z(\bm{s'},\bm{a'})
\label{Bellman_Pred_D_RL}
\end{equation}
where $\bm{s}\in S$, $\bm{a}\in A$, $\bm{p}\in P: S\times A\times S$, $R$ and $\bm{\gamma}\in[0,1]$ are the state vector, action vector, transition probability, immediate reward function and discount rate, respectively, in a tuple Markov Decision Process \cite{puterman2014markov}: $MDP:=\left \langle S,A,P,R,\gamma \right \rangle$. $\pi$ is a stationary policy mapping one state $s \in S$ to one action $a \in A$. In the control setting, a distributional \textit{Bellman optimality operator} $\mathcal{T}$ with quantile approximation is also proposed in \cite{dabney2018distributional}:
\begin{equation}
\begin{aligned}
\mathcal{T} Z(\bm{s},\bm{a})\overset{D}{:=} R(\bm{s},\bm{a})&+\gamma Z(\bm{s'},{\rm{arg}} \underset{a'}{max}\underset{\bm{p},R}{\mathbb{E}}[Z(\bm{s'},\bm{a'})])\\
Z_\theta(\bm{s},\bm{a})&:=\frac{1}{N}\sum\limits_{i=1}^{N} \delta_{q_i(\bm{s},\bm{a})}
\label{Bellman_Control_D_RL}
\end{aligned}
\end{equation}
where $Z_\theta \in Z_Q$ is a quantile distribution mapping one state-action pair $(s,a)$ to a uniform probability distribution supported on $q_i$. $Z_Q$ is the space of quantile distribution within $N$ supporting quantiles. $\delta_z$ denotes a Dirac with $z \in \mathbb{R}$. The state-action value $Q$ is then approximated by $Q_{j|K} \overset{D}{:=}\frac{1}{K}\sum\limits_{k=(j-1)K+1}^{(j-1)K+K}q_k(s,a)$. These quantile approximations -i.e., $\left\{q_i\right\}$ - are operated based on Quantile Huber Loss.

The Wasserstein Metric, also known as the Mallows metric, is a true probability metric with no disjoint support issues, and therefore suitable to calculate the metric distance between the target $\mathcal{T}^{\pi} Z$ and the prediction $Z$. A contraction is proved in \cite{dabney2018distributional} over the Wasserstein Metric:
\begin{equation}
\overset{-}{d}_{\infty}(\Pi_{W_1}\mathcal{T}^{\pi}Z_1,\Pi_{W_1}\mathcal{T}^{\pi}Z_2)\leq \overset{-}{d}_{\infty}(Z_1,Z_2)
\label{contraction_WM}
\end{equation}
where $W_p$, $p\in[1,\infty]$ denotes the $p$-Wasserstein distance. $\overset{-}{d}_{p}:={\rm{sup}}W_p(Z_1,Z_2)$ denotes the maximal form of the  $p$-Wasserstein metrics. $\Pi_{W_1}$ is a quantile approximation under the minimal 1-Wasserstein distance $W_1$.

The aim of our proposed QuaDRED is to track the trajectory reference $\bm{x}_{m,t}$ generated from Kino-JSS \cite{wang2022kinojgm} accurately,  therefore the immediate reward $r_{t+1}$ is defined as:
\begin{equation}
r_{t+1}=-(\bm{x}_{t}-\bm{x}_{m,t})^{\rm{T}}H_1 (\bm{x}_{t}-\bm{x}_{m,t})-\bm{u}_{t}^{\rm{T}}H_2\bm{u}_{t}
\label{reward}
\end{equation}
where $H_1$ and $H_2$ are positive definite matrices. Then we use DDPG architecture \cite{lillicrap2015continuous} for the continuous and high-dimensional disturbance estimation.
\begin{comment}
At each time step, the parameter set $\theta^Q$ of the critic $Q$ are updated to minimize as:
\begin{equation}
\frac{1}{2}(r_{t+1}+\gamma Q(s_{t+1},\mu_{t+1})-Q(s_t,a_t))^2
\label{update_DDPG}
\end{equation}
where $r_t \in R$ is the immediate reward at time step $t$. $\mu$ denotes the actor, parameterized by $\theta^\mu$, in DDPG.
\end{comment}

\begin{comment}
Referred to \cite{zhang2019quota}, the policy gradient is defined as:
\begin{equation}
{\nabla_a Q_{j|K}(s_t,a)|_{a=\mu^{j}(s_t)}}{\nabla_{\theta^\mu}} \mu^{j}(s_t)
\label{pg_DDPG}
\end{equation}
\end{comment}
\begin{comment}
\subsection{Control Parameterization}
To get a convex optimization problem, we construct the predicted control policy as an affine function of future disturbances \cite{goulart2008input}.

The affine feedback prediction control policy was proposed in \cite{lofberg2003approximations} as:
\begin{equation}
\bm{u}_{i|t}=\sum\limits_{j=0}^{i-1} \bm{M}_{i,j|t}\bm{w}_{j|t}+\bm{v}_{i|t}
\label{affine_fb_first}
\end{equation}
where feedback gains $\bm{M}_{i,j|t}\in\mathbb{R}, i\in\mathbb{N}_{[1,N-1]}, j\in\mathbb{N}_{i-1}$ and the open-loop control sequence $\bm{v}_{i|t} \in \mathbb{R}, i\in\mathbb{N}_{N-1}$ are decision variables at each time step $t$.
\end{comment}

\subsection{Control Parameterization}
A SADF \cite{zhang2021stochastic}, as shown in \autoref{affine_fb_SADF}, is an equivalent and tractable formulation of the original affine feedback prediction control policy proposed in \cite{jiang2001input}. More importantly, the SADF has fewer decision variables which can decrease computational complexity and improve calculation efficiency.
\begin{equation}
\bm{u}_{i|t}=\sum\limits_{k=0}^{i-1} \bm{M}_{i-k|t}\bm{w}_{k|t}+\bm{v}_{i|t}
\label{affine_fb_SADF}
\end{equation}
where the $\bm{M}_{t}$ is a lower block diagonal Toeplitz structure. $i\in\mathbb{N}_{[1,N-1]}, j\in\mathbb{N}_{i-1}$ and the open-loop control sequence $\bm{v}_{i|t} \in \mathbb{R}, i\in\mathbb{N}_{N-1}$ are decision variables at each time step $t$.
\begin{comment}
\begin{equation}
\begin{aligned}
&\bm{v}_t := \begin{bmatrix} v_{0|t}^{\rm{T}} \\
v_{1|t}^{\rm{T}} \\
v_{2|t}^{\rm{T}} \\
\vdots \\
v_{N-1|t}^{\rm{T}} \end{bmatrix},
\bm{M}_t := \begin{bmatrix} 0 & 0 & \hdots & 0 \\
M_{1|t} & 0 & \hdots & 0 \\
M_{2|t} & M_{1|t} & \hdots & 0 \\
\vdots & \ddots & \vdots & \vdots \\
M_{N-1|t} & \hdots & M_{1|t} & 0 \end{bmatrix} \\
\end{aligned}
\label{affine_decision_variables}
\end{equation}
\end{comment}

According to \cite{goulart2008input}, the predicted cost can be transformed as:
\begin{equation}
\begin{aligned}
&\mathcal{L}(\bm{x}_t,\bm{u}_t)=\mathcal{L}_N(\bm{x}_t,\bm{M}_t,\bm{v}_t)\\
&=\left \|H_{x}x+H_{u}\bm{v}\right \|^{2}_2+\mathbb{E}[\left \|(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{w}\right \|^{2}_2]
\end{aligned}
\label{predicted_cost}
\end{equation}
where $H_{x}$ and $H_{u}$ are coefficient matrices which are constructed from \autoref{predicted_cost}. $\mathcal{G}:=I_N \bigotimes G$ denotes Kronecker product of matrices $I_N$ and $G$. For the convexity guarantee, the matrix $(H_{u}{\bm{M}}\mathcal{G}+H_{w})$ is positive semidefinite.
\begin{comment}
\begin{equation}
\begin{aligned}
&P:=Q+A^{\rm{T}}PA-K^{\rm{T}}(R+B^{\rm{T}}PB)K\\
&K:=-(R+B^{\rm{T}}PB)^{-1}B^{\rm{T}}PA\\
\end{aligned}
\label{algebraic_Riccati}
\end{equation}
\end{comment}
Thus, the optimal control problem, reformulated by SADF (\autoref{affine_fb_SADF}) is as follows:
\begin{equation}
\begin{aligned}
\underset{\bm{M}_t,\bm{v}_t}{\rm{min}} \mathcal{L}_N(\bm{x}_t,\bm{M}_t&,\bm{v}_t), {\rm{s.t.}} \forall w_{i|t}\in \mathbb{W}, \forall i\in \mathbb{N}_{N-1}\\
\rm{subject}\quad \rm{to}\quad &\bm{x}_{t} = \bm{A}x_{0|t} + \bm{B}\bm{u}_{t} + \bm{G}\bm{w}_{t}\\
&\bm{u}_{i|t}=\sum\limits_{k=0}^{i-1} \bm{M}_{i-k|t}\bm{w}_{k|t}+\bm{v}_{i|t}\\
&H_{u}{\bm{M}}\mathcal{G}+H_{w} \geq 0\\
&(\bm{x}_t,\bm{u}_t)\in \mathbb{Z}\\
&\bm{x}_{N|t}\in \mathbb{X}_{f}\\
&\bm{x}_{0|t}=\bm{x}_t\\
\end{aligned}
\label{SMPC_Formulation}
\end{equation}

The optimal control problem is a strictly convex quadratic program or second-order cone program (SOCP) if $\mathbb{W}$ is a polytope or ellipsoid when $\mathbb{Z}$ and $\mathbb{X}_{f}$ are polytopic \cite{goulart2008input}. In this case, this problem can be seen as deterministic MPC with nonlinear constraints, which can be solved by some nonlinear MPC solvers, e.g., CasADi \cite{andersson2019casadi} and ACADOS \cite{verschueren2018towards}.

\subsection{Quantile-approximation-based Distributional Reinforced-disturbance-estimation for SMPC}
The objective of this work is to design a quadrotor controller achieving accurate aerodynamic effect estimation, which we define as combined wind estimation and aerodynamic disturbance estimation, for tracking the reference state $\bm{x}_{ref}$ of the nominal model (\autoref{quadrotor_dynamics}). The overall control framework for the quadrotor is shown in Fig.~\ref{RL_Control_framework}. The SADF in stochastic MPC and QuaDRED are shown in Algorithm~{\autoref{SADF_SMPC}} and Algorithm~{\autoref{DRL_QuaDRED}}, respectively.
\begin{figure}[]
  \centering
  \includegraphics[scale=0.54]{RL_Control_framework.jpg}
  \caption{QuaDRED-SMPC}
  \label{RL_Control_framework}
\end{figure}

\begin{algorithm}[t]
\caption{SADF-SMPC}
\begin{algorithmic}[1]
\STATE \textbf{Get}:\\
    - the reference data $\bm{x}_{ref}$ from the quadrotor trajectory planning and generation module, i.e., Kino-JSS \cite{wang2022kinojgm}\\
    - the measurement state $\bm{x}_k$ from on-board sensors\\
    - the wind estimation $\bm{e}_{fk}$ from VID-Fusion \cite{ding2020vid}\\
\STATE \textbf{Initialize}:\\
    - the parameters $\theta^\mu$ and $\theta^Q$ for the actor $\mu$ and the critic $Q$, respectively\\
    - the decision variables $\bm{M}_{0}$ and $\bm{v}_{0}$ in \autoref{affine_fb_SADF}\\
    - the initial state $\bm{s}_{0}$\\
\FOR{each sampling timestamps $k$}
\STATE \textbf{Repeat}
\STATE $\bm{s}_k \gets \left [ \bm{x}_k, \bm{e}_{fk} \right ]$
\STATE Select an action vector $\bm{w}_k \gets [w_{0|k}^{\rm{T}}, w_{1|k}^{\rm{T}}, ..., w_{N|k}^{\rm{T}}]^{\rm{T}}$ from $\bm{w}_k=\bm{a}_k \gets \mu(\bm{a}_k|\bm{s}_k)$ in QuaDRED (Algorithm~\autoref{DRL_QuaDRED})
\STATE $\bm{u}_{i|k} \gets \sum\limits_{l=0}^{i-1} \bm{M}_{i-l|k}\bm{w}_{l|k}+\bm{v}_{i|k}$
\STATE $\bm{u}_{k} \gets \bm{u}_{0|k}$, $\bm{w}_{k} \gets \bm{w}_{0|k}$
\STATE $\bm{x}_{k} \gets \bm{A}x_{0|k} + \bm{B}\bm{u}_{k} + \bm{G}\bm{w}_{k}$
\STATE Solve optimization problem \autoref{SMPC_Formulation} with nonlinear MPC solver\\
\STATE \textbf{Until} convergence
\STATE $\bm{u}_{k} \gets \bm{v}_{0|k}$
\STATE $\bm{x}_{k+1}$, $\bm{e}_{fk+1} \gets {\rm{Real Quadrotor}}(\bm{u}_{k})$
\STATE $\bm{s}_{k+1} \gets [\bm{x}_{k+1}, \bm{e}_{fk+1}]$
\STATE $\bm{x}_{ref} \gets$ Kino-JSS
\STATE $k \gets k+1$
\ENDFOR
\end{algorithmic}
\label{SADF_SMPC}
\end{algorithm}

\begin{algorithm}[]
\caption{QuaDRED}
\hspace*{0.02in} {\bf Input:}
$\bm{s}_k$, $\bm{s}_{k+1}$, $\bm{u}_{k}$, $\theta^\mu$, $\theta^Q$\\
\hspace*{0.02in} {\bf Output:}
$\bm{a}_k$
\begin{algorithmic}[1]
\STATE \textbf{Initialize}:\\
    - $\theta^{\mu^t} \gets \theta^\mu$, $\theta^{Q^t} \gets \theta^Q$ update the target parameters from the predicted parameters\\
    - the replay memory $D \gets D_{k-1}$\\
    - the batch $B$, and its size\\
    - a small threshold $\xi \in \mathbb{R_+}$ \\
    - the random option selection probability $\epsilon$
    - the option termination probability $\beta$ \\
    - quantile estimation functions $\left \{ q_i   \right \}_{i=1,...,N}$\\
\STATE \textbf{Repeat}
\FOR{each sampling step from $D$}
\STATE Select a candidate option $z_k$ from $\left \{z^0, z^1, ..., z^M   \right \}$
\STATE $z_k \gets
\begin{cases}
z_{k-1} & {w.p.\;1-\beta} \\
{\rm{random \; option}} & {w.p.\;\beta\epsilon} \\
{\rm{arg max}}_{z}Q(\bm{s}_k,z) & {w.p.\;\beta(1-\epsilon)}
\end{cases}$
\STATE Execute $w_k$, get reward $r_k$ and the next state $\bm{s}_{k+1}$
\STATE $D.\bf{Insert}([\bm{s}_{k}, \bm{u}_{k}, r_k, \bm{s}_{k+1}])$
\STATE $B \gets D.\bf{sampling}$
\STATE $y_{k,i} \gets \rho_{\tau_i}^{\mathcal{K}}(r_k+\gamma q_i'(\bm{s}_{k+1},w_k^*)$
\STATE $J_{\theta^\mu} \gets \frac{1}{N}\sum\limits_{i=1}^{N}\sum\limits_{i'=1}^{N}[y_{k,i'}-q_i(\bm{s}_k,w_k)]$
\STATE $y \gets \beta {\rm{arg max}}_{z'}Q(\bm{s}_{k+1},z')+(1-\beta)Q(\bm{s}_{k+1}.z_k)$
\STATE $J_{\theta^Q} \gets (r_t+\gamma y-Q(\bm{s}_t,z_t))^2$
\STATE $\theta^\mu \gets \theta^\mu - l_{\mu}\nabla_{\theta^\mu} J_{\theta^\mu}$
\STATE $\theta^Q \gets \theta^Q - l_{\theta} \nabla_{\theta^Q} J_{\theta^Q}$
\ENDFOR
\STATE \textbf{Until} convergence, i.e., $J_Q^\theta < \xi$
\end{algorithmic}
\label{DRL_QuaDRED}
\end{algorithm}

\section{Properties of QuaDRED-SMPC}
In this section, the properties of the proposed control framework QuaDRED-SMPC are analyzed, including convergence of QuaDRED and stability guarantees of the Controller SADF-SMPC.
\subsection{Convergence Analysis of QuaDRED}
We present the following \textit{Proposition 3} and \textit{Proposition 4} on the convergence analysis for the Distributional RL (QuaDRED) in Section \uppercase\expandafter{\romannumeral4}-A.

\textit{Lemma 2} (\cite{bellemare2017distributional}): The \textit{Bellman operator} $\mathcal{T}^\pi$ is a $p$-contraction under the $p$-Wasserstein metric $\overset{-}{d}_{p}$.
\begin{comment}
\textit{Proof}: This proof can be found in \textit{Lemma 3} of \cite{bellemare2017distributional}.
\end{comment}

\textit{Lemma 2} suggests that an effective way in practice to minimize the Wasserstein distance between a distribution $Z$ and its Bellman update $\mathcal{T}^{\pi}Z$ can be found in \autoref{Bellman_Pred_D_RL}, which attempts iteratively to minimize the $L2$ distance between $Z$ and $\mathcal{T}^{\pi}Z$ in Temporal Difference learning.

\textit{Proposition 3 (Policy Evaluation)}: Let $\Pi_{W_1}$ be a quantile approximation under the minimal 1-Wasserstein distance $W_1$, $\mathcal{T}^\pi$ be the \textit{Bellman operator} under a deterministic policy $\pi$ and $Z_{k+1}(\bm{s},\bm{a})=\Pi_{W_1}\mathcal{T}^{\pi} Z_k(\bm{s},\bm{a})$. The sequence $Z_{k}(\bm{s},\bm{a})$ converges to a unique fixed point $\overset{\sim}{Z_\pi}$ under the maximal form of $\infty$-Wasserstein metric $\overset{-}{d}_{\infty}$.

\textit{Proof}:
\autoref{contraction_WM} implies that the combined operator $\Pi_{W_1} \mathcal{T}^\pi$ is an $\infty$-contraction \cite{dabney2018distributional}. We conclude using Banach’s fixed point theorem that $\mathcal{T}^\pi$ has a unique fixed point, i.e., $\overset{\sim}{Z_\pi}$. Furthermore, \autoref{Bellman_Control_D_RL} implies that all moments of $Z$ are bounded. Therefore, we conclude that the sequence $Z_{k}(\bm{s},\bm{a})$ converges to $\overset{\sim}{Z_\pi}$ in $\overset{-}{d}_{\infty}$ for $p\in[1,\infty]$. \hfill $\blacksquare$

\textit{Proposition 4 (Policy Improvement)}: Let $\bm{\pi}_{\bm{old}}$ be an old policy, $\bm{\pi}_{\bm{new}}$ be a new policy and $Q(s, a)=\mathbb{E}[Z(s,a)]$ in \autoref{Bellman_Control_D_RL}. There exists $Q^{\bm{\pi}_{\bm{new}}}(s, a) \geq Q^{\bm{\pi}_{\bm{old}}}(s, a)$, $\forall s\in \mathcal{S}$ and $\forall a\in \mathcal{A} $.

\textit{Proof}: Based on \autoref{Bellman_Control_D_RL}, there exists:
\begin{equation}
\begin{aligned}
V^\pi(s_{t})&=\mathbb{E}_\pi {Q^\pi(s_{t},\pi(s_{t}))}\\
&\leq \underset{a\in \mathcal{A}}{\rm{max}} \mathbb{E}_\pi {Q^\pi(s_{t},a)} \\
&=\mathbb{E}_{\pi'}{Q^\pi(s_{t},{\pi'}(s_{t}))}
\end{aligned}
\label{Q_V}
\end{equation}
where $\mathbb{E}_\pi[\cdot]=\sum_{a\in A}\bm{\pi}(a|s)[\cdot]$, and $V^\pi(s)=\mathbb{E}_\pi \mathbb{E}[Z_k(s,a)]$ is the value function. According to \autoref{Q_V} and \autoref{Bellman_Control_D_RL}, it yields:
\begin{equation}
\begin{aligned}
Q^{\bm{\pi}_{\bm{old}}} &=Q^{\bm{\pi}_{\bm{old}}}(s_{t},\bm{\pi}_{\bm{new}}(s_{t})) \\
&= r_{t+1}+\gamma \mathbb{E}_{s_{t+1}} \mathbb{E}_{\bm{\pi}_{\bm{old}}} Q^{\bm{\pi}_{\bm{old}}}(s_{t+1},{\bm{\pi}_{\bm{old}}}(s_{t+1}))\\
&\leq r_{t+1}+\gamma \mathbb{E}_{s_{t+1}} \mathbb{E}_{\bm{\pi}_{\bm{new}}}{Q^{\bm{\pi}_{\bm{old}}}(s_{t+1},{\bm{\pi}_{\bm{new}}}(s_{t+1}))}\\
&\leq r_{t+1}+\mathbb{E}_{s_{t+1}} \mathbb{E}_{\bm{\pi}_{\bm{new}}} [\gamma r_{t+2} \\ &+ {\gamma^2} \mathbb{E}_{s_{t+2}}{Q^{\bm{\pi}_{\bm{old}}}(s_{t+2},{\bm{\pi}_{\bm{new}}}(s_{t+2}))}|]\\
&\leq r_{t+1}+\mathbb{E}_{s_{t+1}} \mathbb{E}_{\bm{\pi}_{\bm{new}}} [\gamma r_{t+2} + {\gamma^2}r_{t+3} + ...]\\
&= r_{t+1}+\mathbb{E}_{s_{t+1}} V^{\bm{\pi}_{\bm{new}}}(s_{t+1})\\
&=Q^{\bm{\pi}_{\bm{new}}}
\end{aligned}
\label{policy_improvement}
\end{equation} \hfill $\blacksquare$

Given \textit{Proposition 3} and \textit{Proposition 4}, we can now analyze the convergence of the QuaDRED.

\textit{Theorem 5 (Convergence)}: Let $\bm{\pi}^{\bm{i}}$ be the policy in the $i$-th policy improvement, $i=1,2,...,\infty$, and $\bm{\pi}^{\bm{i}} \rightarrow \pi^{*}$ when $i\rightarrow\infty$. There exists $Q^{\bm{\pi}^{*}}(s, a) \geq Q^{\bm{\pi}^{\bm{i}}}(s, a)$, $\forall s\in \mathcal{S}$ and $\forall a\in \mathcal{A}$.

\textit{Proof}: Since \textit{Proposition 4} suggests $Q^{\bm{\pi}_{\bm{i+1}}}(s, a) \geq Q^{\bm{\pi}_{\bm{i}}}(s, a)$, the sequence $Q^{\bm{\pi}_{\bm{i}}}(s, a)$ is monotonically increasing where $i \in \mathbb{N}$ is a the policy iteration step. Furthermore, \textit{Lemma 2} implies that the the state-action distribution $Z$ over $\mathbb{R}$ has bounded $p$-th moment, so the first moment of $Z$, i.e., $Q^{\bm{\pi}_{\bm{i}}}(s, a)$, is upper bounded. Therefore, the sequence $Q^{\bm{\pi}_{\bm{i}}}(s, a)$ converges to an upper limit $Q^{\bm{\pi}_{*}}(s, a)$ with $\forall s\in \mathcal{S}$ and $\forall a\in \mathcal{A}$. \hfill $\blacksquare$

\subsection{Stability Guarantee of the Controller}
In this subsection, the closed-loop stability of QuaDRED-SMPC control framework will be demonstrated. The closed-loop stability is analyzed under the Lipschitz Lyapunov function \cite{rifford2000existence} to guarantee ISS. Before the closed-loop stability analysis, the convexity and Lipschitz continuity of the cost function $\mathcal{L}_{\bm{N}}(\bm{x}_t,\bm{M}_t,\bm{v}_t)$ are introduced in \textit{Proposition 6} and \textit{Proposition 7}, respectively. Since the output of QuaDRED is non-zero-mean and bounded values, which are different from the assumption of zero-mean disturbances in most previous work (\cite{zhang2021stochastic} and \cite{goulart2008input}), the following proofs are all based on the non-zero-mean and bounded disturbances.

We first define an optimal control policy based on the affine disturbance feedback control law:
\begin{equation}
(\bm{M}^{*}(x),\bm{v}^{*}(x)):=\underset{(\bm{M},\bm{v})\in \mathcal{V}_N}{\rm{min}} \mathcal{L}_N(\bm{x},\bm{M},\bm{v})
\label{optimal_control_policy}
\end{equation}

where $\mathcal{V}_N$ is the set of feasible policies, and $(\bm{M}^{*}(x),\bm{v}^{*}(x))$ is a optimal control policy group. The optimal value function $\mathcal{L}^{*}_N(x)$ under the affine disturbance feedback control law is defined as:
\begin{equation}
\mathcal{L}^{*}_N(x):=\underset{(\bm{M},\bm{v})\in \mathcal{V}_N}{\rm{min}} \mathcal{L}_N(\bm{x},\bm{M},\bm{v})
\label{optimal_value}
\end{equation}

Then we demonstrate that the optimal value function $\mathcal{L}_N(x)$ is convex (see \textit{Proposition 6}), so that  \autoref{optimal_value} can be operated as a convex optimization problem.

\textit{Proposition 6}: The function $\mathcal{L}_N(\bm{x},\bm{M},\bm{v})$ is convex.

\textit{Proof}: In \autoref{predicted_cost}, the second term  $\mathbb{E}[\left \|(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{w}\right \|^{2}_2]$, i.e., the expected value of a quadratic form with respect to the vector-valued random variable $\bm{w}$, is equal to:
\begin{equation}
\begin{aligned}
&\mathbb{E}[\left \|(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{w}\right \|^{2}_2]=\mathbb{E}[{\rm{tr}}((H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{w} \bm{w}^T)]\\
&={\rm{tr}}((H_{u}{\bm{M}}\mathcal{G}+H_{w})\mathbb{E}[\bm{w} \bm{w}^T])\\
&={\rm{tr}}((H_{u}{\bm{M}}\mathcal{G}+H_{w})({\rm{Cov}}(\bm{w})+\bm{\mu} \bm{\mu}^T))\\
&={\rm{tr}}(\bm{C}^{\frac{1}{2}}_{\bm{w}}(H_{u}{\bm{M}}\mathcal{G}+H_{w})^{T}(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{C}^{\frac{1}{2}}_{\bm{w}})\\
&+\bm{\mu}^T(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{\mu}
\end{aligned}
\label{convex_simplied_w}
\end{equation}
where $\rm{tr}(\cdot)$ denotes the trace of a square matrix. $\bm{\mu}=\mathbb{E}(\bm{w})$ is the expected value of $\bm{w}$, and $\bm{C}_{\bm{w}}={\rm{Var}}(\bm{w})$ is the variance-covariance matrix of $\bm{w}$. Therefore, $\mathcal{L}_N(x)$ can be written as:
\begin{equation}
\begin{aligned}
&\mathcal{L}_N(x)=\left \|H_{x}x+H_{u}\bm{v}\right \|^{2}_2+\left \|\bm{\mu}\right \|^{2}_{(H_{u}{\bm{M}}\mathcal{G}+H_{w})}\\
&+{\rm{tr}}(\bm{C}^{\frac{1}{2}}_{\bm{w}}(H_{u}{\bm{M}}\mathcal{G}+H_{w})^{T}(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{C}^{\frac{1}{2}}_{\bm{w}})
\end{aligned}
\label{convex_l_n}
\end{equation}
where $\left \|x\right \|_{P}$ denotes weighted $2$–norm of the vector $x$. \autoref{convex_l_n} is convex since it consists of convex functions of vector and matrix norms.\hfill $\blacksquare$

\textit{Proposition 7}: The function $\mathcal{L}^{*}_N(\bm{x},\bm{M},\bm{v})$ is Lipschitz continuous.

\textit{Proof}: The cost function $\mathcal{L}_N(\bm{x},\bm{M},\bm{v})$ is proved to be convex in \textit{Proposition 6} so that $\mathcal{L}^{*}_N(\bm{x},\bm{M},\bm{v})$ is convex if $\mathcal{V}_N$ has a non-empty interior (\textit{Proposition 1} of \cite{goulart2008input}). $Z$ is a compact (closed and bounded) set so that $\mathcal{L}^{*}_N(\bm{x},\bm{M},\bm{v})$, defined under the compact space $Z$, is piecewise quadratic (\textit{Corollary 4.6} of \cite{goulart2007affine}). Therefore $\mathcal{L}^{*}_N(\bm{x},\bm{M},\bm{v})$ is a Lipschitz continuity function. \hfill $\blacksquare$
\begin{comment}
We can now give a reasonable condition for ISS under Lipschitz continuity. \textit{Lemma 9} is an ISS criterion for the Lipschitz continuous functions $\mathcal{L}^{*}_N(\bm{x},\bm{M},\bm{v})$ if $\mathbb{W}$ is a polytope or ellipsoid, and $\mathbb{Z}$ and $\mathbb{X}_{f}$ are polytopic.

\textit{Lemma 9 (See \cite{goulart2008input})}: Let $\mathcal{X} \in\mathbb{R}$ have the origin in its interior. The dynamic system (\ref{quadrotor_nonlinear_discrete_system}) is ISS if there exists $\mathcal{H}_{\infty}$-functions $\alpha_1(\cdot)$, $\alpha_2(\cdot)$ and $\alpha_3(\cdot)$, and a Lipschitz continuous function $V$: $\mathcal{X}\rightarrow\mathbb{R_+}$ such that
\begin{subequations}
\begin{align}
\alpha_1(&\left \| x \right\|) \leq V(x) \leq \alpha_2(\left \| x \right\|), \forall x\in\mathbb{R}^n \label{ISS_Lipschitz_continuity_a}\\
&V(f(x,0)) - V(x) \leq -\alpha_3(\left \| x \right\|) \label{ISS_Lipschitz_continuity_b}
\end{align}
\end{subequations}
The function $V(\cdot)$ is an ISS-Lyapunov function, and $f(x,w)$ is Lipschitz continuous on $\mathcal{X} \times \mathcal{W}$.

\textit{Proof}: Let $\mathcal{L}_v$ be the Lipschitz constant of $V(\cdot)$, and  $\mathcal{L}_f$ be the Lipschitz constant of $f(\cdot,\cdot)$. We have $||V(f(x,w)) - f(x,0))|| \leq \mathcal{L}_v||f(x,w)-f(x,0)|| \leq \mathcal{L}_v\mathcal{L}_f||w|| $. Furthermore we have $V(f(x,0)) - V(x) \leq -\alpha_3(\left \| x \right\|)+\sigma (\left \| x \right\|)$ where $V(x)$ is an ISS-Lyapunov function. Therefore there exists $\sigma(x)=\mathcal{L}_v\mathcal{L}_f x$ such that the dynamic system (\ref{quadrotor_nonlinear_discrete_system}) is ISS. \hfill $\blacksquare$
\end{comment}

The above results lead directly to our final result:

\textit{Theorem 9}: Let $\mathcal{W}$, $\mathcal{Z}$ and $\mathcal{X}_{f}$ be polytopes. The closed-loop system (\autoref{linear_discrete_system}) under the SADF control law $\bm{u}_{i|t}$ (in \autoref{affine_fb_SADF}) is ISS. The ISS is also guaranteed in such cases: the stochastic disturbance $\bm{w}_t$ is i.i.d. a bounded and non-zero-mean distribution, i.e., $\mathbb{E}{(w_k)}\neq0$.

\textit{Proof}: According to \textit{Proposition 6} and \textit{Proposition 7}, we first state that the optimal value function $\mathcal{L}^{*}_N(x)$ is a Lipschitz continuity function. The key is then to prove that there exists a Lipschitz continuous function, i.e., $\mathcal{L}^{*}_N(x)$, to satisfy the Lipschitz-ISS criterion (\textit{Proposition 4.15} in \cite{goulart2008input}).

According to \textit{Proposition 1}, there exists a baseline control law $\bm{u}_b$ ensuring ISS under zero-mean distribution disturbance. Let $V_b (x)=V^*_{Nb}(x)-V^*_{Nb}(0)$ be the Lipschitz continuous Lyapunov function \cite{goulart2008input}, where $V^*_{Nb}(x)$ is the optimal value function under the baseline control law $\bm{u}_b$. There exists:
\begin{subequations}
\begin{align}
\alpha_1(\left \| x \right\|) \leq &V_b(x) \leq \alpha_2(\left \| x \right\|) \label{ISS_baseline_control_law_a}\\
V_b(f(x,0)) - &V_b(x) \leq -\alpha_3(\left \| x \right\|) \label{ISS_baseline_control_law_b}
\end{align}
\end{subequations}

Let $V(x)=\mathcal{L}^*_{N}(x)-\mathcal{L}^*_{N}(0)$, where $\mathcal{L}^*_{N}(x)$ is optimal value function under the affine disturbance feedback control law with bounded and non-zero-mean distribution disturbance (in \autoref{optimal_value}). According to \autoref{convex_l_n}, $L^*_{N}(x)$ is shown as:
\begin{equation}
\begin{aligned}
&\mathcal{L}^*_{N}(x)={\rm{min}} \left.\{\mathcal{L}_{N}(x) \right.\}\\
&={\rm{min}}\left.\{ \left \|H_{x}x+H_{u}\bm{v}\right \|^{2}_2 +\left \|\bm{\mu}\right \|^{2}_{(H_{u}{\bm{M}}\mathcal{G}+H_{w})} \right.\\
&\left.+{\rm{tr}}(\bm{C}^{\frac{1}{2}}_{\bm{w}}(H_{u}{\bm{M}}\mathcal{G}+H_{w})^{T}(H_{u}{\bm{M}}\mathcal{G}+H_{w})\bm{C}^{\frac{1}{2}}_{\bm{w}}) \right.\}\\
&=V^*_{Nb}(x)+{\rm{min}}\left.\{\left \|\bm{\mu}\right \|^{2}_{(H_{u}{\bm{M}}\mathcal{G}+H_{w})} \right.\}
\end{aligned}
\label{ISS_opyimal_value}
\end{equation}
where $\bm{\mu}$ is the expected value of disturbances, which is independent with $\bm{v}$ and $\bm{M}$. $V(x)=\mathcal{L}^*_{N}(x)-\mathcal{L}^*_{N}(0)=V^*_{Nb}(x)-V^*_{Nb}(0)=V_b (x)$. Hence, there exists $\mathcal{H}_{\infty}$-functions $\alpha_1(\cdot)$, $\alpha_2(\cdot)$ such that \autoref{ISS_baseline_control_law_a} holds with $V_b(\cdot)=V(\cdot)=\mathcal{L}^*_{N}(\cdot)-\mathcal{L}^*_{N}(0)$.

To prove $V(\cdot)$ satisfying \autoref{ISS_baseline_control_law_b}, note that $V_b(f(x,0)) - V_b(x)=[V^*_{Nb}(f(x,0))-V^*_{Nb}(0)]-[V^*_{Nb}(x)-V^*_{Nb}(0)]=V^*_{Nb}(f(x,0))-V^*_{Nb}(x)$, so that $V^*_{Nb}(f(x,0))-V^*_{Nb}(x) \leq -\alpha_3(\left \| x \right\|)$. It follows that:
\begin{equation}
\begin{aligned}
V&(f(x,0)) - V(x) \\
&= [\mathcal{L}^*_{N}(f(x,0))-\mathcal{L}^*_{N}(0)]-[\mathcal{L}^*_{N}(x)-\mathcal{L}^*_{N}(0)]\\
&=\mathcal{L}^*_{N}(f(x,0))-\mathcal{L}^*_{N}(x)
\end{aligned}
\label{ISS_V_second_inequality}
\end{equation}
where both $\mathcal{L}^*_{N}(f(x,0))$ and $V^*_{Nb}(f(x,0))$ have $w=0$. The only difference between two control laws is zero-mean or non-zero-mean disturbance distributions so that $\mathcal{L}^*_{N}(f(x,0))=V^*_{Nb}(f(x,0))$. Hence, combining with \autoref{ISS_opyimal_value}, \autoref{ISS_V_second_inequality} can be rewritten as:
\begin{equation}
\begin{aligned}
V&(f(x,0)) - V(x) \\
&=V^*_{Nb}(f(x,0))-V^*_{Nb}(x)-{\rm{min}}\left.\{ \left \|\bm{\mu}\right \|^{2}_{(H_{u}{\bm{M}}\mathcal{G}+H_{w})} \right.\}
\end{aligned}
\label{ISS_V_second_inequality1}
\end{equation}

According to \textit{Proposition 6}, we have ${\rm{min}}\left.\{\left \|\bm{\mu}\right \|^{2}_{(H_{u}{\bm{M}}\mathcal{G}+H_{w})}\right.\} \geq 0$. Then we have:
\begin{equation}
\begin{aligned}
V&(f(x,0)) - V(x) \\
&\leq V^*_{Nb}(f(x,0))-V^*_{Nb}(x) \leq -\alpha_3(\left \| x \right\|)
\end{aligned}
\label{ISS_V_second_inequality2}
\end{equation}

\autoref{ISS_V_second_inequality2} above shows that there exists $\mathcal{H}_{\infty}$-functions $\alpha_3(\cdot)$ such that \autoref{ISS_baseline_control_law_b} holds with $V(\cdot)=\mathcal{L}^*_{N}(\cdot)-\mathcal{L}^*_{N}(0)$. Therefore, $V(x)=\mathcal{L}^*_{N}(x)-\mathcal{L}^*_{N}(0)$ is a Lipschitz continuous Lyapunov function, and the ISS of the closed-loop system (\autoref{linear_discrete_system}) is guaranteed with bounded and non-zero-mean distribution disturbances, i.e., $\mathbb{E}{(w_k)}\neq0$. \hfill $\blacksquare$

\section{Numerical Example}
The performance of our proposed QuaDRED-SMPC framework is evaluated in RotorS \cite{furrer2016rotors}, a UAV software simulator. Based on the benchmark \cite{wang2022kinojgm,zhang2019quota}, the parameters of our proposed framework are summarized in \autoref{model_parameters}.

\begin{table}[!t]

    \caption{Parameters of QuaDRED-SMPC}
    \label{model_parameters}

    \begin{tabular}{c c c}
    \toprule
    \textbf{Parameters} & \textbf{Definition} & \textbf{Values} \\
    \hline
    $l_{\mu}$ & Learning rate of actor & 0.001 \\
    $l_{\theta}$ & Learning rate of critic & 0.001 \\
    $\mu$ & \makecell[c]{Actor neural network: fully connected with two \\hidden layers (128 neurons per hidden layer)} & - \\
    $\theta$ & \makecell[c]{Critic neural network: fully connected with two \\hidden layers (128 neurons per hidden layer)} & - \\
    $D$ & Replay memory capacity & $10^6$ \\
    $B$ & Batch size & 256 \\
    $\gamma$ & Discount rate & 0.998 \\
    - & Training episodes & 1000 \\
    $T_s$ & MPC Sampling period & 50ms \\
    $N$ & Time steps & 20 \\
    \bottomrule
    \end{tabular}

\end{table}

\subsection{Comparative performance of QuaDRED Training}
In the training process, the quadrotor system operates with aerodynamic effects in the horizontal plane in the range [-3,3] ($m/s^2$). The quadrotor state $\bm{x}$ is recorded at 16 Hz. The training process is occurs over 1000 iterations. The matrices $H_1$ and $H_2$ in \autoref{reward} are chosen as $H_1=diag\{ 2.5e^{-2}, 2.5e^{-2}, 2.5e^{-2}, 1e^{-3}, 1e^{-3}, 1e^{-3}, 2.5e^{-3}, \\ 2.5e^{-3}, 2.5e^{-3}, 2.5e^{-3}, 1e^{-5}, 1e^{-5}, 1e^{-5}\}$ and $H_2=diag\{ 1.25e^{-4}, 1.25e^{-4}, 1.25e^{-4}, 1.25e^{-4}\}$, respectively. The learning curves are displayed in Fig.~\ref{learning_curves}, where we show the training performance of DDPG \cite{lillicrap2015continuous}, C51-DDPG \cite{bellemare2017distributional} and our proposed QuaDRED. The performance shows that, although all the three algorithms converge to a long-term return eventually, the two distributional RL approaches, C51-DDPG and QuaDRED, outperform the traditional DDPG RL approach. More importantly, in the training process, our proposed QuaDRED achieves the largest return whilst maintaining the highest convergence speed.

\begin{figure}[]
  \centering
  \includegraphics[scale=0.18]{learning_curves.jpg}
  \caption{Learning curves of three RL algorithm: DDPG, C51-DDPG and QuaDRED. The simulated speed in Gazebo is set as $0.6$.}
  \label{learning_curves}
\end{figure}

\begin{figure}[]
  \centering
  \includegraphics[scale=0.20]{simulation_test.jpg
  \caption{Specific scenarios results: Wind estimation $\bm{e}_{fk}$, position error $(m)$ and Y control input $T_y$ (expressed in body frame).}
  \label{specific_scenarios_data}
\end{figure}

\subsection{Comparative performance of QuaDRED-SMPC under variable aerodynamic effects}
We compare our QuaDRED-SMPC against a state-of-the-art trajectory tracking algorithm, Torrente \cite{torrente2021data}, and interactive approaches, DDPG and C51,  with variable aerodynamic forces added to our simulated environment. The experiments are based on the trained QuaDRED model described in Section \uppercase\expandafter{\romannumeral4}-A. We first set the aerodynamic forces as [0.0, 2.0, 0.0] ($m/s^2$).  Fig.~\ref{specific_scenarios_data} shows the tracking position errors and control input with two opposite heading aerodynamic forces (with the same force [0.0, 2.0, 0.0] ($m/s^2$)), where our proposed QuaDRED-SMPC has the smallest tracking position error, and reacts to the sudden aerodynamic effects sufficiently.

Two larger and more complex forces, i.e., [-2.0, 2.0, 0.0] and [-3.0, 3.0, 0.0] ($m/s^2$), are then used in the scenario shown in Fig.~\ref{simulation_env}. In \autoref{Comparison_of_tracking}, the second approach, `DDPG + MPC', is a combination of DDPG \cite{lillicrap2015continuous} and MPC; and the third approach, `C51 + SMPC', is a combination of C51 \cite{bellemare2017distributional}, DDPG and SMPC. Our results show that interactive approaches are not always better than non-interactive approaches. For example, `DDPG + MPC' has a lower success rate than GP-MPC, whilst there is little difference in the operation time and cumulative tracking error with relatively small aerodynamic forces. However, compared with GP-MPC, our proposed QuaDRED-SMPC achieves improvements of $42\%$, $58\%$ and $71\%$ in operation time, and $86\%$, $75\%$ and $66\%$ in cumulative tracking errors, respectively.

\begin{figure}[]
  \centering
  \includegraphics[scale=0.68]{simulation_env.jpg}
  \caption{The simulation scenario in RotorS: both reference trajectories with/without external forces are generated by Kino-JSS \cite{wang2022kinojgm}.}
  \label{simulation_env}
\end{figure}

\begin{table}[]
\caption{Comparison of Trajectory Tracking under Programmatic External forces}
\label{Comparison_of_tracking}
\begin{center}
\setlength{\tabcolsep}{1.7mm}{
\begin{tabular}{|l|l|c|c|c|}
\hline
Ex. forces                        & Method      & Succ. Rate & Time (s) & Err. (m) \\ \hline
\multirow{4}{*}{{[}0.0, 2.0, 0.0{]}}  & GP-MPC        & 86\%              & 15.41     & 8.22                        \\ \cline{2-5}
                                      & DDPG + MPC & 81\%               & 13.57     & 6.89                         \\ \cline{2-5}
                                      & C51 + SMPC & 90\%             & 9.35   & 1.48                        \\
                                      \cline{2-5}
                                      & QuaDRED-SMPC    & \bf{91\%}             & \bf{9.04}     & \bf{1.19}
                                      \\ \hline
\multirow{4}{*}{{[}-2.0, 2.0, 0.0{]}} & GP-MPC         & 74\%              & 26.11     & 21.80                        \\ \cline{2-5}
                                      & DDPG + MPC & 78\%               & 18.26     & 15.33                        \\ \cline{2-5}
                                      & C51 + SMPC & 89\%              & 12.25    & 6.68                        \\
                                      \cline{2-5}
                                      & QuaDRED-SMPC    & \bf{91\%}             & \bf{10.97}     & \bf{5.60}
                                      \\ \hline
\multirow{4}{*}{{[}-3.0, 3.0, 0.0{]}} & GP-MPC         & 17\%              & 41.98     & 35.27                        \\ \cline{2-5}
                                      & DDPG + MPC & 46\%               & 30.43     & 27.24                        \\ \cline{2-5}
                                      & C51 + SMPC & 75\%               & 15.31    & 14.53                        \\
                                      \cline{2-5}
                                      & QuaDRED-SMPC    & \bf{82\%}             & \bf{12.22}     & \bf{11.64}
                                      \\ \hline
\end{tabular}
}
\end{center}
\end{table}

\section{Conclusions}
In this paper, we propose an accurate trajectory tracking framework, QuaDRED-SMPC, for quadrotors operating in environments with variable aerodynamic forces. The QuaDRED-SMPC combines aerodynamic disturbance estimation and stochastic optimal control to address the aerodynamic effects on quadrotor tracking. A quantile-approximated distributional RL, QuaDRED, is developed to improve the accuracy of aerodynamic effect estimation, and its convergence is analyzed. Using SADF for control parameterization to guarantee convexity, an SMPC is used to avoid conservative control returns and significantly improves the accuracy of quadrotor tracking. The aerodynamic disturbances are considered to have non-zero mean in the entire QuaDRED-SMPC framework. We demonstrate that our proposed approach can track aggressive trajectories accurately under complex aerodynamic effects whilst guaranteeing both the convergence of QuaDRED and the stability of the whole control framework.
In future works, we will implement QuaDRED-SMPC in real-world flight tests. We will also analyze and evaluate its performance with a greater variety of aerodynamic forces.

\bibliographystyle{unsrt}


\section{Introduction}
The discovery of gravitational-wave radiation has become one of the most important events of recent decades. The existence of gravitational waves (GW) and their direct detection confirmed the consequences of the general relativity (GR) \cite{gw}. Moreover, the data obtained by GW detectors serve to test modified theories of gravity. One way to use GW data for testing the gravitational theories is a comparison of the polarization content predicted by the theory with observational data. The event GW170814 allows to LIGO collaboration to determine that purely tensor polarizations are preferable to scalar and vector ones \cite{gw1}. The data obtained from this event can be used for accurate quantitative test of gravitational theories \cite{liu}. And first of all, it is worth paying attention to the number of polarization modes predicted by the theory.

In GR the gravitational wave has two tensor polarization states, the plus and cross ones. There may exist up to six polarization modes  in alternative metric theories of gravity. For null plane gravitational waves, the six polarizations are classified by the little group E(2) of the Lorentz group developed by Eardley \cite{eardley, eardley1}. According to this classification, the massless scalar field propagates with the transverse breathing mode, and the second scalar mode, the longitudinal one, appeares only in the presence of all the other five modes. However, it was later calculated by Maggiore and Nicolis that in the presence of a massive scalar field, a longitudinal mode appears in addition to the transverse breathing mode, while vector modes are absent \cite{maggiore}. And just recently, the question about the discrepancy between the number of degrees of freedom contained in the theory and the number of polarization modes in the case of massive scalar-tensor theories has been raised \cite{hougong}. This question was considered on the example of Horndeski theory \cite{horndeski}. This model is the most general scalar-tensor theory with second-order field equations. The authors showed that  the polarization state excited by the scalar field form a mix polarization of the
breathing and the longitudinal states and, as a result, it corresponds to one degree of freedom \cite{hougong}.

This work is devoted to the study of the polarization content of hybrid metric-Palatini f(R)-gravity \cite{hybrid1, hybrid2}. The approach of f(R)-gravity allows to extend GR in the simplest way. The f(R)-gravity is based on generalization of the gravitational part of the action as an arbitrary function of the Ricci scalar. Such models have become widespread after f(R)-gravity was successfully applied in inflation theories \cite{starobinsky}. The f(R)-gravity is attractive because the accelerated expansion of the Universe is natural consequence of the gravitational theory. The family of f(R)-theories is divided into two classes: metric one and Palatini one. In the metric approach metric is the only dynamical variable and the action is varied with respect to it only. The Palatini method is based on the idea of considering the connection defining the Riemann curvature tensor to be a priori independent of the metric.
However, some shortcomings of f(R)-theory are manifested both in the metric and in the Palatini approaches. One of the main problems of metric f(R)-gravity is difficulties with passing standard tests in the Solar System \cite{chiba}. Nevertheless, a limited class of viable models in the metric approach exists \cite{odintzov}. The feature of such models is chameleon mechanism. All Palatini f(R)-models lead to microscopic matter instabilities and to unacceptable features in the evolution patterns of cosmological perturbations \cite{koivisto}.

Recently, a theory constructed as a superposition of the Einstein-Hilbert metric Lagrangian with the Palatini f(R)-term has been proposed \cite{hybrid1}. The model was called the hybrid metric-Palatini f(R)-gravity. It unites all positive results of metric and Palatini approaches but lacks their shortcomings. The hybrid f(R)-gravity allows to describe a cosmological large-scale structure, without affecting the Solar System dynamics. An attractive feature of the theory is that it allows non-standard scalar-tensor representation in terms of a dynamical scalar field (unlike the Palatini models), which does not have to be very massive to be consistent with the data obtained from laboratory experiments and Solar System observations. In this theory the scalar field can play an active role in cosmology, without entering the conflict with local experiments \cite{hybrid3}.

The hybrid f(R)-gravity was studied on astrophysical scales from stars to galaxy clusters. It was shown that the virial mass discrepancy in clusters of galaxies can be explained via the geometric terms appearing in the generalized virial theorem \cite{hybrid4}. The hybrid f(R)-gravity also allows to explain the rotational velocities of test particles gravitating around galaxies \cite{hybrid5}. This approach allows to avoid introducing of a huge amount of dark matter. Besides, wormhole \cite{hybrid6} and black hole solutions \cite{bronnikov} were derived, and physical properties of neutron, Bose-Einstein Condensate and quark stars were considered \cite{danila}. In addition, hybrid f(R)-gravity was tested in the Solar System using parametrized post-Newtonian formalism (PPN). The analytical expressions for $\gamma$ and $\beta$ parameters were obtained, and it was proved that other 8 PPN parameters identically equal to zero. It was shown that the light scalar field in hybrid f(R)-gravity does not contradict the experimental data based on all set of PPN parameters \cite{meppn, meshapiro}. Moreover, hybrid f(R)-gravity was tested on the binary pulsars observational data. In addition, the change of the orbital period due to gravitational radiation was obtained in the quasicircular case. Also, for the first time the restriction on the scalar field mass in hybrid f(R)-gravity was found \cite{memnras}. After that, the previous results were generalized to the case of an orbit with eccentricity. The hybrid f(R)-theory was tested using a parameterized post-Keplerian formalism. As a result, it was shown that hybrid f(R)-gravity predicts a wider range of possible masses for the components of binary pulsars in comparison with GR \cite{meppk}.

The hybrid f(R)-gravity can be represented as a scalar-tensor theory with a massive scalar field. Previously, the polarization content has already been studied in the framework of hybrid f(R)-gravity \cite{polarization}. It has been found that the theory contains four polarization states: plus, cross, transverse breathing and longitudinal ones. However, earlier it was shown that the theory has only three degrees of freedom \cite{koivisto1}. Thus, the contradiction that arose earlier in Horndeski theory, exist also in hybrid f(R)-gravity. The purpose of this work is to resolve this inconsistency.

Another effective way to test theories of gravity is to determine the propagation speed of the GW signal predicted by the theory. The supernova, event GW170817/SHB170817A \cite{gw2}, allows to obtain a limit on the propagation velocity of gravitational waves, which led to the reduction or closure of a large number of modified theories of gravity \cite{dead}. This test was also applied to massive scalar-tensor theories. It was shown that only the Brans-Dicke theory with a massive scalar field completely passes this test. In this work, we also calculate the propagation velocity of gravitational waves in hybrid f(R)-gravity and compare obtained result with observational data.

The structure of the paper is the following. In section~\ref{sec2} we consider an action and  field equations of the hybrid metric-Palatini theory in a general form and in a scalar-tensor representation.  In section~\ref{sec3}, we calculate polarization content of hybrid f(R)-gravity. In section~\ref{sec4} we test hybrid f(R)-gravity using the speed of gravitational waves. We conclude in section~\ref{sec:conclusions} with a summary and discussion.

		Throughout this paper the Greek indices $(\mu, \nu,...)$ run over $0, 1, 2, 3$ and the signature is  $(-,+,+,+)$. We will be working in units $h = c = k_B = 1$ throughout the paper.

		\section{Hybrid f(R)-gravity}\label{sec2}

	The action of hybrid f(R)-gravity has the form~\cite{hybrid1,hybrid3}:
	\begin{equation}\label{act}
	S=\frac{1}{2k^2}\int d^4x\sqrt{-g}\left[R+f(\mathfrak{R})\right]+S_m,
	\end{equation}
	where $k^2=8\pi G$, $G$ is the Newtonian gravitational constant, $R$ and $\mathfrak{R}=g^{\mu\nu}\mathfrak{R}_{\mu\nu}$ are the metric and Palatini curvatures respectively, $g$ is the metric determinant, $S_m$ is the matter action. Here the Palatini curvature $\mathfrak{R}$ is defined as a function of $g_{\mu\nu}$ and the independent connection $\hat\Gamma^\alpha_{\mu\nu}$:
	\begin{equation}\label{re}
	\mathfrak{R}=g^{\mu\nu}\mathfrak{R}_{\mu\nu}	=g^{\mu\nu}\bigl(\hat\Gamma^\alpha_{\mu\nu,\alpha}-\hat\Gamma^\alpha_{\mu\alpha,\nu}+\hat\Gamma^\alpha_{\alpha\lambda}\hat\Gamma^\lambda_{\mu\nu}-\hat\Gamma^\alpha_{\mu\lambda}\hat\Gamma^\lambda_{\alpha\nu}\bigr).
	\end{equation}

	Like in the pure metric and Palatini cases, the hybrid f(R)-gravity~(\ref{act}) can be rewritten in a scalar-tensor representation (for details see~\cite{hybrid1,hybrid3}):
	\begin{equation}\label{stact1}
	S=\frac{1}{2k^2}\int d^4x\sqrt{-g}\biggl[(1 + \phi)R + \frac{3}{2\phi}\partial_\mu \phi \partial^\mu \phi - V(\phi)\biggr]+S_m,
	\end{equation}
	where $\phi$ is a scalar field and $V(\phi)$ is a scalar field potential. This potential  $V(\phi)$ is related to the Palatini curvature $\mathfrak{R}$ by the following expression $V_{\phi}=\mathfrak{R}$, where $V_{\phi}$ is the derivative of the potential $V(\phi)$ with respect to the scalar field $\phi$ (for details see~\cite{hybrid1,hybrid3}). Here and further we use the Jordan frame.

	Then the metric and scalar field equations take the following forms~\cite{hybrid1,hybrid3}:
\begin{wide}
	\begin{eqnarray}
 	&&(1+\phi)R_{\mu\nu}=\frac{k^2}{c^4}\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right)-\frac{3}{2\phi}\partial_\mu\phi\partial_\nu\phi +\frac{1}{2}g_{\mu\nu}  \biggl[V(\phi)+\nabla_\alpha\nabla^\alpha\phi  \biggr]+\nabla_\mu\nabla_\nu\phi,\label{feh}\\
 	&&\nabla_\mu\nabla^\mu\phi-\frac{1}{2\phi}\partial_\mu\phi\partial^\mu\phi-\frac{\phi[2V(\phi)-(1+\phi)V_\phi]}{3}=-\frac{k^2}{3c^4}\phi T_{\rm m},\label{fephi}
	 \end{eqnarray}
	 \end{wide}
	where $T_{\mu\nu}$ and $T$ are the energy-momentum tensor and its trace respectively.

	\section{Polarization of gravitational waves}\label{sec3}
	Gravitational radiation propagates by plane waves  in a vacuum. To obtain the linearized field equations in vacuum ($T_{\mu\nu}=0$) we consider the following perturbations of a scalar field and metric tensor:
	\begin{equation}\label{decompos}
	\phi=\phi_0+\varphi,\qquad\ g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},
	\end{equation}
	where $\phi_0$ is the asymptotic background value of the scalar field far away from a source, $\eta_{\mu\nu}$ is the Minkowski background, $h_{\mu\nu}$ and $\varphi$ are small perturbations of tensor and scalar fields respectively. Here we consider $\phi_0$ as a constant. The scalar potential $V(\phi)$ could be expanded in a Taylor series around the background value of the scalar field $\phi_0$ like
	\begin{equation}\label{V}
	V(\phi)=V_0+V'\varphi+\frac{V''\varphi^2}{2!}+\frac{V'''\varphi^3}{3!}...,
	\end{equation}
	hence its derivative with respect to $\varphi$ will be the following: $V_\phi=V'+V''\varphi+V'''\varphi^2/2$.

	The field equation for the scalar field~(\ref{fephi}) in the leading perturbation order takes the form:
	\begin{equation}\label{fephi2}
	\left(\nabla^2-m_\varphi^2\right)\varphi=0,
	\end{equation}
		where we denote $m_\varphi^2=[2V_0-V'-(1+\phi_0)\phi_0V'']/3$ as a scalar field mass.

	The linearized equations for the metric are given by \cite{polarization}
	\begin{equation}\label{feh002}
	\overline{G}_{\mu\nu}=\cfrac{1}{1+\phi_0}(\nabla_\mu\nabla_\nu\varphi-\eta_{\mu\nu}\Box\varphi),
	\end{equation}
where $\overline{G}_{\mu\nu}$ is the linearized Einstein tensor:
\begin{align}\label{gmn}
	\overline{G}_{\mu\nu}&=\overline{R}_{\mu\nu}-\cfrac{1}{2}g_{\mu\nu}\overline{R},\\
\overline{R}_{\mu\nu}&=\cfrac{1}{2}(\partial_\mu\partial_\alpha h^\alpha_\nu+\partial_\nu\partial_\alpha h^\alpha_\mu-\partial_\mu\partial_\nu h-\Box h_{\mu\nu}),\\
\overline{R}&=\partial_\mu\partial_\nu h^{\mu\nu}-\Box h.
	\end{align}
Further we apply the transformations:
\begin{align}\label{theta}
\theta_{\mu\nu}&=h_{\mu\nu}-\cfrac{1}{2}\eta_{\mu\nu}h-\cfrac{1}{1+\phi_0}\eta_{\mu\nu}\varphi,\\
\theta&=-h-4\cfrac{1}{1+\phi_0}\varphi.
\end{align}
The choice of the transverse gauge $\partial_\mu \theta^{\mu\nu}=0$ reduces the field equation (\ref{feh}) to the following form:
\begin{equation}\label{fet}
\Box\theta_{\mu\nu}=0.
\end{equation}

Equations (\ref{fet}) and (\ref{fephi2}) describe the null  gravitational wave in hybrid metric-Palatini f(R)-gravity. The field $\theta_{\mu\nu}$ is responsible for the familiar massless graviton and tensor part contains two polarization states: the plus and cross
ones. The scalar field $\phi$ is massive, and it decouples from the massless tensor field. Suppose the massless and the massive modes both propagate in the +z direction with the wave vectors
\begin{equation}
k^\mu=(\Omega,0,0,\Omega), \ \ \ \ p^{\mu}=(p_t,0,0,p_z),
\end{equation}
respectively, where the dispersion relation for the scalar field is $p_t^2-p_z^2=m^2$. Then the propagation speed of the massive scalar field is $v=\sqrt{p_t^2-m^2}/p_t$. The plane wave solutions of equations (\ref{fet}) and (\ref{fephi2})  take the following forms:
\begin{align}\label{sol}
\theta_{\mu\nu}&=q_{\mu\nu}\exp^{-ik_\alpha x^\alpha},\\
\varphi&=\phi_0\exp^{-ip_\alpha x^\alpha},\ \nonumber
\end{align}
where $\phi_0$ and $q_{\mu\nu}$ are the amplitudes of the waves with $k^\nu q_{\mu\nu}=0$ and $\eta^{\mu\nu}q_{\mu\nu}=0$.

The main purpose of this work is to determine the number of polarization states in hybrid f(R)-gravity. These polarization states are contained in the perturbation tensor $h_{ij}$. For a wave traveling in the z-direction, these polarizations become
\begin{equation}
h_{ij} = \left(
\begin{array}{ccc}
h_b+h_+ & h_\times & h_x\\
 h_\times & h_b-h_+& h_y\\
 h_x & h_y & h_L
 \end{array}
 \right),
\end{equation}
where $h_+$ is plus polarization mode (tensor),  $h_\times$ is cross polarization mode (tensor), $h_b$ is transverse breathing polarization mode (scalar), $h_L$ is longitudual polarization mode (scalar).

 The gravitational waves can influence the distance between the freely moving test particles. This influence is described by the geodesic deviation equation. Assuming that the distance $x^i<<\lambda$, where $\lambda$ is the wavelength of the GWs, and the motion of test particles is slow, the geodesic deviation equation becomes the approximate form $d^2x^i/dt^2 = -R_{0i0j}x^j$, where $R_{0i0j}$ is the electric part of Riemann tensor generated by the GWs. Also The GW field $h_{ij}$ is defined by the Riemann tensor, like $\partial^2 h_{ij}/\partial t^2 = -2R_{0i0j}$. For further calculations it is necessary to use the linear order of the Riemann tensor
\begin{equation}
\overline{R}_{\mu\nu\alpha\beta}=\cfrac{1}{2}(\partial_\mu\partial_\beta h_{\alpha\nu}+\partial_\nu\partial_\alpha h_{\mu\beta}-\partial_\mu\partial_\alpha h_{\nu\beta}-\partial_\nu\partial_\beta h_{\mu\alpha}).
\end{equation}

Calculation of the electric part of Riemann tensor according to the expressions (\ref{sol}, \ref{theta}) gives the following relations
\begin{wide}
\begin{equation}
R_{0i0j} = \left(
\begin{array}{ccc}
-\frac{1}{2(1+\phi_0)}p_t^2\varphi+\frac{1}{2}\Omega^2\theta_{xx}& \frac{1}{2}\Omega^2\theta_{xy} & 0\\
\frac{1}{2}\Omega^2\theta_{xy} & -\frac{1}{2(1+\phi_0)}p_t^2\varphi-\frac{1}{2}\Omega^2\theta_{xx}& 0\\
0 & 0 &-\frac{1}{2(1+\phi_0)}m^2\varphi
 \end{array}
 \right).
\end{equation}
\end{wide}

To study the polarizations caused by the scalar field, one sets $\theta_{\mu\nu}=0$. In the massless case $(m = 0)$, the component of Riemann tensor $R_{tztz}$ is equal to zero, longitudual mode vanishes, and the scalar field excites only the transverse breathing polarization. In the massive case, it is possible to perform a Lorentz boost such that $p_z = 0$, and work in the rest frame of the scalar field. In this rest frame, $p_t^2=m^2$ and $R_{txtx} = R_{tyty} = R_{tztz} \neq 0$. Moreover, the longitudinal and breathing modes coincide in this case. Thus, we obtain that in the chosen reference frame, the longitudinal mode coincides with the breathing one. Due to the absence of a preferred frame, we can conclude that the scalar mode is a mixture of two modes in other  reference frames, which does not contradict the conclusion about the number of degrees of freedom of the theory.

\section{Speed of gravitational waves}\label{sec4}
The GW170817 event gave the possibility to determine the propagation velocity of gravitational waves with unprecedented accuracy \cite{gw2}. As a result, many theories were closed or restricted, including Horndeski model \cite{dead}. Within the framework of this theory, an expression for the propagation velocity of tensor perturbations  was obtained \cite{kase, latosh}. It depends on the general functions of the scalar field included in the structure of Horndeski theory:
\begin{equation}
c_t^2=\cfrac{2G_4-(\dot\phi)^2G_{5,\phi}-(\dot\phi)^2\ddot\phi G_{5,X}}{2G_4-2(\dot\phi)^2G_{4,X}+(\dot\phi)^2G_{5,\phi}-H(\dot\phi)^3G_{5,X}}.
\end{equation}
This expression describes the speed of tensor perturbations $c_t$ propagating over a cosmological background, $X=-1/2\partial_\mu\phi\partial^\mu\phi$, $H$ is Hubble constant.

Hybrid f(R)-gravity is a special case of Horndeski theory. The transition functions were obtained earlier \cite{memnras}:
\begin{equation}\label{hybrid}
	G_2=-\cfrac{3X}{G\phi}-V(\phi),\ G_3=0,\ G_4=\cfrac{1+\phi}{G},\ G_5=0.
	\end{equation}
	Using these expressions, it is possible to obtain the speed of tensor perturbations in hybrid f(R)-gravity:
	\begin{equation}
c_t^2=\cfrac{2\frac{1+\phi}{G}}{2\frac{1+\phi}{G}}=1.
\end{equation}
Thus, gravitational waves in hybrid f(R)-gravity propagate at the speed of light. Thus, the hybrid f(R)-theory is in full agreement existing experimental data.

	\section{Conclusions and discussion}\label{sec:conclusions}

	This work is devoted to the issue of determining the number of polarization states predicted by hybrid f(R)-gravity. This theory belongs to the f(R)-family of theories and unites  metric and Palatini approaches. Hybrid f(R)-gravity can be represented as a massive scalar-tensor theory \cite{hybrid1, hybrid3}. Earlier, the question about the discrepancy between the number of polarization modes in massive scalar-tensor theories and degrees of freedom was raised \cite{hougong}. This issue has been studied in detail on the example of Horndeski theory \cite{hougong}. It was shown that the theory has two tensor  and one scalar degrees of freedom. The scalar degree of freedom is responsible for the mixture of transverse breathing and longitudinal scalar modes. In addition, it was shown that the Eardley classification is not applicable to scalar-tensor theories with a massive scalar field.

For hybrid  f(R)-gravity, the number of polarization states was previously obtained \cite{polarization}. It was concluded, in addition to tensor cross and plus modes, there are transverse breathing and longitudinal modes in the theory. However, it was determined that hybrid f(R)-gravity has three degrees of freedom \cite{koivisto1}. Therefore, the main purpose of this article was to explain the existing discrepancy.

For obtaining the polarization content we calculate the electrical part of the Riemann tensor. Also, we showed that both transverse breathing and longitudinal modes are present. However, in the rest frame of a massive field, it turns out that both modes become identically equal. This implies the conclusion that these two modes are manifestations of the same degree of freedom and represent a mixture of polarization states. Thus, we can conclude there is no contradiction between the number of degrees of freedom and the number of polarization states of the gravitational wave in hybrid f(R)-gravity.

One of the most effective tests of modified theories of gravity  is the calculation of the GWs speed. This test made it possible to cut off a large number of theories, including the scalar-tensor ones. Moreover, some authors argue that only the Brans-Dicke theory with a massive scalar field is able to pass this test \cite{dead}. In this paper, we calculated the propagation velocity of tensor perturbations in hybrid f(R)-gravity and showed that the velocity is identically equal to the speed of light. Thus, hybrid f(R)-gravity is in full agreement with the available experimental data.

This work is the first step in the study of gravitational-wave radiation in the framework of hybrid f(R)-gravity. For a realistic test of the theory, it is necessary to obtain gravitational wave templates and compare them with the available LIGO observational data. Such investigation is supposed to become an important step in development of a new powerful formalism for testing other alternative theories of gravity in the strong field regime of coalescing binary systems.

The author thanks N. A. Avdeev and B.N.Latosh for discussions and comments on the topics of this paper. The work was supported by the Foundation for the Advancement of Theoretical Physics and Mathematics ``BASIS''.