# Distributed Computing for Scalable Optimal Power Flow in Large Radial
Electric Power Distribution Systems with Distributed Energy Resources

Rabayet Sadnan∗ and Anamika Dubey
Department of Electrical Engineering and Computer Science, Washington State
University
Email: {rabayet.sadnan∗<EMAIL_ADDRESS>
###### Abstract

Solving the non-convex optimal power flow (OPF) problem for large-scale power
distribution systems is computationally expensive. An alternative is to solve
the relaxed convex problem or linear approximated problem, but these methods
lead to sub-optimal or power flow infeasible solutions. In this paper, we
propose a fast method to solve the OPF problem using distributed computing
algorithms combined with a decomposition technique. The full network-level OPF
problem is decomposed into multiple smaller sub-problems defined for each
decomposed area or node that can be easily solved using off-the-shelf
nonlinear programming (NLP) solvers. Distributed computing approach is
proposed via which sub-problems achieve consensus and converge to network-
level optimal solutions. The novelty lies in leveraging the nature of power
flow equations in radial network topologies to design effective decomposition
techniques that reduce the number of iterations required to achieve consensus
by an order of magnitude.

###### Index Terms:

Distributed Computing, Distributed Optimization, Optimal Power Flow, Power
Distribution Systems.

##  I Introduction

The nature and the requirements of the power systems, specially at the
distribution-level is changing rapidly with a large-scale integration of
controllable distributed energy resources (DERs). The continued proliferation
of DERs, that include Photovoltaic (PV) systems, battery energy storage units
(BESS), and controllable loads such as Electric Vehicles (EVs) is leading to a
drastic increase in the number of active nodes at the distribution-level that
need to be controlled/managed optimally for an efficient and resilient grid
operations. Traditionally, grid operations are centrally managed upon solving
an optimal power flow (OPF) problem where centralized optimization techniques
are used to solve the resulting difficult non-linear non-convex OPF problem
[1, 2]. Unfortunately, the computational challenges, primarily posed by the
non-convex power flow constraints in OPF formulation, increases drastically
with the size of the distribution systems motivating computationally efficient
approaches [3].

Existing methods manage the computational challenges using convex relaxation
or linear approximation techniques [4, 5]. The primary drawbacks of the convex
relaxed models are the possibilities of inexact/infeasible power flow
solutions [6]. The linear approximated models may lead to NLP infeasible
solutions and high optimality gap depending upon the problem type. Moreover,
methods based on both approximation and relaxation techniques use a
centralized paradigm that may lead to scalability challenges as the problem
size increases. With a majority of DER integration happening at the secondary
feeder level, the OPF problem will need to solve even larger feeder with
thousands of secondaries. For example, the largest IEEE test feeder is a 8500
node test system that terminates at the secondary transformer level and does
not include secondary feeders. If each service transformer is expanded to a 20
node secondary feeder, it will lead to a total of 22000 secondary nodes added
to the problem formulation. Such problem complexities motivate the move
towards a distributed computing paradigm. In general, the distributed
computing is facilitated by the rapid growth in high-performance computing
using many-core machines. Fortunately, the radial operational topology of
power distribution systems make it highly conducive for parallelization and
distributed computing, In this paper, we develop a distributed computing
approach for OPF problems for distribution systems, that can scale for very
large distribution feeders and converges using fewer message passing among
distributed computing nodes thus significantly reducing the overall compute
time.

Within this context, existing literature includes numerous approaches on the
application of distributed optimization algorithms for power distribution
systems [7, 8]. In general, these methods adopt the traditional distributed
optimization techniques to model a distributed optimal power flow (D-OPF)
problem [7, 9, 10, 8]. A D-OPF formulation decomposes the OPF into several
smaller subproblem that require multiple micro- and macro-iterations for
convergence. Micro-iterations involve solving the distributed sub-problems in
parallel. And macro-iterations involve exchanging the solutions or more
specifically the updated boundary variables obtained from the distributed
subproblems. Both micro and macro-iterations together decide the time-of-
convergence (ToC) for the algorithm. Unfortunately, the exiting distributed
optimization algorithms require a very large number of macro-iterations to
converge for medium-scale distribution grids [11, 12, 8, 13]. A practical
implementation of such algorithms requires a very fast communication among
distributed computing nodes to reach a converged solution within a reasonable
time. A large number of communication rounds/message-passing events among
distributed agents is not preferred since this leads to significant delays in
decision-making. Lately, to address some of these challenges, real-time
feedback based online distributed algorithms have been explored in the related
literature for network optimization [14, 15, 16, 17, 18, 19]. Generally, these
algorithms do not wait to optimize for a time-step but asymptotically arrive
at an optimal decision over several steps of real-time decision-making.
However, these algorithms also take hundreds of iterations to converge/track
the optimal solution for a mid-size feeder. This raises further challenges to
the performance of the algorithm for larger feeders, especially during the
fast-varying phenomenon and slow communication channels.

To address these challenges, recently we have developed a distributed OPF
formulation for the radial distribution systems based on the equivalence of
networks principle [20, 21]. In this paper, we further expand the previously
proposed method to solve OPF for very large notional distribution test feeders
(with 10,000 nodes) for several different problem objectives. The proposed
approach solves the original non-convex optimal power flow problem for power
distribution systems using a novel decomposition technique combined with
distributed computing approach. The distributed subproblems are related via
the boundary variables shared by the neighboring nodes. First, the low-compute
distributed OPF sub-problems are locally solved. The consensus of the boundary
variables is achieved using a Fixed-Point Iteration (FPI) algorithm. Upon
consensus, the solutions converge to network-level OPF solutions. The proposed
approach leverages the radial topology of the power distribution system and
the associated unique power flow properties to develop message passing
routines that reduces the number of message-passing among distributed agents
by an order of magnitude. We demonstrate the proposed approach for three
problem objectives (1) loss minimization (2) DER generation maximization and
(3) voltage deviation minimization using a single-phase equivalent of
8500-node test feeder (with 2500 nodes) and a balanced synthetic 10,000 node
distribution feeder. The proposed approach is shown to scale for all problem
objectives while most centralized formulation can’t be solved for more than
2000 nodes using off-the-shelf optimization solvers such as Artelys Knitro. To
our knowledge, this is the first paper to demonstrate an approach that solves
such a large-scale D-OPF on a regular CPU without the use of any high-
performance computing (HPC) machines.

##  II Centralized OPF Model

In this paper, $(\cdot)^{*}$ represents the complex-conjugate; $(\cdot)^{T}$
represents matrix transpose; $|~{}.~{}|$ symbolizes the absolute value of a
number or the cardinality for a discrete set; $(\cdot)^{n}$ represents the
$n^{th}$ iteration; and $j=\sqrt{-1}$ in a complex number; $\overline{(.)}$
and $\underline{(.)}$ denotes the maximum and minimum limit of a given
quantity.

###  II-A Network and DER models

Let us consider a balanced radial power distribution network – represented by
the directed graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$, where
$\mathcal{N}$ be the set of all nodes in the system and $\mathcal{E}$ denotes
the set of all distribution lines connecting the ordered pair of buses $(i,j)$
i.e., from node $i$ to node $j$. Also, $r_{ij}+jx_{ij}$ is the series
impedance $\forall\\{ij\\}\in\mathcal{E}$. Let, for node $j$, $k$ be the set
of all children nodes; thus, in $k:j\rightarrow k$, $k$ represents the set of
children nodes for the node $j$. Next, we denote
$v_{j}=|V_{j}|^{2}=V_{j}{V_{j}}^{*}$ as the squared magnitude of voltage at
node $j$. Let $l_{ij}$ be the squared magnitude of current flow in branch
$\\{ij\\}$. We denote $P_{ij},Q_{ij}$ as the sending-end active and reactive
power flows for branch ${ij}$, and complex power
$S_{L_{j}}=p_{L_{j}}+jq_{L_{j}}$ is the load connected and
$S_{D_{j}}=p_{Dj}+jq_{Dj}$ is the power output of DER connected at node $j$.
The network is modeled using the branch flow equations [22] defined for each
line $\\{ij\\}\in\mathcal{E}$ and $\forall j\in\mathcal{N}$ in (1a).

$\displaystyle\small P_{ij}-r_{ij}l_{ij}-p_{L_{j}}+p_{Dj}=\sum_{k:j\rightarrow
k}P_{jk}$ (1a) $\displaystyle
Q_{ij}-x_{ij}l_{ij}-q_{L_{j}}+q_{Dj}=\sum_{k:j\rightarrow k}Q_{jk}$ (1b)
$\displaystyle
v_{j}=v_{i}-2(r_{ij}P_{ij}+x_{ij}Q_{ij})+(r_{ij}^{2}+x_{ij}^{2})l_{ij}$ (1c)
$\displaystyle v_{i}l_{ij}=P_{ij}^{2}+Q_{ij}^{2}$ (1d)

The DERs are modeled as Photovoltaic modules (PVs) interfaced using smart
inverters, capable of two-quadrant operation. If the reactive power
generation, $q_{Dj}$, is controllable and modeled as the decision variable for
the optimal operation, then the real power generation by the DER, $p_{Dj}$, is
assumed to be known (measured). Let the rating of the DER connected at node
$j$ be $S_{DRj}$, then the limits on $q_{Dj}$ are given by (2).

$-\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}\leq q_{Dj}\leq\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}$
(2)

On the contrary, if the active power generation, $p_{Dj}$, is modeled as the
decision variable, then $q_{Dj}$ is set to $0$, and $p_{Dj}$ can vary between
$0$ and $S_{DRj}$, see (3).

$0\leq p_{Dj}\leq S_{DRj}$ (3)

###  II-B Centralized OPF problems

To optimize the network for some cost function, we define a centralized OPF
problem defined by a network-level problem objective, the power flow models in
(1a), and the operating constraints on the power flow variables. In this
paper, we formulate three different optimal power flow problems for the power
distribution grids, (i) active power loss minimization, (ii) DER generation
maximization, and (iii) Voltage deviation ($\Delta$V) minimization. The
corresponding OPF problems are detailed below.

####  II-B1 Loss Minimization

The problem objective is to reduce the network losses by controlling the
reactive power output from DERs ($q_{Dj}$). Let
$X_{lm}=[P_{ij},~{}Q_{ij},~{}l_{ij},~{}v_{j},~{}q_{Dj}]^{T}$ be the problem
variables $\forall j\in\mathcal{N}$, and $\forall\\{ij\\}\in\mathcal{E}$. Note
that, if node $j$ doesn’t have any DER, then $q_{D_{j}}=0$. Also, let
$F_{lm}(X_{lm})$ denote the objective function representing the total power
loss in the given distribution system. Note that $F_{lm}(X_{lm})$ is a
function of both the power flow variables and decision variables. Then, the
OPF problem is defined as the following in (C1).

$\displaystyle\text{{(C1)}}\hskip 11.38092pt\min\hskip
5.69046ptF_{lm}(X_{lm})=\sum_{\\{ij\\}\in\mathcal{E}}l_{ij}r_{ij}\hskip
17.07182pt$ (4a) s.t.(1a) and (2) (4b) $\displaystyle\underline{V}^{2}\leq
v_{j}\leq\overline{V}^{2}\hskip 19.91684pt;\forall j\in\mathcal{N}$ (4c)
$\displaystyle l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}\hskip
14.22636pt;\forall\\{ij\\}\in\mathcal{E}$ (4d)

where, $\overline{V}=1.05$ and $\underline{V}=0.95$ are the limits on bus
voltages, and $(I^{rated}_{ij})^{2}$ is the thermal limit for the branch
$\\{ij\\}$.

####  II-B2 DER Maximization

In DER maximization problem objective, the DER active power generation is
maximized without violating the operational limits of the distribution system.
This is achieved by maximizing the active power output from DERs ($p_{Dj}$).
Let $X_{dm}=[P_{ij},~{}Q_{ij},~{}l_{ij},~{}v_{j},~{}p_{Dj}]^{T}$ be the
problem variables. Here, the objective function is denoted by
$F_{dm}(X_{dm})$, representing the total active power generation by DERs.
Then, this DER maximization OPF problem is defined as the following in (C2).
Similar to the previous formulation, if any node $j$ doesn’t have any DER,
then we set $p_{D_{j}}=0$.

$\displaystyle\text{{(C2)}}\hskip 11.38092pt\max\hskip
5.69046ptF_{dm}(X_{dm})=\sum_{j\in\mathcal{N}}p_{Dj}\hskip 17.07182pt$ (5a)
s.t.(1a) and (3) (5b) $\displaystyle\underline{V}^{2}\leq
v_{j}\leq\overline{V}^{2}\hskip 19.91684pt;\forall j\in\mathcal{N}$ (5c)
$\displaystyle l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}\hskip
14.22636pt;\forall\\{ij\\}\in\mathcal{E}$ (5d)

Kindly note that DER maximization problem is also known as PV hosting problem
if DERs are modeled as PV modules.

####  II-B3 $\Delta$V Minimization

In this specific network level optimization problem, we try to keep the nodal
voltages as close as possible to a pre-specified reference, $V_{ref}$. The
problem objective is to minimize the nodal voltage deviations from the
reference value by controlling the reactive power output from DERs ($q_{Dj}$).
The problem variables are denoted by
$X_{dv}=[P_{ij},~{}Q_{ij},~{}l_{ij},~{}v_{j},~{}q_{Dj}]^{T}$, $\forall
j\in\mathcal{N}$ and $\forall\\{ij\\}\in\mathcal{E}$. Also, the cost function,
$F_{dv}(X_{dv})$, represents the total two-norm distances of nodal voltages,
$v_{j}$, from reference voltage $v_{ref}$. Mathematically
$F_{dv}(X_{dv})=\sqrt{\sum(v_{j}-v_{ref})^{2}}$, $\forall j\in\mathcal{N}$.
The OPF problem is defined as the following in (C3). Here in this paper, we
used $V_{ref}=1.00$ as the bus reference voltage.

$\displaystyle\text{{(C3)}}\hskip 11.38092pt\min\hskip
5.69046ptF_{dv}(X_{dv})=\sqrt{\sum_{\forall
j\in\mathcal{N}}(v_{j}-v_{ref})^{2}}\hskip 17.07182pt$ (6a) s.t.(1a) and (2)
(6b) $\displaystyle\underline{V}^{2}\leq v_{j}\leq\overline{V}^{2}\hskip
19.91684pt;\forall j\in\mathcal{N}$ (6c) $\displaystyle
l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}\hskip
14.22636pt;\forall\\{ij\\}\in\mathcal{E}$ (6d)

Assumption 1: The loads in the network for all three OPFs are modeled as
constant power loads; i.e., in ZIP load model, $(Z,I,P)=(0,0,1)$.

In the next section, we detail the method on how to decompose the optimization
problems for large scale distribution grids into several sub-problems, solve
in parallel, and converge into the final solution.

##  III Decomposition of the Central OPF Problem

The OPF problems described in the previous section are formulated as a
centralized optimization problem for the radial power distribution systems.
For a large scale distribution system with thousands of nodes and decision
variables, solving the NLP OPF is computationally expensive and difficult to
converge for very large-scale distribution systems. Since the power
distribution system is operated radially, the OPF problems defined in
(C1)-(C3) are naturally decomposable into multiple sub-problems defined for
the connected areas. The details of the proposed problem decomposition
technique and the resulting distributed OPF problem are discussed next.

Figure 1: An example of a two-area system.

###  III-A Decomposition of the OPF Problem

First, we decompose the whole distribution grid into $N$ smaller areas. Let
$\mathcal{A}_{R}=\\{A_{1},~{}A_{2},…,~{}A_{N}\\}$, be the set of all
decomposed areas. Also, let each area, $A_{m}\in\mathcal{A}_{R}$, be defined
as a directed graph $A_{m}=\mathcal{G}(\mathcal{N}_{m},\mathcal{E}_{m})$.
Here, each area $A_{m}$ has a maximum number of nodes/variables, so that the
respective OPF sub-problems for that area can easily be solved using off-the-
shelf NLP solvers. The coupling/complicating variables among these smaller
sub-problems, associated with respective areas, are directed by the structure
of the network. Since, the grid was radial to begin with, the decomposed
areas, or the sub-trees of the networks are also connected radially with each
other. This specific structure of the network helps to identify the unique
parent area and the child areas for any area $A_{m}$, which in turns
associates the complicating/shared variables – exchanged among sub-problems to
solve the overall master problem. For this decomposition method, the
complicating variables are the shared bus voltages and power flows in the
shared bus. Computationally, sub-problems associated with each area is solved
in parallel by assuming a fixed voltage at the shared bus with the unique
parent area, and a constant loads at the shared buses with child areas. After
solving the sub-problems, the respective complicating variables, i.e., the
total power requirements in that area is shared with sub-problem for the
parent area and the shared bus voltages are shared with sub-problems
associated with child areas. This exchange of complicating variables is called
as macro-iteration here. After this macro-iteration, sub-problems are solved
again, till convergence for the complicated variables are achieved.

Specifically, lets assume the network is decomposed into $2$ areas – area
$A_{1}$ and $A_{2}$; each with their own, purely local variables – $x_{1}$ and
$x_{2}$ (see Fig. 1). Area $A_{1}$ is the parent area of area $A_{2}$. Let
$\textbf{Y}=[y_{1},y_{2}]^{T}$ be the complicating variable that couples the
sub-problems for the two areas. Here, $y_{1}$ and $y_{2}$ are the bus voltage
magnitude ($v$) and the complex power flow through the bus ($S=P+jQ$) shared
between $A_{1}$ and $A_{2}$, respectively; i.e.,
$[y_{1},y_{2}]^{T}=[v,S]^{T}$. If the set of all local variables for $A_{1}$
and $A_{2}$ is denoted by $X_{1}$ and $X_{2}$, respectively, then
$X_{1}=\\{x_{1},y_{1}\\}{\hskip 5.69046pt}\text{\&{\hskip
5.69046pt}}X_{2}=\\{x_{2},y_{2}\\}$. Let $X=X_{1}\cup X_{2}$ be the set of
central OPF problem variables, and $\mathcal{S}$ is the set of constraints for
the overall centralized optimization problem. If $F$ is a decomposable cost
function, then the problem can be decomposed and written as (7), where,
$\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are the set of constraints on local
variables for decomposed area $A_{1}$ and $A_{2}$, respectively. Also,
$f_{1}$, $f_{2}$ are the cost functions for the respective local sub-problems.

$\small\min_{X\in\mathcal{S}}{\hskip
2.84544ptF(X)}=\min_{X_{1}\in\mathcal{S}_{1},X_{2}\in\mathcal{S}_{2}}{\hskip
2.84544pt}f_{1}(X_{1},y_{2})+f_{2}(X_{2},y_{1})$ (7)

The original problem defined in (7) can be readily decomposable into the
following two sub-problems (see (8)), associated for respective decomposed
areas; i.e., equation (8a) and (8b) for area $A_{1}$ and $A_{2}$,
respectively.

$\displaystyle\text{For $A_{1}$}\hskip 8.5359pt:\min_{X_{1}\in S_{1}}{\hskip
5.69046pt}f_{1}(X_{1},y_{2})$ (8a) $\displaystyle\text{For $A_{2}$}\hskip
8.5359pt:\min_{X_{2}\in S_{2}}{\hskip 5.69046ptf_{2}(X_{2},y_{1})}$ (8b)

Remark 1: Please note that, the decomposition of the OPFs also works for any
maximization problem, such as (C2).

Remark 2: The decomposition method described here can easily be extended for a
network, where multiple area decomposition is required to make each sub-
problem small enough to make it solvable efficiently. Similar to the 2-area
distributed OPF, the optimization problem can be decomposed into several
smaller sub-problems, representing one decomposed area of the network.

Remark 3: The decomposition approach can be further extended to nodal
decompositions, where each node represents one area.

###  III-B Consensus for the Decomposed Sub-problems

After decomposing the optimization problem into several smaller sub-problems,
the proposed distributed algorithm solves the sub-problems individually to
obtain respective local and complicating variables. Here, at each boundary
among decomposed areas, the complicating variable $y_{2}$ and $y_{1}$ are kept
fixed to solve sub-problem (8a) and sub-problem (8b), respectively. Then the
solved $y_{1}$ by sub-problem (8a) and solved $y_{2}$ by sub-problem (8b) are
exchanged again between areas. After each macro-iteration, the update step of
complicated variable, Y, is performed using Fixed Point Iteration (FPI)
method, described by (9) for $n^{th}$ macro-iteration. Here, instead of a
constant value, $\alpha$ can be made adaptive as well. The exchange process is
repeated until the change in the all complicating variables for all decomposed
boundaries over macro-iterations are within tolerance $\epsilon_{tol}$, see
(10).

$\textbf{Y}^{(n)}:=\frac{\textbf{Y}^{(n)}+\alpha\textbf{Y}^{(n-1)}}{1+\alpha}$
(9)

$\Big{|}\textbf{Y}^{(n)}-\textbf{Y}^{(n-1)}\Big{|}\leq\epsilon_{tol}$ (10)

In the next section, we will discuss the distributed approach to solve large
scale OPF problems for radial power distribution networks using the
decomposition method described in this section.

##  IV Distributed OPF for Scalability

In this section, we describe the distributed method to solve large scale OPF
problems. For that, we use the decomposition technique, that we proposed
previously. Specifically, we will detail the distributed algorithm to solve
previously developed central OPFs, i.e., (C1), (C2), and (C3). First, we
discuss the formulation of the sub-problems, and then we describe the
algorithm.

###  IV-A Distributed Sub-Problems

For a system decomposed into multiple areas, the sub-problems are defined for
each $A_{m}\in\mathcal{A}_{R}$, where
$A_{m}=\mathcal{G}(\mathcal{N}_{m},\mathcal{E}_{m})$. Please note, while
decomposing the network, it is ensured that the number of nodes/variables in
each decomposed areas do not exceeds a certain number, that might cause
computation complexities. Here, the power flow model is defined in
(11b)-(11d), and are used by the corresponding sub-problem for area $A_{m}$ –
defined $\forall j\in\mathcal{N}_{m}$ and
$\forall\\{ij\\},\\{jk\\}\in\mathcal{E}_{m}$. Also, we let $C_{h}$ be the set
of buses, that is shared by area $A_{m}$ with its child areas. The sub-
problems for (i) loss minimization, (ii) DER maximization, and (iii) $\Delta$V
minimization is detailed next. For these OPF objectives, we use the same OPF
formulation as central problem, except only for the corresponding area,
$A_{m}$. Also, the shared bus voltage and power flow (complicated variables)
are updated using (9), and kept fixed for $n^{th}$ iteration, as shown in
equation (11e)-(11g). The sub-problem for loss minimization is described
below, (11a).

$\displaystyle\text{{(D1)}}\hskip 11.38092pt\min\hskip
5.69046ptf_{m}=\sum_{\\{ij\\}\in\mathcal{E}_{m}}l_{ij}r_{ij}\hskip 17.07182pt$
(11a) $\displaystyle\text{s.t.}\hskip
5.69046ptP_{ij}-r_{ij}l_{ij}-p_{L_{j}}+p_{Dj}=\sum_{k:j\rightarrow k}P_{jk}$
(11b) $\displaystyle Q_{ij}-x_{ij}l_{ij}-q_{L_{j}}+q_{Dj}=\sum_{k:j\rightarrow
k}Q_{jk}$ (11c) $\displaystyle
v_{j}=v_{i}-2(r_{ij}P_{ij}+x_{ij}Q_{ij})+(r_{ij}^{2}+x_{ij}^{2})l_{ij}$ (11d)
$\displaystyle\vspace{1pt}v_{o}=v_{o^{\prime}};$ (11e)
$\displaystyle\vspace{1pt}P_{jk}=p_{k^{\prime}}\hskip
5.69046pt;\forall\\{jk\\}\in\mathcal{E}_{m}\text{, where }k\in C_{h}$ (11f)
$\displaystyle\vspace{1pt}Q_{jk}=q_{k^{\prime}}\hskip
5.69046pt;\forall\\{jk\\}\in\mathcal{E}_{m}\text{, where }k\in C_{h}$ (11g)
$\displaystyle\vspace{2pt}-\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}\leq
q_{Dj}\leq\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}$ (11h)
$\displaystyle\underline{V}^{2}\leq v_{j}\leq\overline{V}^{2}\hskip
19.91684pt;\forall j\in\mathcal{N}_{m}$ (11i) $\displaystyle
l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}\hskip
14.22636pt;\forall\\{ij\\}\in\mathcal{E}_{m}$ (11j)

Here, area $A_{m}$ shares bus ’$o$’ with its parent area, and $v_{o^{\prime}}$
is the solved bus voltage by that parent area in the previous iteration.
Similarly, $\forall k\in C_{h}$, $,~{}p_{k^{\prime}},\&~{}q_{k^{\prime}}$ is
the solved shared power flows by the child areas of $A_{m}$ in the previous
iteration. Note that, the symbol $|.|^{\prime}$ depicts that the variable is
solved by other areas, and not by the area that is associated with the
corresponding sub-problem. Further, the sub-problems for DER maximization and
$\Delta$V minimization is formulated in (12a) and (13a), respectively.

$\displaystyle\text{{(D2)}}\hskip 11.38092pt\max\hskip
5.69046ptf_{m}=\sum_{j\in\mathcal{N}_{m}}p_{D_{j}}\hskip 17.07182pt$ (12a)
s.t.(11b) - (11g), (3), (11i) - (11j) (12b) $\displaystyle\text{{(D3)}}\hskip
11.38092pt\min\hskip 5.69046ptf_{m}=\sqrt{\sum_{\forall
j\in\mathcal{N}_{m}}(v_{j}-v_{ref})^{2}}\hskip 17.07182pt$ (13a) s.t.(11b) -
(11g), (2), (11i) - (11j) (13b)

###  IV-B Algorithm

For completeness, now we discuss the full distributed algorithm that
decomposes the OPFs for large scale power distribution systems, and solves
iteratively to reach the global solution. Here, we use the decomposition
technique that we developed in Section III, and solve sub-problems for
different network level objectives, i.e., (D1)-(D3) until convergence. We use
the same notation for network variables that we denoted for decomposed areas
in Section IIIA and IIIB. We use tolerance of $\epsilon_{tol}=0.001$ to meet
the convergence criterion. The algorithm is detailed below (see Algorithm 1).
To better understand the distributed computing of the OPFs, the sub-routine in
step 5 of Algorithm 1 is described in Algorithm 2,

1 Decompose the network into $N$ areas, so that each area has a maximum
specified node numbers

2 Initialize complicating variables, ${\textbf{Y}}^{0}\in\mathcal{S}$; error,
$e=1$; and macro-iteration count $n=0$

3 If $|e|\leq\epsilon_{tol}$, stop iteration count, and go to step 10

4 Else, increase iteration count $n$: $n\leftarrow n+1$

5 Solve $\Phi_{m}$ in parallel using Algorithm 2, for all decomposed areas,
where, $\Phi_{m}$ depicts the solution of sub-problems; $\Phi_{m}:\hskip
2.84544ptX_{m}^{(n)}:=\underset{X_{m}\in S_{m}}{\text{argmin/argmax}}\hskip
5.69046ptf_{m}\Big{(}X_{m},{y}_{m^{^{\prime}}}^{(n-1)}\Big{)}$

6 Update all the complicating variables, Y, using (9), where $\alpha$ can be
constant or adaptive

7 Check residual vector
$\mathcal{R}^{(n)}=\left[\begin{array}[]{c}\textbf{Y}^{(n)}-\textbf{Y}^{(n-1)}\end{array}\right]$

8 $e=\max\hskip 3.0pt|\mathcal{R}^{(n)}|$

9 Go to step 2

Return Global Minimizer: $X^{*}=\\{X_{m}^{n}~{}|~{}m=1,2,...,N\\}$

Algorithm 1 Distributed Algorithm for Scaled OPFs

Sub-Problem : For decomposed area $A_{m}\in\mathcal{A}_{R}$

Macro-Iteration step : n

Complicated Variables : $y_{m^{^{\prime}}}^{(n-1)}$, variables that are used
for coupling sub-problems

Optimization Variable : $X_{m}$

Steps :

1 Assume the complicating variable for area $A_{m}$ as constant; i.e.,
$y_{m^{^{\prime}}}^{(n-1)}$ is set to either constant voltage (if it has
parent area), or constant loads (if it has child areas), or both – depending
on the position of the area $A_{m}$ (See equation (11e)-(11g))

2 Solve distributed sub-problems of minimizing or maximizing the decomposed
cost function $f_{m}$, e.g., (D1), (D2), etc. by assuming the constant
complicating variables by off-the-shelf NLP solvers

Store the local minimizer in the variable $X_{m}^{(n)}$

Algorithm 2 Sub-routine to Solve Sub-Problems at Step-5

##  V Numerical Simulation

To show the efficacy and validate the proposed decomposition approach to solve
the large scale OPFs for power distribution systems, we simulate our algorithm
for a very large scaled, balanced synthetic 10,000 node distribution system
and medium-scale balanced IEEE-8500 node test system with 2500 nodes. All
experiments were simulated in Matlab 2018b on a machine with 8GB memory and
Core i7-8700 CPU @3.19 GHz. The NLP sub-problems of the distributed method is
solved using fmincon of Matlab using ’sqp’ algorithm. However, given the NP-
hard nature of the centralized OPF problems, for bench-marking against
centralized OPF, we use a commercial NLP solver Artelys Knitro with ’active-
set’ algorithm, that scales relatively well with the problem size [23].

###  V-A Simulated System

The simulations are conducted using the following two test systems: (1)
Synthetic 10,000 node distribution system with different DER penetration
levels, and (2) Balanced IEEE-8500 node test system with 100% DER penetration
for nodal decomposition. Please note, the % penetration means the percentage
of DER nodes compare to load nodes. The Synthetic 10,000 node system is shown
in Fig. 2. The distribution system is comprised of 1 main-feeder, 20 laterals,
where each lateral supplies 20 neighborhoods. It is assumed that each
neighborhood is comprised of 20 households. Thus, each lateral supplies a
total of 400 houses. Also, in between 2 laterals, we assume 4 nodes in the
main-feeder that represent the distributed loads. Every load in this
distribution network is set to consume a total of $S_{L}=0.1+0.01j$ pu, and
the line impedance of all the branches is assumed to be $z=0.07+0.01j$ pu. The
base voltage for the network is 12.47 kV ($V_{LL}$) and base kVA is 1000. For
loss minimization and $\Delta$V minimization objectives, each DER in the
system can generate 7 kW of real power, with nominal rating of 8.4 kVA. For
DER maximization problem, the rating of the DERs are increased to at-most 10
times to stress the system. We decompose the distribution system in multiple
areas where each area is composed of 100 nodes (see Fig. 2).

For the IEEE 8500 node test system, the DER sizes are chosen randomly with a
rating ranging from 1.3 to 5.8 kVA. We use this medium-scale distribution
system to further decompose the problem into nodal level, i.e., each node is
considered as an area. Please note that, this is a balanced, single-phase
equivalent distribution network of the test system, that has 2522 nodes. We
used the same base values, i.e., 12.47 kV as base voltage and 1000 kVA as base
kVA for this test system. The proposed decomposition technique is next
simulated for various DER penetration with different network objectives to
solve the OPFs for these scaled-networks.

Figure 2: Synthetic 10,000 Node System

###  V-B Loss Minimization Objective: (D1)

First, we solve the loss minimization problem (D1) for 10,000 node test system
with varying DER penetration levels. The active power loss in the network is
minimized by the reactive power generations of the DERs. For loss minimization
OPFs, we have used the nominal case for synthetic 10,000 node system. That is
all the loads are at nominal value ($S_{L}$). The kVA rating of these DERs are
120% of the nominal active power generation. We have simulated (i) 100%, (ii)
50%, and (iii) 10% DER penetration levels for loss minimization objective with
grid voltage of 1.00 pu. The result of this OPF is detailed next. We have used
$\alpha=0$ for FPI update in (9).

The converged solution of the decomposed central OPF and the convergence
properties of the proposed method for the loss minimization problem is shown
in Fig. 3. We can see that the converged voltage does not violate any voltage
constraints, i.e., the voltage is not outside of the specified limit of
0.95-1.05 pu bound (Fig. 3(a)). Also, with increased penetration, the nodal
voltages over the network has less standard deviation (S.D.), that results in
lower active power loss in the network. The nodal voltages has higher S.D. for
10% DER penetration, i.e., more spreaded for less DER penetration. Fig. 3(b)
shows the objective value of the OPF problem w.r.t. macro-iterations. 100%
penetration can reduce the line losses to 4.5 kW, but with lower DER
penetrations, the line losses increases. The convergence properties for this
case is shown in Fig. 3(c). For all the cases, it takes around 11 macro-
iterations to meet the convergence criterion. Besides, the time taken at each
iteration for this case is plotted in Fig. 3(d). This time represents the
highest time taken to solve any sub-problem at each iteration. It only takes
$\sim 30$ seconds to solve the OPF by decomposing the problems into several
sub-problems for all the DER penetration levels (see Table I).

(a) Nodal voltages

(b) Objective value

(c) Convergence at the boundary

(d) Simulation time

Figure 3: Numerical Results for Loss Minimization Objective for Synthetic
10,000 Node System

(a) Nodal voltages

(b) Objective value

(c) Convergence at the boundary

(d) Simulation time

Figure 4: Numerical Results for DER Maximization Objective for Synthetic
10,000 Node System

###  V-C DER Maximization Objective: (D2)

In this section, we present the result for DER maximization OPF problem, which
is also equivalent to DER curtailment problem for the power distribution
networks. Here, similar to the previous OPF problem, we solve the decomposed
problem (D2) for 10,000 node test system with different DER penetration
levels. The active power generation of the DERs are maximized while
maintaining the operation limits of the network, such as voltage limits. For
this optimization problem, we have used various load and generation multiplier
to stress the system at max level. We simulate 3 different cases – (i) 50% DER
penetrations where the active power generation capacities of DERs are 21 kW
and loads are nominal, (ii) 20% DER penetrations with 28 kW of active power
generation capacities for each DERs and load multiplier is 0.5, and (iii) 10%
DER penetration levels with max of 70 kW generation capabilities and load
multiplier is 0.5. The grid is assumed to be operating at 1.05 pu. The result
of this OPF is discussed next. We have used $\alpha=2.33$ for FPI update in
(9).

From the Fig. 4, we can see the converged solution of DER maximization OPF,
that has been solved using proposed decomposition method. Similar to the
previous objective, we can see that the nodal voltage does not violate the
voltage limits (Fig. 4(a)). The voltages are near the upper bound shows that
the systems were highly stressed for these different simulated cases. With
increased penetrations, more nodes have voltages that are closer to the upper
limits of 1.05 pu. Fig. 4(b) shows the normalized objective value of the OPF
problem w.r.t. macro-iterations. Here, the objective value is scaled w.r.t.
the converged/final cost as the orders of the final costs are different. The
actual values of the objective function upon solving OPFs using distributed
approach is shown in table I. Note that, even though the number of DERs in 10%
penetration case is lower than 20% case, but individual DERs has higher
capacity in 10% penetration case than later, and thus total generation is
higher in 10% penetration case than 20% penetration case. It takes 7, 10 and
14 iterations to converge for 50, 20 and 10% DER penetrations (Fig. 4(b),4(c)
). The simulation time per macro-iteration is shown in Fig. 4(d). The total
simulation time to solve the OPF using proposed decomposition approach is
reported in Table I.

TABLE I: Results Summary OPF Problem | DER% | Converged Objective Value | Time (s)
---|---|---|---
Loss Min | 100 | 4.5 kW | 34
50 | 21.58 kW | 35
10 | 44.05 kW | 30
DER Max | 50 | 1.03 MW (Total capacity 1.05 MW) | 120
20 | 0.55 MW (Total capacity 0.56 MW) | 240
10 | 0.66 MW (Total capacity 0.70 MW) | 300
$\Delta$V Min | 100 | 2.65 pu | 30
50 | 6.68 pu | 15

###  V-D $\Delta$V Minimization Objective: (D3)

Now we show the result for the third objective function, i.e., $\Delta$V
minimization problem. Here, we solve the decomposed problem (D3) for 10,000
node test system with different DER penetration levels, but with nominal
values of DER generation and loads. The reactive power generation of the DERs
are controlled to make the nodal voltages closer to a reference value of
$V_{ref}=1.00pu$, which is the substation node voltage. For this optimization
problem, simulated 2 different cases – (i) 100% DER penetrations and (ii) 50%
DER penetrations. Again, we set $\alpha=0$ for FPI update.

The optimal result is shown in Fig. 5, where Fig. 5(a) shows that the nodal
voltages after optimization. The higher penetrations of DERs results in closer
node voltages to the substation voltage. Also, for both of the cases, it only
takes 11 macro-iterations to reach convergence (Fig. 5(b)). The objective
value of this cost function is shown in Table I with the total solution time.
The time taken at each iteration is almost same for all the iterations, and
thus not shown in this paper.

(a) Nodal voltages

(b) Convergence at the boundary

Figure 5: Numerical Results for $\Delta$V Minimization Objective for Synthetic
10,000 Node System

(a) Loss minimization objective

(b) DER maximization objective

(c) $\Delta$V minimization objective

Figure 6: Solution Time for Central OPF Problems for Various Objective
Functions

###  V-E Failure of Central Solution

We have solved the large scale OPFs by decomposing the central problem into
several smaller sub-problems, and exchanging the complicating variables. In
this section we solve the same central OPF problems for increasing node
numbers, and show how poorly the central problem scales with the increasing
node numbers. Also, we show when the central OPF fails to solve the OPFs for
all the previously simulated cases. From the Fig. 6, we can see with
increasing DER penetration, it takes higher time to solve the NLP OPFs, and
fails earlier. For example, in case of loss minimization, the NLP solver can
solve the OPF problems for 800 nodes for 50% DER penetration case (Fig. 6(a)).
Similarly, for DER maximization objective, with 20% DER penetration, the
central problem can be solved for no more than 1100 nodes. It is clear that
for any OPFs, the NLP problem cannot be solved for a distribution network with
more than 2000 nodes. Please note, all of the cases have been solved using
knitro with active-set algorithm.

Figure 7: Result for Nodal Decomposition for IEEE-8500 Node Test System

###  V-F Nodal Decomposition

To show the efficacy of the algorithm, we further decomposed the problem at
individual node level, and solved the OPF problem with loss minimization
objective. For this case, we simulated the balanced IEEE-8500 bus test system
with 2522 nodes. While the commercial NLP solvers, such as knitro fails to
optimize this system, the nodal decomposition of the problem speeds up the
whole process significantly. It only takes $\sim 5$ seconds to solve the NLP
OPF with nodal decomposition of the network. This method decomposes the OPF
into 2522 sub-problems, but significantly reduces the variables in each sub-
problem – only 5 variables to solve. It takes 275 macro-iterations to
converge, however, since each sub-problem takes microseconds to solve the
small NLP problems, the overall solution time is only $\sim 5$ seconds. The
convergence is shown in Fig. 7.

##  VI Conclusion

In this paper we have solved the large scale OPF problems for power
distribution systems by decomposing the networks and distributing the problem
into smaller sub-problems. The proposed approach is a generalized
decomposition for radial power distribution system that scales very well for
all general classes of distributed OPF problems. The proposed distributed
approach converges within fewer iterations and in a short-period of time for
large feeders even for the cases where the centralized OPF either takes
significant amount of time or simply fails to converge. We have demonstrated
the successful application of the proposed approach for a synthetic 10,000
node distribution test system with a total $\sim 50,000$ variables, on a
regular CPU, and within reasonable time. To the best of our knowledge, this is
the first work to demonstrate the application of distributed algorithms to
solve OPF problem for the selected very large distribution feeder without
requiring HPC machines. Furthermore, the decomposition is amenable for
implementation on many-core machines; the fast convergence and fewer
communication requirements for the proposed algorithm will lead to significant
advancement in solving large-scale OPF problems for active power distribution
systems. Authors are working on extending the approach to the three-phase
unbalanced D-OPF problems.

## References

  * [1] J. A. Momoh, R. Adapa, and M. El-Hawary, “A review of selected optimal power flow literature to 1993. i. nonlinear and quadratic programming approaches,” IEEE transactions on power systems, vol. 14, no. 1, pp. 96–104, 1999.
  * [2] A. Castillo and R. P. O’Neill, “Survey of approaches to solving the ACOPF,” US Federal Energy Regulatory Commission, Tech. Rep, 2013.
  * [3] R. R. Jha, A. Dubey, C.-C. Liu, and K. P. Schneider, “Bi-level volt-var optimization to coordinate smart inverters with voltage control devices,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 1801–1813, 2019\.
  * [4] L. Gan and S. H. Low, “Convex relaxations and linear approximation for optimal power flow in multiphase radial networks,” in 2014 Power Systems Computation Conference, pp. 1–9, IEEE, 2014.
  * [5] R. R. Jha and A. Dubey, “Network-level optimization for unbalanced power distribution system: Approximation and relaxation,” IEEE Transactions on Power Systems, vol. 36, no. 5, pp. 4126–4139, 2021.
  * [6] R. R. Jha and A. Dubey, “Exact distribution optimal power flow (d-opf) model using convex iteration technique,” in 2019 IEEE Power Energy Society General Meeting (PESGM), pp. 1–5, 2019.
  * [7] D. K. Molzahn, F. Dörfler, H. Sandberg, S. H. Low, S. Chakrabarti, R. Baldick, and J. Lavaei, “A survey of distributed optimization and control algorithms for electric power systems,” IEEE Transactions on Smart Grid, vol. 8, no. 6, pp. 2941–2962, 2017.
  * [8] Q. Peng and S. H. Low, “Distributed optimal power flow algorithm for radial networks, i: Balanced single phase case,” IEEE Transactions on Smart Grid, vol. 9, no. 1, pp. 111–121, 2016.
  * [9] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine learning, vol. 3, no. 1, pp. 1–122, 2011\.
  * [10] W. Zheng, W. Wu, B. Zhang, H. Sun, and Y. Liu, “A fully distributed reactive power optimization and control method for active distribution networks,” IEEE Transactions on Smart Grid, vol. 7, no. 2, pp. 1021–1033, 2015.
  * [11] T. Erseghe, “Distributed optimal power flow using admm,” IEEE transactions on power systems, vol. 29, no. 5, pp. 2370–2380, 2014.
  * [12] E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal power flow for smart microgrids,” IEEE Transactions on Smart Grid, vol. 4, no. 3, pp. 1464–1475, 2013.
  * [13] W. Lu, M. Liu, S. Lin, and L. Li, “Fully decentralized optimal power flow of multi-area interconnected power systems based on distributed interior point method,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 901–910, 2017.
  * [14] G. Cavraro and R. Carli, “Local and distributed voltage control algorithms in distribution networks,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 1420–1430, 2017.
  * [15] A. Bernstein and E. Dall’Anese, “Real-time feedback-based optimization of distribution grids: A unified approach,” IEEE Transactions on Control of Network Systems, vol. 6, no. 3, pp. 1197–1209, 2019.
  * [16] N. Bastianello, A. Ajalloeian, and E. Dall’Anese, “Distributed and inexact proximal gradient method for online convex optimization,” arXiv preprint arXiv:2001.00870, 2020.
  * [17] G. Qu and N. Li, “Optimal distributed feedback voltage control under limited reactive power,” IEEE Transactions on Power Systems, vol. 35, no. 1, pp. 315–331, 2019.
  * [18] X. Hu, Z.-W. Liu, G. Wen, X. Yu, and C. Li, “Branch-wise parallel successive algorithm for online voltage regulation in distribution networks,” IEEE Transactions on Smart Grid, vol. 10, no. 6, pp. 6678–6689, 2019.
  * [19] S. Magnússon, G. Qu, and N. Li, “Distributed optimal voltage control with asynchronous and delayed communication,” IEEE Transactions on Smart Grid, 2020.
  * [20] R. Sadnan and A. Dubey, “Real-time distributed control of smart inverters for network-level optimization,” in 2020 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm), pp. 1–6, IEEE, 2020.
  * [21] R. Sadnan and A. Dubey, “Distributed optimization using reduced network equivalents for radial power distribution systems,” IEEE Transactions on Power Systems, 2021.
  * [22] M. E. Baran and F. F. Wu, “Optimal capacitor placement on radial distribution systems,” IEEE Transactions on power Delivery, vol. 4, no. 1, pp. 725–734, 1989.
  * [23] J. Nocedal, “Knitro: An integrated package for nonlinear optimization,” in Large-Scale Nonlinear Optimization, pp. 35–60, Springer, 2006.


# Swashplateless-elevon Actuation for a Dual-rotor Tail-sitter VTOL UAV

Nan Chen∗, Fanze Kong∗, Haotian Li, Jiayuan Liu, Ziwei Ye,
Wei Xu, Fangcheng Zhu, Ximin Lyu, and Fu Zhang ∗These two authors contribute
equally to this work.N. Chen, F. Kong, H. Li, J. Liu, Z. Ye, F. Zhu, W. Xu,
and F. Zhang are with Department of Mechanical Engineering, The University of
Hong Kong<EMAIL_ADDRESS>{xuweii,fuzhang} @hku.hk. X. Lyu is with School of intelligent systems
engineering, Sun Yat-sen University, Shenzhen<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a novel swashplateless-elevon actuation (SEA) for
dual-rotor tail-sitter vertical takeoff and landing (VTOL) unmanned aerial
vehicles (UAVs). In contrast to the conventional elevon actuation (CEA) which
controls both pitch and yaw using elevons, the SEA adopts swashplateless
mechanisms to generate an extra moment through motor speed modulation to
control pitch and uses elevons solely for controlling yaw, without requiring
additional actuators. This decoupled control strategy mitigates the saturation
of elevons’ deflection needed for large pitch and yaw control actions, thus
improving the UAV’s control performance on trajectory tracking and disturbance
rejection performance in the presence of large external disturbances.
Furthermore, the SEA overcomes the actuation degradation issues experienced by
the CEA when the UAV is in close proximity to the ground, leading to a
smoother and more stable take-off process. We validate and compare the
performances of the SEA and the CEA in various real-world flight conditions,
including take-off, trajectory tracking, and hover flight and position steps
under external disturbance. Experimental results demonstrate that the SEA has
better performances than the CEA. Moreover, we verify the SEA’s feasibility in
the attitude transition process and fixed-wing-mode flight of the VTOL UAV.
The results indicate that the SEA can accurately control pitch in the presence
of high-speed incoming airflow and maintain a stable attitude during fixed-
wing mode flight. Video of all experiments can be found in
youtube.com/watch?v=Sx9Rk4Zf7sQ

##  I INTRODUCTION

Vertical takeoff and landing (VTOL) unmanned aerial vehicles (UAVs) can take
off and land vertically like multi-rotor UAVs and achieve efficient long-range
and high-speed flight similar to fixed-wing UAVs [1]. VTOL UAVs can be
implemented using different configurations, such as tilt-rotor [2, 3], tilt-
wing [4, 5], and tail-sitter [6, 7, 8]. Among these, the tail-sitter UAV is
advantageous due to its ability to transition its attitude to enter fixed-wing
flight mode without additional tilting mechanisms, resulting in a simpler and
more compact structure [9].

Tail-sitter UAVs can be controlled using various actuator configurations, such
as quad-rotor type [10, 11] and dual-rotor type with additional control
surfaces (e.g., elevons) [12, 13, 14]. Compared with the quad-rotor type, the
dual-rotor type uses fewer motors and propellers, making it less expensive,
lighter, and more portable. However, the limited deflection range and small
size of the additional control surfaces may lead to saturation, which
constrain the achievable performance or even cause instability [15]. In
particular, when the UAV hovers in multi-rotor mode, wind disturbances can
significantly affect the UAV’s attitude due to its large wing area. These
disturbances often act on the pitch and yaw directions of the dual-rotor-type
VTOL UAVs, which are controlled by the elevons (call it conventional elevon
actuation (CEA)). One problem of the CEA is that large control efforts in both
directions can induce elevons’ deflection sharply, causing actuation
saturation that further degrades the control performance or even destabilizes
the system. One approach to mitigate this problem is to place elevons at the
high-speed airflow region of the propellers by reducing the distance between
rotors and elevons [16]. However, this design also decreases the airflow speed
flowing through the main wing, resulting in lift loss. Overall, wind-
resistance is one of the major challenges in moving the dual-rotor tail-sitter
VTOL UAV towards practical use.

Figure 1: Dual-rotor tail-sitter VTOL UAV, named Hong He. (a) The UAV hovers
at the multi-rotor mode. Main components on the front side of the UAV are
labeled. (b) The UAV flies at the fixed-wing mode. Main components on the back
side of the UAV and the definition of coordinate frame are shown.

Another problem of CEA is encountered in the take-off process. Due to the
constraints of the UAV footprint and landing stability, the landing gear
cannot be very high, resulting in a short distance between the elevons and the
ground. Before take-off, the downwash airflow generated from the propellers
can be easily affected by the ground, leading to a decreased or even
disappeared control moment of the elevons [17]. As a result, the pitch
attitude may become unstable, potentially causing the UAV to crash.

Motivated by these problems, in this study, we propose a swashplateless-elevon
actuation (SEA) which employs a swashplateless mechanism with improved
structure to dual-rotor tail-sitter VTOL UAVs. The swashplateless mechanism is
a passive structure mounted on the motors. It generates lateral moment by
controlling the high-frequency acceleration and deceleration of the motor
without requiring additional actuators [18, 19]. For VTOL UAVs, it can produce
an extra control moment in the pitch direction, which can dramatically
alleviate elevons’ deflection saturation caused by large disturbances applied
in both pitch and yaw directions or by large moments required by controllers,
leading to a better disturbance rejection performance and higher
maneuverability. Furthermore, the moment generated by the swashplateless
mechanism is independent of the propeller’s airflow, which enables a more
stable take-off process for this type of UAV. We compare the SEA with the CEA
in various experiments, which are described in the following sections.

##  II SYSTEM DESIGN

###  II-A UAV Structure and Avionics

We designed and manufactured a dual-rotor tail-sitter VTOL UAV named Hong He
as shown in Fig. 1. For realizing large-scale scanning and mapping, it is
equipped with a Livox AVIA LiDAR and an onboard computer Manifold 2-C. The
LiDAR is capable of achieving long-distance scanning up to 400 m, making it
suitable for outdoor large-scale scenes. To meet the requirements of long-
range flight and control performance, we choose a canard wing layout and
optimize the aircraft’s aerodynamic shape design. For this prototype, our
mechanical design concept emphasizes ease of assembly and maintenance, and
hence we use formed wings and 3D printed materials (e.g., PLA, Nylon) to
manufacture the UAV, making it simple to assemble.

The avionics control system is centered around a Pixhawk mini 4 flight
controller attached to the fuselage. It is connected to two Electronic Speed
Controllers (ESCs) (model T-MOTOR F35A 3-6S) for driving two motors (model
T-MOTOR MN5006 KV450), two servos (model KST DS215MG V3.0) to control the
deflection of two elevons, and two magnetic encoders (model AS5600) to measure
the rotor angles of two motors. The final designed aircraft has a wingspan of
107 cm, weighs about 2.25 kg, and has a cruising speed of approximately 9 m/s.
Detailed parameters can be found in Table I.

As shown in Fig. 1 (b), the origin of UAV’s coordinate frame (i.e. body frame)
is attached to the UAV’s center of mass. Beside the $x$ axis is perpendicular
to the plane of main wing, both $y$ and $z$ axes are located in the plane of
main wing, but the $y$ axis is parallel with the main wing while the $z$ axis
points to the tail of UAV. The roll, pith, and yaw are defined in UAV’s multi-
rotor mode, and hence their rotation are along with the $x$, $y$, and $z$
axes, respectively. All the experimental data in this paper follow this
definition.

###  II-B Swashplateless Mechanism

The swashplateless mechanism, originally proposed in [20], can realize
functions similar to the cyclic blade pitch control mechanism of swashplates
that has been widely used in helicopters, providing both thrust and moment.
The moment of the swashplateless mechanism comes from the unbalanced thrust of
blades, induced by cyclic blade pitch changes. Unlike traditional swashplates
driven by additional servos, the swashplateless mechanism is entirely passive.
As shown in Fig. 2, two passive hinges connect side hubs to the central hub
asymmetrically. Through periodic acceleration and deceleration of the motor
(i.e., motor speed modulation), the unsymmetrical hinges rotate due to blade
inertia, leading to different pitch angle changes of positive and negative
blades. The blade with an increased pitch angle produces more thrust, while
the blade with a decreased pitch angle produces less thrust, resulting in a
net moment being generated. Detailed working principles are detailedly
introduced in [20, 21, 22].

The propulsion system shown in Fig. 2 mainly composes of a brushless DC motor,
an improved swashplateless mechanism mounted on top of the motor’s rotor, a
pair of propeller blades, and a magnetic encoder mounted on the bottom of the
motor’s shaft. Compared with the original design, the improved structure of
swashplateless mechanism includes additional ball bearings and pressure
bearings to reduce the friction caused by the high-frequency rotation of the
hinges, providing a smoother and more linear moment output and a faster
response.

TABLE I: Specifications of the VTOL UAV Hong He. Item | Parameter
---|---
Wing span | 107 cm
Total weight (with battery) | 2.25 kg
Thrust weight ratio | 2.5
Cruise speed | 9 m/s
Cruise current | 9 A
Hover current | 12 A
Battery | 6S, 2600 mAh
Figure 2: Components of one set of Hong He’s propulsion system, mainly
including a motor, a swashplateless mechanism, and two blades. Figure 3:
Control system overview. All modules in the control part are inherited from
the standard PX4 and the mixer module of actuation part is additionally
implemented into the PX4.

##  III Control and Actuation

###  III-A Overview of Control System

The overview of Hong He’s control system with the proposed SEA is shown in
Fig. 3. The external inputs of the position & velocity controller can be
switched to two sources: remote controller (i.e., manual mode) and onboard
computer (i.e., autonomous mode). In the manual mode, the remote controller
directly generates the desired velocity. In the autonomous mode, the onboard
computer produces desired position as the input of the position controller and
the desired velocity & acceleration as the feed-forward terms.

The control section (i.e., the orange zone) is the standard realization of PX4
while the actuation section (i.e., the green zone) are specifically designed
for actuating the swashplateless mechanism and hence should be focused on. The
operating frequency of the SEA mixer is dependent on the output frequency of
the angular velocity controller (i.e., the measurement frequency of the
inertial measurement unit (IMU)). Therefore, we set the measurement frequency
of the IMU to 1 kHz, allowing the mixer run at the same frequency for
processing the 910-Hz measurement from the two magnetic encoders without
dropping any measurement. The outputs of the mixer are motor commands and
servo commands, which are sent to the ESCs using the DShot600 communication
protocol and to the servo by a standard 50-Hz PWM signal. The DShot600
protocol has a very short communication delay of 26.7 us, making it suitable
for high-frequency speed modulation of the motor.

###  III-B Actuation Principles of SEA and CEA

The SEA and CEA exhibit identical actuation principles with respect to thrust
generation and the control of roll and yaw moments. The two propellers, which
are driven by motors, contribute to the total UAV thrust. Differential thrust
produced by the propellers induces the roll control moment, while the airflow
deflection by the elevons create the yaw control moment. However, there exists
a difference between the SEA and CEA in terms of pitch control moment
generation. The elevons contribute to the pitch control moment in the CEA,
whereas the SEA employs swashplateless mechanisms for this purpose. This
actuation principle facilitates decoupling of the attitude control of pitch
and yaw.

In the CEA, the mixer blends the pitch and yaw moments to regulate the
elevons’ deflection through servo commands. Conversely, in the realization of
SEA (i.e., the green zone in Fig. 3), the pitch and yaw controls are directed
towards separate actuators (i.e., the motors and the elevons). The SEA can
mitigate elevon saturation by decoupling the generation of pitch control
moment, thereby enhancing control and disturbance rejection performance. This
feature is of particular significance, given that disturbances are frequently
applied in the pitch and yaw directions due to the UAV’s large main wing
surface area.

###  III-C Design of SEA Mixer

####  III-C1 Cyclic speed control of motor

Figure 4: Coordinate frame definition for the swashplateless mechanism of the
two motors shown in the top view of the UAV. The UAV coordinate frame is also
labeled for reference.

Before introducing the mixer, the cyclic speed control of motor is needed to
be described since it defines a part of the output variables of the mixer. To
prevent undesired vibration caused by motor speed modulation, a sinusoidal
signal is employed for two motors. The total motor throttle $U_{i},i=1,2$ is
designed as

$\displaystyle\left[\begin{matrix}U_{1}\\\
U_{2}\end{matrix}\right]=\left[\begin{matrix}C_{1}\\\
C_{2}\end{matrix}\right]+\left[\begin{matrix}A_{1}&0\\\
0&A_{2}\end{matrix}\right]\left[\begin{matrix}\vspace{0.1cm}\cos(\theta_{1}-\phi_{1}+\gamma_{0})\\\
\cos(\theta_{2}-\phi_{2}-\gamma_{0})\end{matrix}\right],$ (1)

where $C_{i}$, $A_{i}$, $\theta_{i}$, and $\phi_{i}$ are the nominal throttle,
sinusoidal amplitude, motor’s rotor angle measured by the magnetic encoder,
and the moment direction in the UAV coordinate frame of motor $i$, $i=1,2$,
respectively, the $\gamma_{0}$ is a positive constant for compensating the
angle delay caused by blades’ inertia and can be calibrated in advance through
experiment data from a test stand of the swashplateless mechanism.

Despite the fact that the motors rotate in opposite directions, we can define
the same coordinate frame for both motors, as shown in Fig. 4. The absolute
angle of the rotor, $\theta_{i}$, is defined as the angle between the $x$-axis
and the positive blade, regardless of the rotation direction. In the motor
coordinate frame, the moment direction is represented as $\phi_{i}+\pi/2$.
However, due to the 90∘ rotation between the UAV frame and the motor frame,
$\phi_{i}$ can represent the actual moment direction in the UAV frame. The
sign of $\gamma_{0}$ is determined by the rotation direction since the delay
direction is opposite to the rotation direction. As only the pitch angle is
controlled by the moment of the swashplateless mechanism, $\phi_{i}$ can be
simply set as $0$ or $\pi$, depending on the direction of the desired pitch
moment.

####  III-C2 Mixer mapping

The mixer is designed to calculate all actuator inputs (e.g., motor throttles,
servo angles) according to the desired thrust $f_{T,d}^{\mathcal{B}}$ and
moment $\bm{\tau}^{\mathcal{B}}_{d}$. Given the geometry of the mechanical
structure of the UAV, the mapping from actuator outputs to the body thrust and
moment are

$\displaystyle\left[\begin{matrix}\vspace{0.1cm}f_{T}^{\mathcal{B}}\\\
\vspace{0.1cm}\tau_{x}^{\mathcal{B}}\\\
\vspace{0.1cm}\tau_{y}^{\mathcal{B}}\\\
\tau_{z}^{\mathcal{B}}\end{matrix}\right]=\left[\begin{matrix}1&1&0&0&0&0\\\
-L&L&0&0&0&0\\\ 0&0&1&1&0&0\\\ 0&0&0&0&1&1\\\
\end{matrix}\right]\left[\begin{matrix}T_{1}\\\ T_{2}\\\ \tau_{s,1}\\\
\tau_{s,2}\\\ \tau_{e,1}\\\ \tau_{e,2}\end{matrix}\right],$ (2)

where $T_{i}$ and $\tau_{s,i}$ are the thrusts and moments generated by the
motor $i,i=1,2$, respectively. $\tau_{e,j}$ are the moments generated by the
elevon $j,j=1,2$, respectively. $L$ is the body-$y$-axis distance between
UAV’s center of mass and the motor.

Assuming that the actuator outputs are approximately linear to the actuator
inputs, leads to

$\displaystyle\left[\begin{matrix}T_{1},T_{2},\tau_{s,1},\tau_{s,2},\tau_{e,1},\tau_{e,2}\end{matrix}\right]^{T}=$
(3)
$\displaystyle\text{diag}(K_{t},K_{t},K_{a},K_{a},K_{e},K_{e})\left[\begin{matrix}C_{1},C_{2},A_{1},A_{2},\delta_{1},\delta_{2}\end{matrix}\right]^{T},$

where $\delta_{j}$ are the angle command of servo $j$, $j=1,2$. $K_{t}$,
$K_{a}$, and $K_{e}$ are proportion coefficients, representing the conversions
from throttle to thrust, from sinusoidal amplitude to swashplateless mechanism
moment, and from servo angle to elevon moment, respectively. Here, the
$C_{1}$, $C_{2}$, $A_{1}$, and $A_{2}$ are required by the cyclic speed
control of motor.

Given the desired thrust $f_{T,d}^{\mathcal{B}}$ and the moment vector
$\bm{\tau}_{d}^{\mathcal{B}}=\left[\begin{matrix}\tau_{x,d}^{\mathcal{B}}&\tau_{y,d}^{\mathcal{B}}&\tau_{z,d}^{\mathcal{B}}\end{matrix}\right]^{T}$,
combining (2) and (3) and assuming the equal distribution of moments, we can
finally determine the desired actuator inputs as

$\displaystyle\left[\begin{matrix}C_{1,d}\\\ C_{2,d}\\\ A_{1,d}\\\ A_{2,d}\\\
\delta_{1,d}\\\
\delta_{2,d}\end{matrix}\right]=\left[\begin{matrix}\vspace{0.1cm}\frac{1}{2K_{t}}&\frac{-1}{2LK_{t}}&0&0\\\
\vspace{0.1cm}\frac{1}{2K_{t}}&\frac{1}{2LK_{t}}&0&0\\\
\vspace{0.1cm}0&0&\frac{1}{2K_{a}}&0\\\
\vspace{0.1cm}0&0&\frac{1}{2K_{a}}&0\\\
\vspace{0.1cm}0&0&0&\frac{1}{2K_{e}}\\\
\vspace{0.1cm}0&0&0&\frac{1}{2K_{e}}\end{matrix}\right]\left[\begin{matrix}\vspace{0.1cm}f_{T,d}^{\mathcal{B}}\\\
\vspace{0.1cm}\tau_{x,d}^{\mathcal{B}}\\\
\vspace{0.1cm}\tau_{y,d}^{\mathcal{B}}\\\
\tau_{z,d}^{\mathcal{B}}\end{matrix}\right],$ (4)

where the $\delta_{j,d}$ are directly sent to the servo $j$ by PWM signal. The
$C_{i,d}$ and $A_{i,d}$ are used to generate the total motor throttle
$U_{i,d}$ sent to the ESCs. In (4), the desired moments of pitch and yaw are
decoupled into different actuators (i.e., the motors’ swashplateless
mechanisms and elevons).

##  IV Experiments

This section presents experimental results comparing the performance of the
SEA and the CEA in several aspects, including take-off, trajectory tracking,
and disturbance rejection. Additionally, tests of attitude transition and
fixed-wing-mode flight are also performed using the SEA. To ensure a fair
comparison, the controllers’ parameters of angular rate, attitude, and
position are tuned carefully in the best effort to optimize the performance
for both SEA and CEA. The tuning process are conducted by analyzing the logged
actual flight data to ensure fast responses while keeping small overshoots.
All experiments can be found in an accompanying video uploaded in
youtube.com/watch?v=Sx9Rk4Zf7sQ

###  IV-A Comparison of Take-off Performance

First, three take-off experiments are conducted on the UAV with attitude
control (i.e., temporarily turn off the position controller). They are CEA
with ground take-off, CEA with take-off on a lifted pedestal, and SEA with
ground take-off. The results are shown in Fig. 5 (a)-(c), respectively. In
Fig. 5 (a), the pitch angle unexpectedly increases to a maximum value of 8.6∘,
although the desired attitude is still zero during take-off. The large pitch
error during ground take-off is caused by the small distance between elevons
and the ground. Constrained by both the footprint and landing stability, the
landing gear cannot be higher, making the elevons easily affected by the near-
ground airflow and deteriorating their control effect. In contrast, Fig. 5 (b)
shows the same take-off process on a lifted pedestal, and the pitch error
significantly decreases to a maximum value of 2.1∘. The lifted pedestal
alleviates the effect caused by the near-ground airflow, enabling the elevons
to maintain a control effect during the whole take-off process. In the case of
SEA shown in Fig. 5 (c), because the pitch is fully controlled by the
swashplateless mechanism, which is not affected by the near-ground airflow,
the pitch error during take-off is minimal (1.4∘).

By turning the position controller off, one can manually control the throttle
and allow the UAV to take-off quickly. In such a case, the pitch error can be
less than 10∘ and does not cause significant disadvantages. However, when the
UAV takes off with position control, the throttle is decided by the
calculation of position controller. If a smooth take-off is expected, the
take-off throttle may increase gradually, resulting in a significant pitch
angle and position error in the CEA. Figure 6 shows the UAV’s position and
attitude during take-off in position control mode. A significant position
error (0.79 m) is caused by the large pitch error (14.7∘) when using the CEA.
For the SEA, the position error is only 0.04 m, which is merely 5.1% of the
CEA. Therefore, the take-off performance of the SEA control is much better
than that of the CEA, leading to a smooth take-off even when the UAV is close
to the ground, no matter the position controller is turned on or off.

Figure 5: Comparison of pitch angle error during take-off in attitude control
mode with different conditions. Figure 6: Comparison of $x$-position error and
pitch angle error during take-off in position control mode with different
actuation approach.

###  IV-B Comparison of Trajectory Tracking Performance

The tracking performance of both SEA and CEA are compared by conducting a
figure-of-eight trajectory tracking experiment, as shown in Fig. 7. The length
and width of the figure-of-eight trajectory are 2 m and 1 m, respectively, and
the trajectory is executed for four cycles. The trajectory period is 5 s, but
in the first and last cycles, the trajectory period is extended to 7.5 s to
achieve smooth acceleration and deceleration. During tracking, the yaw angle
is commanded to maintain the body $x$-axis same as the horizontal component of
the velocity direction, which can generate more attitude control effort to
evaluate the differences between SEA and CEA. The results in Fig. 8 show that
the tracked trajectory of SEA is more consistent than that of CEA. The
absolute errors of position norm and attitude yaw during tracking are shown in
Fig. 9. For the absolute error of position norm, the median value of SEA and
CEA are 14.17 cm and 15.35 cm, respectively, and the maximum values for both
are 38.91 cm and 40.64 cm, respectively, indicating that the SEA has a
slightly better performance than the CEA in the position control. For the
absolute error of attitude yaw, the median value of SEA and CEA were 19.5∘ and
27.1∘, respectively, and the maximum value for both were 87.9∘ and 119.3∘,
respectively, indicating that the SEA has relatively obvious advantage in yaw
control than the CEA during the fast trajectory tracking.

Figure 7: The overlaid snapshots of Hong He when it tracks a 3D figure-of-
eight trajectory with continuous yaw rotation. Figure 8: The reference
trajectory and tracked trajectory in the trajectory tracking experiment.
Figure 9: The errors of position norm and yaw angle during tracking the 3D
figure-of-eight trajectory. The central mark indicates the median, and the
bottom and top edges of the box represent the 25th and 75th percentiles,
respectively. The whiskers extend to the maximum and minimum value.

###  IV-C Comparison of Disturbance Rejection

Three experiments are performed to evaluate the disturbance rejection
performance of the SEA and the CEA: (i) hovering under balanced wind
disturbance, (ii) forward position steps under unbalanced wind disturbance,
and (iii) lateral position steps under unbalanced wind disturbance. The
position information is provided by a motion capture system (Vicon).

####  IV-C1 Hover under balanced wind disturbance

As shown in Fig. 10 (a), two fans are placed in parallel and located 1 m from
the UAV’s hover position. The fans are off initially and are then turned on to
produce a wind gust of about 4.5 m/s in this distance. After the UAV reaches a
stable pose, the fans are turned off to cancel the wind disturbance. The
results are shown in Fig. 10 (b). Since the wind disturbance is applied on the
body $x$-axis, we can only focus on the errors on $x$-axis position and pitch
angle. For the $x$-axis position, both the median value and maximum value of
SEA are smaller than that of CEA, which are 1.77 cm versus 3.83 cm and 6.60 cm
versus 7.92 cm, respectively. The results show the SEA can maintain the UAV’s
position better. For the pitch angle errors, although the CEA has slightly
smaller median values than the SEA (i.e., 1.09∘ versus 1.32∘), the maximum
value of CEA are relatively larger (i.e., 5.78∘ versus 4.08∘). This shows that
both actuation approaches have good pitch angle control performance under wind
disturbance but the SEA can realize a smaller error band.

Figure 10: Hover under balanced wind disturbance. (a) Experimental setup for
the disturbance rejection test of fan wind. (b) errors of $x$-axis position
and pitch angle of SEA and CEA during a balanced wind disturbance.

####  IV-C2 Forward position steps under unbalanced wind disturbance

Maintaining good response in the presence of external disturbances is more
difficult than without disturbance, as the actuators need to actively suppress
the disturbance, leading to reduced actuator margins for tracking controller
commands. To verify actual performance under disturbances, a position step
response in the $x$ direction is performed with a fan wind disturbance being
applied to a single side of the main wing.

In the experiment, the UAV first hovers 1 m from a fan that has been turned
on. Then, the UAV steps 1 m in the same direction as the wind airflow and
finally steps back to the original hover position. Since the wind area of the
UAV is large, we pay more attention to the attitude response during the wing
disturbance. The results of attitude response is shown in Fig. 11. When the
UAV steps away from the fan, the wind disturbance decreases. In this
condition, and attitude response of the two actuation approaches are similar.
However, when the UAV steps back to the original hover position, the wind
disturbance increases. Because the elevons of CEA need to deflect an angle to
generate control moment both in pitch and yaw directions, they cannot provide
sufficient moment to resist the wind disturbance on these two directions.
Consequently, when using CEA, obvious oscillation occurred in the converging
process of pitch angle and significant shaking in yaw angle (maximum error of
31.2∘). Since large attitude error occurred in pitch and yaw, the wind
disturbance affects roll more easily due to the changing of the windward side,
leading to a maximum error of 8.2∘. In contrast, no oscillation or obvious
error occurs in the attitude response of the SEA during the entire process,
indicating that it can achieve more actuator margin and maintain good attitude
control performance even with external disturbance.

Figure 11: Attitude responses (left panel) and their corresponding errors
(right panel) when applied a step command in body $x$ direction under an
unbalanced wind disturbance on the main wing. Figure 12: Attitude responses
(left panel) and their corresponding errors (right panel) when applied a step
command in body $y$ direction under an unbalanced wind disturbance on the main
wing.

####  IV-C3 Lateral position steps under unbalanced wind disturbance

In addition to the forward and backward position step experiment in the body
$x$ direction, we also conducted a lateral position step experiment in the
body $y$ direction. The experiment setup is same as before. The UAV initially
hovered 1 m away from a fan that had been turned on, then stepped 1 m away
from the fan in the direction which is orthogonal to the direction of the wind
airflow, and finally stepped back to the original hover position. The attitude
response results of this experiment are shown in Fig. 12.

When the UAV stepped away from the fan, the wind disturbance decreased. In
this condition, the attitude control of the two actuation approaches are still
similar. However, when the UAV stepped back to the original hover position,
the wind disturbance increased. Similar to the forward step experiment, the
CEA showed obviously larger yaw angle shaking (maximum error of 22.1∘) than
the SEA (maximum error of 8.4∘), further verifying that the SEA has better
disturbance rejection performance than the CEA.

###  IV-D Transition and Fixed-wing Mode Flight

We also conducted experiments to verify the performance of the proposed SEA in
the fixed-wing mode flight as shown in Fig. 1 (b). The results are shown in
Fig. 13. The definition of UAV coordinate frame in the fixed-wing mode is the
same as that of the multi-rotor mode. The pitch angle of fixed-wing mode is
set at -65∘ (i.e., the angle between main wing and horizontal plane is 35∘).
During the transition, the attitude controller tracked the desired attitude
accurately without any overshoot, showing that the swashplateless mechanism
has a fast response and can generate sufficient moment. In the fixed-wing mode
flight, the UAV accelerated continuously and finally reached a speed of 9.6
m/s. The pitch angle had no vibration and always remained at 65∘, indicating
that the swashplateless mechanism can work well in the fixed-wing mode.
However, the roll and yaw, which are controlled by the motors and elevons,
respectively, had slight vibrations around the desired values. The maximum
errors in roll and yaw appeared when the UAV just completed the transition
process, which may have been caused by the airflow disturbance. Nonetheless,
these errors were not obvious, only 5.3∘ in roll and 3.8∘ in yaw. It is should
be noted that the CEA-based fixed-wing mode flight of a VTOL UAV has been
verified previously. Hence, we did not present a quantitative comparison
experiment here and only verify the feasibility of the fixed-wing mode flight
based on SEA.

Figure 13: Desired attitude, actual attitude, and flight speed during
transition and fixed-wing mode flight using SEA.

##  V Conclusion

In this paper, an actuation approach called SEA is proposed for dual-rotor
VTOL UAV to decouple pitch and yaw control and improve take-off and
disturbance rejection performance. The proposed SEA-based UAV showed reduced
pitch and position errors during take-off compared to the CEA-based UAV, which
had noticeable errors due to ground-distorted airflow. The control performance
of both SEA and CEA are evaluated by tracking a 3D figure-of-eight trajectory
with continuous yaw angle rotation, showing that the SEA has less error both
in position and yaw angle. Disturbance rejection performance was evaluated by
the experiment of hovering in the balanced wind gust produced by two fans. The
SEA exhibited better performance in position and pitch angle, indicating that
SEA is more robust in an environment with wind gust. Step response experiments
under unbalanced wind disturbance showed that the SEA outperformed the CEA
with obviously lower attitude errors. These experiments validate that the SEA
can mitigate actuator saturation by decoupling the actuation of pitch and yaw,
and improve the performances of both tracking control and disturbance
rejection. Finally, we validate the SEA in the transition process and fixed-
wing mode flight in an outdoor environment, demonstrating its capability to
maintain a stable attitude for a VTOL UAV in the presence of high-speed
incoming airflow.

## ACKNOWLEDGMENT

This work is supported by the Grants Committee Early Career Scheme of The
University of Hong Kong under Project 27202219, General Research Fund of Hong
Kong under project 17206920, and a DJI research donation.

## References

  * [1] M. S. Gandhi, L. Whitcher, E. Theodorou, and E. N. Johnson, “Practical system identification for small vtol unmanned aerial vehicle,” in _AIAA Scitech 2019 forum_ , 2019, p. 1982.
  * [2] Y. Chen, J. Zheng, Z. Shen, P. Xiao, M. Tong, L. Guo, and Y. Xin, “Design, analysis and tests of a tri–tiltrotor aircraft,” in _AIAA Aviation 2019 Forum_ , 2019, p. 2883.
  * [3] D. N. Cardoso, S. Esteban, and G. V. Raffo, “A new robust adaptive mixing control for trajectory tracking with improved forward flight of a tilt-rotor uav,” _ISA transactions_ , vol. 110, pp. 86–104, 2021.
  * [4] D. Rohr, T. Stastny, S. Verling, and R. Siegwart, “Attitude and cruise control of a vtol tiltwing uav,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 3, pp. 2683–2690, 2019.
  * [5] L. M. Sanchez-Rivera, R. Lozano, and A. Arias-Montano, “Development, modeling and control of a dual tilt-wing uav in vertical flight,” _Drones_ , vol. 4, no. 4, p. 71, 2020.
  * [6] X. Lyu, H. Gu, Y. Wang, Z. Li, S. Shen, and F. Zhang, “Design and implementation of a quadrotor tail-sitter vtol uav,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 3924–3930.
  * [7] W. Xu, H. Gu, Y. Qing, J. Lin, and F. Zhang, “Full attitude control of an efficient quadrotor tail-sitter vtol uav with flexible modes,” in _2019 International Conference on Unmanned Aircraft Systems (ICUAS)_. IEEE, 2019, pp. 542–550.
  * [8] B. Li, W. Zhou, J. Sun, C. Wen, and C. Chen, “Model predictive control for path tracking of a vtol tailsitter uav in an hil simulation environment,” in _2018 AIAA modeling and simulation technologies conference_ , 2018, p. 1919\.
  * [9] W. Xu and F. Zhang, “Learning pugachev’s cobra maneuver for tail-sitter uavs using acceleration model,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3452–3459, 2020.
  * [10] A. Oosedo, S. Abiko, A. Konno, and M. Uchiyama, “Optimal transition from hovering to level-flight of a quadrotor tail-sitter uav,” _Autonomous Robots_ , vol. 41, no. 5, pp. 1143–1159, 2017.
  * [11] B. Li, W. Zhou, J. Sun, C.-Y. Wen, and C.-K. Chen, “Development of model predictive controller for a tail-sitter vtol uav in hover flight,” _Sensors_ , vol. 18, no. 9, p. 2859, 2018.
  * [12] J. Sun, B. Li, C.-Y. Wen, and C.-K. Chen, “Design and implementation of a real-time hardware-in-the-loop testing platform for a dual-rotor tail-sitter unmanned aerial vehicle,” _Mechatronics_ , vol. 56, pp. 1–15, 2018.
  * [13] E. A. Tal and S. Karaman, “Global trajectory-tracking control for a tailsitter flying wing in agile uncoordinated flight,” in _AIAA AVIATION 2021 FORUM_ , 2021, p. 3214.
  * [14] R. H. Stone, P. Anderson, C. Hutchison, A. Tsai, P. Gibbens, and K. Wong, “Flight testing of the t-wing tail-sitter unmanned air vehicle,” _Journal of Aircraft_ , vol. 45, no. 2, pp. 673–685, 2008.
  * [15] J. Zhong, B. Song, Y. Li, and J. Xuan, “L1 adaptive control of a dual-rotor tail-sitter unmanned aerial vehicle with input constraints during hover flight,” _IEEE Access_ , vol. 7, pp. 51 312–51 328, 2019.
  * [16] W. Wang, J. Zhu, M. Kuang, X. Yuan, Y. Tang, Y. Lai, L. Chen, and Y. Yang, “Design and hovering control of a twin rotor tail-sitter uav,” _Science China Information Sciences_ , vol. 62, no. 9, pp. 1–3, 2019.
  * [17] P. A. Baker, W. G. Schweikhard, and W. R. Young, _Flight evaluation of ground effect on several low-aspect-ratio airplanes_. National Aeronautics and Space Administration, 1970.
  * [18] J. Paulos and M. Yim, “Flight performance of a swashplateless micro air vehicle,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 5284–5289.
  * [19] Paulos, James and Yim, Mark, “Cyclic blade pitch control without a swashplate for small helicopters,” _Journal of Guidance, Control, and Dynamics_ , vol. 41, no. 3, pp. 689–700, 2018.
  * [20] J. Paulos and M. Yim, “An underactuated propeller for attitude control in micro air vehicles,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 1374–1379.
  * [21] N. Chen, F. Kong, W. Xu, Y. Cai, H. Li, D. He, Y. Qin, and F. Zhang, “A self-rotating, single-actuated uav with extended sensor field of view for autonomous navigation,” _Science Robotics_ , vol. 8, no. 76, p. eade4538, 2023.
  * [22] Y. Qin, N. Chen, Y. Cai, W. Xu, and F. Zhang, “Gemini ii: Design, modeling, and control of a compact yet efficient servoless bi-copter,” _IEEE/ASME Transactions on Mechatronics_ , 2022.