# Simulation-Integrated Distributed Optimization for Unbalanced Power
Distribution Systems

Rabayet Sadnan , Nathan Gray , Anjan Bose, and Anamika Dubey, Rabayet Sadnan,
Nathan Gray, Anjan Bose and Anamika Dubey are with the Department of EECS,
WSU, Pullman, WA-99164, USA; e-mail<EMAIL_ADDRESS>
###### Abstract

Distributed optimization methods have been extensively applied for the
optimization of electric power distribution systems, especially for grid-edge
coordination. Existing distributed optimization algorithms applied to power
distribution systems require many communication rounds among the distributed
agents and may pose convergence challenges in difficult nonlinear settings.
The communication network parameters also significantly impact the algorithm’s
performance. In this paper, we propose a scalable, equivalent network
approximation-based, distributed optimization algorithm that employs
simulation within optimization using the system’s digital twin (DT) to solve
the optimal power problems (OPF) for a three-phase unbalanced distribution
system. The proposed approach is implemented using a cyber-physical co-
simulation platform to validate the robustness of the proposed distributed
algorithm under stressed communication. The proposed approach is thoroughly
validated using the IEEE 123-bus test system.

###### Index Terms:

Distributed algorithms, optimal power flow, distributed energy resources
(DERs), co-simulation, digital twin.

##  I Introduction

Optimization of power distribution systems has gained significant attention
with growing penetrations of controllable distributed energy resources (DERs).
Traditionally, the distribution system is centrally managed using an advanced
distribution management system (ADMS) that collects the measurement data and
use it to coordinate DERs and other grid-edge technologies. With the
increasing number of such controllable assets and growing system complexity, a
centrally managed system suffers from several challenges including poor
scalability, and susceptibility to a single-point failures [1, 2, 3]. The
additional nonlinearity introduced by the three-phase and unbalanced power
distribution systems may pose convergence issues for the centralized
optimization problems and may result in failure to attain the optimal solution
within a reasonable time [3]. Besides, a centrally managed system is
impractical for emerging distribution system typically composed of several
geographically separated areas controlled by different entities that may not
want to share model and/or data. To address these challenges, distributed
coordination and optimization paradigm has gained significant attention for
grid-edge coordination.

The existing literature has shown the applicability of the distributed
optimization paradigm in reducing the computational requirements for optimal
power flow (OPF) problems and in mitigating other challenges including
ownership boundaries and avoiding single-point failures [3]. To this end,
general distributed optimization techniques, such as, Alternating Direction
Method of Multipliers (ADMM), Auxiliary Problem Principle (APP), etc. have
been adopted to solve distributed OPF problems for power distribution systems
[4, 5, 6]. Distributed OPF algorithms involve two iterative processes: (1)
sub-problem optimization involving micro-iterations, where distributed agents
individually solve optimization problems specific to their areas, (2)
distributed coordination requiring macro-iterations, where distributed agents
exchange the boundary variables typically with their immediate neighbors. The
overall compute complexity of the algorithm depends upon the total number
micro- and macro- iterations taken to achieve network-level convergence. The
nature of the optimization problem dictates micro-iterations needed to solve
each sub-problem; the added nonlinearity with three-phase OPF may pose
convergence and compute challenges for this step. The distributed coordination
approach dictates the number of macro-iterations needed for network-wide
convergence; traditional algorithms take 100s of rounds of communication
exchange or macro-iterations resulting in a slow converge. This calls for
algorithmic advances in distributed OPF algorithms to realistically apply them
to a real-world distribution system.

Although several recent methods attempted to reduce the computational
requirements for local subproblems [7, 8, 9], challenges with convergence and
speed persist due to added nonlinearities in three-phase unbalanced OPF models
for power distribution system. Earlier optimization with simulation has been
introduced for general complex optimization problems that can leverage the
already developed complex models of the system, and avoid including the
detailed system model in the optimization problem formulation [10, 11, 12].
Although this approach has not been directly applied to solve OPFs for power
systems, the use of high-fidelity simulation models for design optimization
has been reported [13, 14]. This high-fidelity simulation model is referred as
a digital twin (DT) of the physical power system and has been employed to aid
with situational awareness and decision-making for power grid operations [15,
16, 17]. The DT models provide an ability to perform high-fidelity simulations
of the grid that can potentially help reduce the complexity of difficult OPF
problems. However, the use of DT models has been limited to evaluating the
physical system in advance for potential actions; their utility in directly
assisting with difficult power system optimization problems has not been
investigated.

Another challenge with traditional algorithms to solve distributed OPF
algorithms emerge from a large number of macro-iterations ($\geq 100$
iterations) needed to converge for relatively small and/or balanced systems
[18, 19, 3]. Thus, these algorithms require an expensive communication
infrastructure with high bandwidth and low communication delays. Although
real-time feedback-based online distributed algorithms address some of the
challenges [20, 21], they still require thousands of time-steps to track the
optimal solution, resulting in a sub-optimal system operation and constraint
violations. Recently, authors proposed an algorithm based on the structure of
the power flow problem for a radial distribution system that drastically
reduces the number of macro-iterations [22]. However, the sub-problems are
modeled as nonlinear programming (NLP) problems that may pose convergence and
compute challenges. Moreover, the performance of distributed OPF algorithms
has not been evaluated in a realistic cyber-physical environment with
appropriate models for communication systems [22, 18, 21].

The objective of this paper is to develop a fast, computationally tractable,
and communication-efficient distributed OPF (D-OPF) algorithm for unbalanced
power distribution system. Our approach employs optimization with simulation
concepts using DTs (at sub-problem level) and a distributed coordination
algorithm to obtain a computationally efficient and near-optimal D-OPF
approach. Specifically, we optimize a low-fidelity linear model of the sub-
system and project the obtained solution to the high-fidelity simulation model
using DTs resulting in a near-optimal but feasible power flow solution.
Primarily, the use of DTs (at sub-problem level) excludes the necessity of
solving complex nonlinear OPFs and reduces the overall solution time for the
local subproblems. Then the distributed coordination algorithm is employed to
exchange the projected power flow variables with the neighboring distributed
areas. We carefully select shared variables based on the power flow in radial
power distribution systems that help us reduce the required number of macro-
iterations by an order of magnitude. Together these innovations result in
highly scalable and near-optimal algorithms for solving difficult D-OPF
problems. The proposed approach is also extensively evaluated using a
comprehensive cyber-power co-simulation environment.

##  II Network Modeling & Problem Formulation

In this paper, $(\cdot)^{T}$ represents matrix transpose; $(\cdot)^{*}$
represents the complex-conjugate; $(\cdot)^{(n)}$ represents the $n^{th}$
iteration; $\mathbb{Re}$ denotes the real part of the complex number; the
superscript $p$ (without parenthesis) denotes the three-phases, i.e.,
$\\{a,b,c\\}$ of the system.

###  II-A Power Flow Model

This section details the power flow models to formulate the OPF problems for
an unbalanced power distribution system. The unbalanced power flow equations
are based on our prior work [23]. Let us consider an unbalanced radial power
distribution network of $n$ buses where, $\mathcal{N}$ denotes the set of
buses in that system and $\mathcal{E}$ denotes the set of edges identifying
distribution lines that connect the ordered pair of buses $\\{ij\\}$,
$\forall$ $i,j\in\mathcal{N}$. Here, $\phi_{j}$ denotes the set of phases in
the bus $j$. Let $v_{j}^{p}=|V_{j}^{p}|^{2}$ be the squared magnitude of
voltage at bus $j\in\mathcal{N}$ for phase $p\in\phi_{j}$. Define
$\phi_{ij}=\\{pq:p\in\phi_{i}$ and $q\in\phi_{j}$
$\forall\\{ij\\}\in\mathcal{E}\\}$. Let
$l_{ij}^{pq}=(|I_{ij}^{p}||I_{ij}^{q}|)$ be the squared magnitude of the line
current flowing in the phase ${pq}\in\phi_{ij}$ of line $(i,j)$. Also,
$S^{pq}_{ij}=P^{pq}_{ij}+Q^{pq}_{ij}$ and
$z^{pq}_{ij}=r^{pq}_{ij}+jx^{pq}_{ij}$, where $pq\in\phi_{ij}$. Variable
$p_{L,j}^{p}$ and $q_{L,j}^{p}$ denote the active and reactive load
(respectively) connected at node $j$ of phase $p\in\phi_{j}$. Similarly,
subscript $D,j$ denotes the DER generation at node $j$; e.g., $p^{p}_{D,j}$
denotes the active power generation at node $j$ for phase $p$.
$\delta^{pq}_{ij}$ is the angle difference between the phase currents. First,
we define the nonlinear power flow model and then the linearized model is
detailed. For more details, please refer to [23].

####  II-A1 Nonlinear Model

With the approximated phase voltage and branch current angles, in [23]
developed the nonlinear power flow model for an unbalanced radial power
distribution system. The model is defined below in (1a). We model the loads in
the systems as constant power loads.

$\displaystyle
P_{ij}^{pp}-\sum_{q\in\phi_{j}}{l_{ij}^{pq}\left(r_{ij}^{pq}\cos(\delta_{ij}^{pq})-x_{ij}^{pq}\sin(\delta_{ij}^{pq})\right)}$
(1a) $\displaystyle=\sum_{k:j\rightarrow
k}P_{jk}^{pp}+p_{L,j}^{p}-p_{D,j}^{p}$ (1b) $\displaystyle
Q_{ij}^{pp}-\sum_{q\in\phi_{j}}{l_{ij}^{pq}\left(x_{ij}^{pq}\cos(\delta_{ij}^{pq})+r_{ij}^{pq}\sin(\delta_{ij}^{pq})\right)}$
$\displaystyle=\sum_{k:j\rightarrow
k}Q_{jk}^{pp}+q_{L,j}^{p}-q_{D,j}^{p}-q_{C,j}^{p}$ (1c) $\displaystyle
v_{j}^{p}=v_{i}^{p}-\sum_{q\in\phi_{j}}{2\mathbb{Re}\left[S_{ij}^{pq}(z_{ij}^{pq})^{*}\right]}+\sum_{q\in\phi_{j}}{z_{ij}^{pq}l_{ij}^{qq}}$
$\displaystyle+\sum_{q1,q2\in\phi_{j},q1\neq
q2}{2\mathbb{Re}\left[z_{ij}^{pq1}l_{ij}^{q1q2}\left(\angle(\delta_{ij}^{q1q2})\right)(z_{ij}^{pq2})^{*}\right]}$
(1d) $\displaystyle(P_{ij}^{pp})^{2}+(Q_{ij}^{pp})^{2}=v_{i}^{p}l_{ij}^{pp}$
(1e) $\displaystyle(l_{ij}^{pq})^{2}=l_{ij}^{pp}l_{ij}^{qq}$ (1f)

####  II-A2 Linear-approximated Model

The computational complexity augmented by the nonlinear power flow models to
the existing scalability issues associated with large-scale OPF problems can
be reduced using the approximated linearized power flow models. The linearized
power flow models can solve the power flow variables, such as nodal voltages,
within reasonable accuracy and reduce the computation time by the order of
magnitudes making the linear models scalable for large-scale power
distribution systems. A linear-approximated power flow model, three-phase
LinDistFlow model, is defined in (2a). The primary assumption is that the line
loses are negligible compared to the power flow in the system.

$\displaystyle P_{ij}^{pp}=\sum_{k:j\rightarrow
k}P_{jk}^{pp}+p_{L,j}^{p}-p_{D,j}^{p}$ (2a) $\displaystyle
Q_{ij}^{pp}=\sum_{k:j\rightarrow
k}Q_{jk}^{pp}+q_{L,j}^{p}-q_{D,j}^{p}-q_{C,j}^{p}$ (2b) $\displaystyle
v_{j}^{p}=v_{i}^{p}-\sum_{q\in\phi_{j}}{2\mathbb{Re}\left[S_{ij}^{pq}(z_{ij}^{pq})^{*}\right]}$
(2c)

###  II-B DER Modeling

We define the DER models for different network objectives. In this paper, to
achieve minimal cost, either reactive or active power generation of the DERs
is set as the decision variables of the OPF problems. The active power
generation, $p^{p}_{D,j}$ is assumed to measured and known, when the reactive
power generation, $q^{p}_{D,j}$ is the decision variable. If $S^{p}_{DR,j}$ is
defined as the nominal rating of the DER at node $j$ for phase $p\in\phi_{j}$,
then the DER is modeled as (3).

$-\sqrt{(S_{DR,j}^{p})^{2}-(p_{D,j}^{p})^{2}}\leq
q_{D,j}^{p}\leq\sqrt{(S_{DR,j}^{p})^{2}-(p_{D,j}^{p})^{2}}$ (3)

For network objectives, such as DER maximization, the active power generation,
$p^{p}_{D,j}$ is set as decision variables, and we assume $q^{p}_{D,j}=0$. The
DER model is defined in (4). However, if both active and reactive power
generation is considered as decision variables, then the DER is modeled as a
second-order conic equation as defined in (5).

$\small 0\leq p_{D,j}^{p}\leq S_{DR,j}^{p}$ (4)
$\small\left(p_{D,j}^{p}\right)^{2}+\left(q_{D,j}^{p}\right)^{2}\leq\left(S_{DR,j}^{p}\right)^{2}$
(5)

###  II-C Centralized OPF Problem

In this section, we define the centralized OPF problems for an unbalanced
radial distribution system, that can be solved using a central controller;
both nonlinear and linear power flow models have been used for such OPF
problems with minimization and maximization objectives.

####  II-C1 Loss Minimization

First, the active power loss minimization problem has been formulated where
the line loss is reduced by generating optimal $q^{p}_{D,j}$, without
violating operational limits. The optimization variable is
$X=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j},l^{pq}_{ij},q^{p}_{D,j}]^{T}$, $\forall
j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$. The OPF problem (C1) is
defined in (6a) using nonlinear power flow models. Here,
$\left(I^{rated}_{ij}\right)^{2}$ is the thermal limit for the line
$\\{ij\\}\in E$, and $V_{max}$ & $V_{min}$ denotes the nodal voltage limit.

$\displaystyle\text{{(C1)}}\hskip
11.38092pt\text{Minimize}\sum_{p\in\phi_{j},j:i\rightarrow
j}{l^{pp}_{ij}r^{pp}_{ij}}$ (6a) s.t.equation (1a) & (3) (6b) $\displaystyle
l_{ij}^{p}\leq\left(I^{rated}_{ij}\right)^{2}\text{\&}\hskip
8.5359pt{V_{min}}^{2}\leq v_{j}^{p}\leq{V_{max}}^{2}$ (6c)

To take advantage of the computational efficacy and scalability, next we
define the OPF problem with linear power flow models to minimize the line
losses. To further ease the computations, the nodal voltages are approximated
as $1.00$ pu, and thus, the line loss has been approximated by a convex cost
function $\sum\left(P^{pp}_{ij}\right)^{2}+\left(Q^{pp}_{ij}\right)^{2}$ for
all edges $\\{ij\\}\in\mathcal{E}$. Here for the OPF problem with linear power
flow model, the optimization variable is
$X=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j},q^{p}_{D,j}]^{T}$, $\forall
j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$. The OPF problem with
linear power flow model, (C2) is defined in (7a).

$\displaystyle\text{{(C2)}}\hskip
11.38092pt\text{Minimize}\sum_{p\in\phi_{j},j:i\rightarrow
j}{\left(P^{pp}_{ij}\right)^{2}}+{\left(Q^{pp}_{ij}\right)^{2}}$ (7a)
s.t.equation (2a) & (3) (7b) $\displaystyle{V_{min}}^{2}\leq
v_{j}^{p}\leq{V_{max}}^{2}$ (7c)

####  II-C2 DER Maximization

As a maximization problem, DER maximization objective has been formulated in
this paper. Specifically, the network objective is to maximize active power
generations, $p_{D,j}^{p}$, from all the DERs in all of the phases, without
exceeding any operational limits. For the DER maximization OPF, the
optimization variable is
$X=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j},l^{pq}_{ij},p^{p}_{D,j}]^{T}$, $\forall
j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$. The OPF problem (C3) is
defined in (8a) using the nonlinear power flow models. Similar to the previous
OPF problems, $\left(I^{rated}_{ij}\right)^{2}$ is the thermal limit for the
line $\\{ij\\}\in E$, and $V_{max}$ & $V_{min}$ denotes the nodal voltage
limit.

$\displaystyle\text{{(C3)}}\hskip 11.38092pt\text{Maximize}\sum_{\forall
j\in\mathcal{N}}\sum_{p\in\phi_{j}}{p^{p}_{D,j}}$ (8a) s.t.equation (1a) & (4)
(8b) $\displaystyle
l_{ij}^{p}\leq\left(I^{rated}_{ij}\right)^{2}\text{\&}\hskip
8.5359pt{V_{min}}^{2}\leq v_{j}^{p}\leq{V_{max}}^{2}$ (8c)

Similar to the problem (C2), the DER maximization objective problem can be
formulated as a linear OPF model with linear approximated power flow model
(eqn. (2a)). Since the objective function for this problem is already a linear
cost function, thus the whole OPF problem becomes linear optimization problem.
The optimization problem is detailed in (9a) and referred as problem (C4). The
optimization variable for (C4) is
$X=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j},p^{p}_{D,j}]^{T}$, $\forall
j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$.

$\displaystyle\text{{(C4)}}\hskip 11.38092pt\text{Maximize}\sum_{\forall
j\in\mathcal{N}}\sum_{p\in\phi_{j}}{p^{p}_{D,j}}$ (9a) s.t.equation (2a) & (4)
(9b) $\displaystyle{V_{min}}^{2}\leq v_{j}^{p}\leq{V_{max}}^{2}$ (9c)

##  III Distributed Algorithm

The optimization problem for radial power distribution networks– composed of
several connected areas, are naturally decomposable into multiple sub-problems
defined for the areas. In this section, first we detail the decomposition
approach for the OPF problems of three-phase unbalanced systems, and then the
local sub-problems are defined. Lastly, we briefly discuss the developed
cyber-physical co-simulation platform, where both the physical and the
communication layer are simulated to study the impact of the developed
distributed OPF algorithm on a real life distribution network.

###  III-A Decomposition Method & Network Approximations

Let the distribution system be composed of $N$ connected areas/microgrids, and
$\mathcal{A}=\\{A_{1},~{}A_{2},…,~{}A_{N}\\}$ be the set of all areas or
microgrids. Any area $A_{m}\in\mathcal{A}$ is defined as the directed radial
graph $A_{m}=\mathcal{G}(\mathcal{N}_{m},\mathcal{E}_{m})$; $\mathcal{N}_{m}$
and $\mathcal{E}_{m}$ defines the set of nodes and lines in the system. The
set of central optimization variables is defined as
$X=\bigcup_{m=1}^{N}X_{m}$, where $X_{m}$ represents the set of all local
optimization variables for Area $A_{m}$. Similarly
$\mathcal{S}=\bigcup_{m=1}^{N}\mathcal{S}_{m}$ represent the set of
constraints for overall problem, and $\mathcal{S}_{m}$ represents the local
constraints. If central objective $f(X)$ is a decomposable cost function, then
the problem can be decomposed into several local sub-problems and written as
(10). Here, $f_{m}$ is defined as the local objective function for Area
$A_{m}$. The complicating variables, $y_{{m}^{^{\prime}}}$ are shared by the
neighboring areas of $A_{m}$, and are kept constant for solving the local
optimization problem in this primal decomposition-based method.

$\small\min_{X\in\mathcal{S}}{\hskip 2.84544ptf(X)}={\hskip
2.84544pt}{\sum_{m=1}^{N}\min_{X_{m}\in\mathcal{S}_{m}}f_{m}(X_{m},y_{{m}^{^{\prime}}})}$
(10)

As the overall system topology was radial to begin with, the decomposed areas
of the networks are also connected radially and have unique upstream &
downstream sections. This specific structure of the system helps to identify
the unique parent area and the child areas for any area $A_{m}$, which in
turns associates the complicating/shared variables, $y_{{m}^{^{\prime}}}$ for
the local sub-problems – exchanged among the neighboring computing agents to
solve the overall problem. For the proposed decomposition approach, the
complicating variables are the voltages and the power flows at the shared bus.
Let for any Area $A_{m}$, bus $0$ is shared with its parent/Upstream area (UA)
and node $k$ ($\\{jk\\}\in\mathcal{E}_{m}$) is shared with its
child/Downstream Area (DA). Then the complicating variable,
$y_{{m}^{^{\prime}}}=[y^{ua},y^{da}]$, where $y^{ua}=[v_{0_{ua}}]$ , and
$y^{da}=[P_{jk_{da}},Q_{jk_{da}}]$. Here, each area solves related local sub-
problems in parallel by assuming the constant complicating variable, i.e., a
fixed voltage at the shared bus with the unique UA and constant loads at the
shared buses with DAs. In terms of power system’s physics, this approximates
the whole upstream and downstream network segments of an area for radial power
distribution systems.

After solving local sub-problems, the respective complicating variables, i.e.,
the total power requirements in that area and the shared (with DA) bus
voltages area are shared with UA and DAs, respectively. Then, the sub-problems
are solved again with updated shared variables from the neighbors. When a
consensus is attained at the shared boundaries, then the decision variables
are dispatched within the area. The consensus at the boundary can be achieved
using Fixed Point Iteration methods as described in [24].

Figure 1: Use of Digital Twin Model in Distributed Optimization for Power
System

###  III-B Local Distributed OPF Model

In this section, the decomposed local sub-problems at $n^{th}$ macro-
iteration, (D1)-(D4), of the corresponding central OPFs, (C1)-(C4) are
defined. For any Area $A_{m}$, let bus $0$ is shared with its UA and node $k$
($\\{jk\\}\in\mathcal{E}_{m}$) is shared with its DAs. Also, the subscript of
a node represents the area, that has solved the variable; e.g., $v_{0_{m}}$,
$p_{k_{m}}$ are solved by area $A_{m}$, where $v_{0_{ua}}$, $P_{jk_{da}}$ are
solved by UA and DA, respectively. In the proposed model, each local computing
agents solves the corresponding distributed sub-problems, and over macro-
iterations, convergence is achieved.

####  III-B1 Loss Minimization

The OPF problem (C1)-(C2) represent the central formulation of the loss
minimization problem using nonlinear and approximated models, respectively,
where the active power line losses are minimized by controlling the reactive
power generations from the DERs. The corresponding decomposed sub-problems for
Area $A_{m}$ , $\forall\\{ij\\}\in\mathcal{E}_{m}$ and $\forall
j\in\mathcal{N}_{m}$ are defined in (11a)-(12a), respectively. Constraint
(11c) and (12c) represents the approximated UA and DA with the help of
complicating/shared variables, as described in the previous section.

$\displaystyle\text{{(D1)}}\hskip 11.38092pt\min\hskip
5.69046ptf_{m}(X_{m}^{(n)})=\sum_{p\in\phi_{j},j:i\rightarrow
j}{l^{pp}_{ij}r^{pp}_{ij}}$ (11a) s.t. equation (6b) & (6c) (11b)
$\displaystyle v_{0_{m}}^{p\hskip 1.0pt(n)}=v_{0_{ua}}^{p\hskip
1.0pt(n-1)};\hskip 2.84544ptp_{k_{m}}^{p\hskip 1.0pt(n)}=P_{jk_{da}}^{pp\hskip
1.0pt(n-1)};\hskip 2.84544ptq_{k_{m}}^{p\hskip 1.0pt(n)}=Q_{jk_{da}}^{pp\hskip
1.0pt(n-1)}$ (11c)

$\displaystyle\text{{(D2)}}\hskip 11.38092pt\min\hskip
5.69046ptf_{m}(X_{m}^{(n)})=\sum_{p\in\phi_{j},j:i\rightarrow
j}{\left(P^{pp}_{ij}\right)^{2}}+{\left(Q^{pp}_{ij}\right)^{2}}$ (12a) s.t.
equation (7b) & (7c) (12b) $\displaystyle v_{0_{m}}^{p\hskip
1.0pt(n)}=v_{0_{ua}}^{p\hskip 1.0pt(n-1)};\hskip 2.84544ptp_{k_{m}}^{p\hskip
1.0pt(n)}=P_{jk_{da}}^{pp\hskip 1.0pt(n-1)};\hskip 2.84544ptq_{k_{m}}^{p\hskip
1.0pt(n)}=Q_{jk_{da}}^{pp\hskip 1.0pt(n-1)}$ (12c)

####  III-B2 DER Maximization

The OPF problem (C3)-(C4) represent the central OPF formulation of the DER
generation maximization problem using nonlinear and approximated models,
respectively, where the real power generations for all the phases are
maximized within the physical limits and without violating the operational
constraints. Similar to loss minimization objective, the corresponding
decomposed sub-problems for Area $A_{m}$ are defined in (13a)-(14a),
respectively. Similar to loss minimization problem, constraint (13c) and (14c)
represents the approximated UA and DA with the help of complicating/shared
variables.

$\displaystyle\text{{(D3)}}\hskip 11.38092pt\max\hskip
5.69046ptf_{m}(X_{m}^{(n)})=\sum_{\forall
j\in\mathcal{N}_{m}}\sum_{p\in\phi_{j}}{p^{p}_{D,j}}$ (13a) s.t. equation (8b)
& (8c) (13b) $\displaystyle v_{0_{m}}^{p\hskip 1.0pt(n)}=v_{0_{ua}}^{p\hskip
1.0pt(n-1)};\hskip 2.84544ptp_{k_{m}}^{p\hskip 1.0pt(n)}=P_{jk_{da}}^{pp\hskip
1.0pt(n-1)};\hskip 2.84544ptq_{k_{m}}^{p\hskip 1.0pt(n)}=Q_{jk_{da}}^{pp\hskip
1.0pt(n-1)}$ (13c)

$\displaystyle\text{{(D4)}}\hskip 11.38092pt\max\hskip
5.69046ptf_{m}(X_{m}^{(n)})=\sum_{\forall
j\in\mathcal{N}_{m}}\sum_{p\in\phi_{j}}{p^{p}_{D,j}}$ (14a) s.t. equation (9b)
& (9c) (14b) $\displaystyle v_{0_{m}}^{p\hskip 1.0pt(n)}=v_{0_{ua}}^{p\hskip
1.0pt(n-1)};\hskip 2.84544ptp_{k_{m}}^{p\hskip 1.0pt(n)}=P_{jk_{da}}^{pp\hskip
1.0pt(n-1)};\hskip 2.84544ptq_{k_{m}}^{p\hskip 1.0pt(n)}=Q_{jk_{da}}^{pp\hskip
1.0pt(n-1)}$ (14c)

###  III-C Simulation Within Optimization: Use of Digital Twin

To reduce the compute burdens to solve the NLP optimization problem, a method
utilizing simulation within optimization has been proposed by incorporating
the Digital Twin (DT) of the system. DTs are the detailed virtual model of the
real-world system/network, linked to their physical counterpart by the
continuous, real-time information flow. Specific to the power distribution
systems, DTs are the software-based abstractions of complex physical networks
of power distribution systems which are connected via a communication link to
the actual physical system; real-time analysis can be conducted on these
virtual models/DTs for decision support purposes. By using the DT models for
solving optimization problems, we can examine the optimal solutions ahead of
time, and without enacting the decision variables directly in the real
systems, the feasibility of the systems can be ensured – preventing the
potential unwanted events for the network. This unfolds the opportunity to use
approximated linear power flow models for solving OPFs – defined in equation
(2a), instead of using computationally complex nonlinear network models as
described by equation (1a).  Upon implementing the optimal decision variables
in the DT models (e.g., optimal DER generations – obtained from the
approximated models), the corresponding nonlinear solutions of the state
variables of the system, such as nodal voltages, branch power flows etc., of
the system can be attained.

Mathematically, the optimization variable $X$ for the OPF problems is composed
of state variables, $w$, and the optimization decision variables, $z$, i.e.,
$X=[w,z]^{T}$. For example, in the (C1) OPF problem,
$w=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j},l^{pq}_{ij}]$, and $z=q^{p}_{D,j}$
$\forall j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$. Similarly for
approximated (C4) OPF problem, $w=[P^{pp}_{ij},Q^{pp}_{ij},v^{p}_{j}]$, and
$z=p^{p}_{D,j}$ $\forall j\in\mathcal{N}$ & $\forall\\{ij\\}\in\mathcal{E}$.
We use DT models for the simulation within optimization as such: after solving
the approximated OPF problem (equation (15)), the optimal decision variables,
$z_{Lin}^{*}$ are implemented in the DT model – defined by the function
$\Phi_{DT}$, and then the nonlinear solutions are extracted. Thus, using (16),
the corresponding nonlinear state variables of the system, $w_{Nl}^{*}$ are
obtained upon effectuating the linear optimizer $z_{Lin}^{*}$ in the system;
the optimizer $X^{*}$ attained from this simulation within optimization
approach is shown in equation (17). Any nonlinear power-flow model of the
system can be used to represent the functional form for the DT.

$X_{Lin}^{*}=[w_{Lin}^{*},z_{Lin}^{*}]:=\underset{X_{Lin}\in\mathcal{S}}{\text{argmin}}{\hskip
2.84544ptf(X_{Lin})}$ (15) $w_{Nl}^{*}=\Phi_{DT}(z_{Lin}^{*})$ (16)
$X^{*}=[w^{*}_{Nl},z_{Lin}^{*}]^{T}$ (17)

Similarly, in the distributed decision-making process, the use of DTs for
solving the OPF problems helps to build a routine that can leverage the
computational speed and scalability of the approximated optimization problems
with linear power flow models. Each distributed agent has the DT models
($\Phi_{DT}$) of the corresponding local networks. The agents first solve the
local approximated sub-problems (similar to (15)), and then extract the
nonlinear solutions by using (16) for the local system. During the macro-
iteration, i.e., communication among these distributed agents, appropriate
variables from nonlinear solution, $w_{Nl}^{*}$ are exchanged, instead of
using the solutions from linear models, i.e., $w_{Lin}^{*}$. Specifically,
alternate to the direct exchange of linear solutions, the nonlinear solutions
from DT simulations are shared among the distributed agents. This effectively
helps to ensure the feasibility of the optimal solution of the overall problem
at the end of the macro-iterations, even though the computationally simple and
scalable approximated models are being solved.  The block diagram of the
proposed distributed algorithms with simulation within optimization using DT
models is shown in Fig. 1

###  III-D Co-simulation Platform

A realistic evaluation of distributed algorithms require modeling the
associated communication system. To this end, we developed a multi-agent
cyber-physical co-simulation package using HELICS (Heirarchical Engine for
Large-scale Infrastructure Co-Simulation), a generic co-simulation platform
[25]. Our cyber-physical co-simulation package is composed of the following
software: Python 3, GridLAB-D [26], NS-3 [27], which are used to simulate the
computing agents, the power distribution networks, the communication networks,
respectively. HELICS [25] is then employed to synchronize the simulations and
facilitates message passing among the three simulation software (see Fig. 2).
Each area has a controller agent shown at the top which has a single
connection, via HELICS, to the communication network simulated in NS-3. In
this way, the controllers are able to communicate with each other and with the
device agents in the local area. The device agents each have a connection, via
HELICS, to the physical device simulated in GridLAB-D and to the communication
network. The purpose of the device agent is twofold. It formats the raw data
it receives from GridLAB-D so the control agent can understand it and it
parses commands received from the control agent to the GridLAB-D. The control
agents are capable of running their own GridLAB-D simulations using local area
models as digital twins (DT). Each area employs the respective DT to obtain
the actual power flow variables the system observes upon implementing the
approximated OPF control variables.

Figure 2: Co-Simulation Structure for Cyber-Physical System

##  IV Results and Discussions

In this section, a detailed evaluation of the proposed D-OPF algorithm for a
3-phase unbalanced distribution system is presented. First, we have validated
our distributed OPF approach for nonlinear optimization problems by comparing
the solution with the central formulation. Then, we have simulated and
verified the simulation within optimization approach that has been proposed in
this paper with cyber-physical co-simulation platform and uses approximated
OPF problems. Lastly, a time-series simulation has been performed to showcase
the effect of network optimization using the proposed method, and different
stressed communication network topologies have been tested to demonstrate the
efficacy and robustness of the proposed distributed OPF algorithm.

Figure 3: Sectioned IEEE 123 Bus System

###  IV-A Simulated System

To test and validate our proposed method, we have simulated two network
objective problems with three-phase unbalanced IEEE-123 bus test system that
has been assumed to be composed of 4 areas. Three different scenarios of DERs
are simulated in this paper to test the proposed method – (i) 10 three-phase
DERs (at node 15, 23, 30, 37, 49, 50, 51, 67, 78, 107) with capacity of 60 kVA
per phase, (ii) 30 three-phase DERs with capacity of 48 kVA per phase, and
(iii) 30 three-phase DERs with capacity of 60 kVA per phase (Fig. 3). Please
note, scenario (ii) is simulated for loss minimization problem to provide a
more challenging problem for the algorithm in terms of required communication
properties. Further, for the time series simulation, the loads are also varied
randomly between 78% to 98% of the nominal values for 60 minutes. In addition
to that, as a digital twin function, $\Phi_{DT}$, GridLab-D power flow
simulator has been used in this study.

(a) Convergence at Boundary

(b) Node Voltages: Loss Minimization

(c) Node Voltages: DER Maximization

Figure 4: D-OPF Results for Loss Minimization and DER Maximization Objectives

###  IV-B Validation of Decomposition Method

To validate the proposed D-OPF method for unbalanced power distribution
systems, we have simulated two optimization problems, (D1) & (D3) for DER
scenarios (i) & (iii), respectively, and compared with the solutions of (C1) &
(C3) (see Table I). It is evident that the optimal solutions of D-OPF match
with the corresponding Central OPF (C-OPF) solutions for both objectives. For
example, the objective value for the DER maximization problem, i.e., total DER
generation is 5.08 MW & 5.18 MW, when solved using D-OPF & C-OPF method,
respectively. However, D-OPF takes significantly less time to solve, compared
to the C-OPF method – takes almost 10 times or more for some cases. From Fig.
4(a), it can be seen that it only takes 4 macro-iterations to reach the global
solution. In addition to that, the nodal voltages of D-OPF and C-OPF solutions
are compared (Fig. 4(b)-4(c)). Later the nodal voltages of D-OPFs are compared
with the GridLab-D/OpenDSS model – upon enacting the decision variables from
D-OPF, and the voltage difference is in the order of $10^{-4}$ or less (Fig.
5). These comparisons validate the feasibility and the optimality of the
proposed D-OPF solutions for both minimization and maximization objectives.

(a) Loss Minimization

(b) DER Maximization

Figure 5: Validation for Decomposition Approach TABLE I: Validation of D-OPF method OPF Problem | Objective | Time (s)
---|---|---
D-OPF | C-OPF | D-OPF | C-OPF
Loss Min.: Scenario (i) | 43.976 kW | 43.975 kW | 20 | 520
DER Max.: Scenario (iii) | 5.08 MW | 5.18 MW | 20 | 200

###  IV-C Co-simulation & Use of DT

Next, to showcase the robustness and the practicability of the proposed
algorithm, the D-OPF method is implemented on a cyber-physical co-simulation
platform. For this case, distributed problems (D2) & (D4) are solved to
leverage the simplicity of the linear network model for edge computing use-
cases to solve the OPFs. In this paper, the DTs of the local areas are modeled
using GridLab-D power flow simulator. For this paper, it is assumed that the
distributed agents/local controllers communicate with neighbor every
2-seconds. And also, after the convergence is attained, it takes around $0.2$
seconds to dispatch the decision variables in the system.

TABLE II: Loss minimization Result Comparison Method | Loss (kW) | Communication-round
---|---|---
DER Scenario | Scenario (i) | Scenario (ii) | Scenario (i) | Scenario (ii)
No OPF | 63.132 | 53.338 | 0 | 0
Linear C-OPF | 44.875 | 26.508 | 0 | 0
Linear D-OPF | 44.959 | 26.546 | 5 | 5
Linear D-OPF+DT | 44.959 | 27.238 | 5 | 5
Nonlinear C-OPF | 43.975 | 27.160 | 0 | 0

####  IV-C1 Validation

Now we validate the co-simulation platform using the linear approximate models
by comparing the resulting objective values. Table II and Table III, show the
comparisons between D-OPF and C-OPF results for different cases. Table II
shows the loss minimization case for DER scenarios (i) and (ii), where real
power generations are fixed at $50$ and $20$ kW per phase, respectively. For
both DER scenarios, line losses matches with the nonlinear solutions. For
example, in scenario (i), the active power loss is $44.959$ kW for linear
model with the DT model, and the nonlinear solution is $43.975$ kW. However,
for scenario (ii), it can be seen that the use of DT model has be potential to
get closer solutions to the nonlinear solution. In Table II, the line loss for
scenario (ii) is $26.546$ kW when DT models are not being used, however, if
the DT models are used, the objective value is $27.23$ kW – that is closer to
the nonlinear solution $27.16$ kW. The use of DT with approximated model does
not necessarily improve the results, but it guarantees the feasibility of the
local solutions, before exchanging the boundary variables with the neighbors.
Also, it takes around 5 communication-rounds among the agents to reach the
consensus at the boundary. Similarly, Table III validates the algorithm on the
co-simulation platform for DER maximization problem that uses DER scenario
(iii). When no OPF is performed, there were voltage limit violation in the
system. However, if operational limits are considered, the optimal generation
is $5.135$, $5.163$, and $5.18$ MW for linear D-OPF, linear D-OPF+DT, and
nonlinear C-OPF, respectively. Thus the use of the DT model for DER
maximization problem also brings the solution closer to the actual nonlinear
solutions. Also, using DT model reduces the number of communication rounds
from 7 to 5. Thus, Table II and Table III validate the co-simulation platform
that uses computationally simple linear approximated OPF model with the DT
model of the system.

TABLE III: DER Maximization Results Comparison Method | Generation (MW) | Communication-round
---|---|---
DER Scenario | Scenario (iii) | Scenario (iii)
No OPF | 5.400 | 0
Linear C-OPF | 5.135 | 0
Linear D-OPF | 5.135 | 7
Linear D-OPF+ DT | 5.163 | 5
Nonlinear C-OPF | 5.180 | 0

####  IV-C2 Time-series Simulation

In Fig. 6, we show how the voltages in the system are impacted by the loss
minimization (scenario (ii)) and DER maximization (scenario (iii)) objectives
respectively. In the loss-minimization case, we expect that all of the bus
voltages will be shifted towards the upper limit. In the DER maximization case
we expect that bus voltages that exceed the voltage limit will be reduced to
be within the voltage limits. These expectations are confirmed in Fig. 6.

(a) Loss Minimization

(b) DER Maximization

Figure 6: Voltage Distributions

####  IV-C3 Communication Network Stress Test

Next we test the robustness of the proposed algorithm, when subjected to poor
communication network conditions; it is tested with two different
communication network topologies (ideal and ring topology) for loss
minimization problem (D2) with scenario (ii) and for DER maximization problem
(D4) with scenario (iii). In the ideal topology, each area controller has a
direct link with its neighbors and with devices in its own area. In the ring
topology all devices and controllers are connected in a single large loop
without regard to physical location. Each link between network nodes is a
point-to-point link with a negligible delay and has different bandwidths
(3kbps, 2kbps, or 1kbps) to show the worst-case performance of the algorithm
(see Fig. 7). The lines marked with circles indicate results run with the ring
topology and lines marked with a star indicate the ideal topology. The blue
line shows the results if no communication is allowed between controllers and
the controllers have perfect communication links with local devices. There are
several effects of poor communications that we have observed:

  * •

If a controller doesn’t receive data from either UA or DA because of data
delays, it will falsely assume convergence and dispatch the inverters
prematurely. When the delayed data does arrive it will continue iterating
again.

  * •

Delayed communication between areas does not prevent convergence except in
extreme cases.

  * •

The quality of the resulting solution as measured after convergence is not
impacted by the communication delays.

(a) Loss Minimization

(b) DER Maximization

Figure 7: Communication Stress Test

####  IV-C4 Advantages of Using DT with LP

The LP model can estimate the nodal voltages with moderate accuracy; however,
for corner cases, such as when the system is operating at the voltage limits,
the approximation error can lead to a higher optimality gap and even
infeasible solutions. In Fig. 8(a), we have shown a snippet of voltage at
phase C for a high loading condition; using a DT model yields more accurate
nonlinear nodal voltages than LP solutions that approximate higher voltages.
If the voltage lower bound constraint is active, an overestimation of the
voltage can lead to infeasible space upon implementing the LP solutions to the
actual system. However, by using the DT model, we can estimate the nonlinear
voltages more accurately and can preemptively tackle any infeasibility. Again
for the DER maximization case, when the voltage upper bound is active, the
higher approximated voltage can lead us to sub-optimal DER generation
solutions. Here, a low loading condition scenario has been simulated, and a
snippet of phase A voltages are shown in Fig. 8(b). Upon implementing the LP
solution in the OpenDSS, we obtain lower feeder voltages indicating that more
DER can be integrated before the actual feeder limit is reached. Thus, LP
model results in sub-optimal solutions. However, the LP+DT model does not
overestimate the voltages and hence optimized for the maximum generation that
can be integrated without causing any voltage violations. We observed that the
optimal solutions for the two models differ by $300$ kW. This case further
demonstrates the usability of the DT model in accurately solving the OPF
problems without augmenting any complicated nonlinear formulation, but still
resulting in more accurate optimal solutions.

(a) Loss Minimization

(b) DER Maximization

Figure 8: Voltage Approximation Error

##  V Conclusions

In this paper, a novel distributed optimal power flow method for radial three-
phase, unbalanced power distribution system has been proposed. The proposed
approach uses simulation within optimization to improve the scalability and
reduce solution time for each local subproblem. We employ detailed DT models
of the local network in conjugation with a linear-approximated OPF to ensure a
power-flow feasible and near-optimal solution for multiple OPF problems for
power distribution system. Additionally, the proposed distributed coordination
method achieves a converged network-level optimal solution using significantly
fewer communication rounds. We also evaluated the performance of the proposed
algorithm using a cyber-power co-simulation environment with various
communication network parameters, to validate the robustness of the proposed
distributed OPF algorithm under stressed communication. For various cases, it
is shown that the proposed method successfully converges to the network-level
optimal solution within reasonable time, even when the communication network
is facing significant stress, such as, communication delays and a low
bandwidth communication infrastructure.

## 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] 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.
  * [4] B. H. Kim and R. Baldick, “A comparison of distributed optimal power flow algorithms,” IEEE Transactions on Power Systems, vol. 15, no. 2, pp. 599–604, 2000.
  * [5] 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.
  * [6] 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\.
  * [7] G. Qu and N. Li, “Harnessing smoothness to accelerate distributed optimization,” IEEE Transactions on Control of Network Systems, vol. 5, no. 3, pp. 1245–1260, 2017.
  * [8] S. Mhanna, G. Verbič, and A. C. Chapman, “Adaptive admm for distributed ac optimal power flow,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 2025–2035, 2018.
  * [9] 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.
  * [10] J. April, M. Better, F. Glover, and J. Kelly, “New advances and applications for marrying simulation and optimization,” in Proceedings of the 2004 Winter Simulation Conference, 2004., vol. 1, IEEE, 2004.
  * [11] H. Nakayama, M. Arakawa, and R. Sasaki, “Simulation-based optimization using computational intelligence,” Optimization and Engineering, vol. 3, no. 2, pp. 201–214, 2002.
  * [12] A.-T. Nguyen, S. Reiter, and P. Rigo, “A review on simulation-based optimization methods applied to building performance analysis,” Applied energy, vol. 113, pp. 1043–1058, 2014.
  * [13] M. Zhou, J. Yan, and D. Feng, “Digital twin framework and its application to power grid online analysis,” CSEE Journal of Power and Energy Systems, vol. 5, no. 3, pp. 391–398, 2019.
  * [14] H. Pan, Z. Dou, Y. Cai, W. Li, X. Lei, and D. Han, “Digital twin and its application in power system,” in 2020 5th International Conference on Power and Renewable Energy (ICPRE), pp. 21–26, IEEE, 2020.
  * [15] C. Brosinsky, D. Westermann, and R. Krebs, “Recent and prospective developments in power system control centers: Adapting the digital twin technology for application in power system control centers,” in 2018 IEEE International Energy Conference, pp. 1–6, IEEE, 2018.
  * [16] X. He, Q. Ai, R. C. Qiu, and D. Zhang, “Preliminary exploration on digital twin for power systems: Challenges, framework, and applications,” arXiv preprint arXiv:1909.06977, 2019.
  * [17] R. Darbali-Zamora, J. Johnson, A. Summers, C. B. Jones, C. Hansen, and C. Showalter, “State estimation-based distributed energy resource optimization for distribution voltage regulation in telemetry-sparse environments using a real-time digital twin,” Energies, vol. 14, no. 3, p. 774, 2021.
  * [18] T. Erseghe, “Distributed optimal power flow using admm,” IEEE transactions on power systems, vol. 29, no. 5, pp. 2370–2380, 2014.
  * [19] 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.
  * [20] 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.
  * [21] S. Magnússon, G. Qu, and N. Li, “Distributed optimal voltage control with asynchronous and delayed communication,” IEEE Transactions on Smart Grid, 2020.
  * [22] R. Sadnan and A. Dubey, “Distributed optimization using reduced network equivalents for radial power distribution systems,” IEEE Transactions on Power Systems, vol. 36, no. 4, pp. 3645–3656, 2021.
  * [23] 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\.
  * [24] R. Sadnan, T. Asaki, and A. Dubey, “Online distributed optimization in radial power distribution systems: Closed-form expressions,” in 2021 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm), pp. 51–56, IEEE, 2021.
  * [25] B. Palmintier, D. Krishnamurthy, P. Top, S. Smith, J. Daily, and J. Fuller, “Design of the helics high-performance transmission-distribution-communication-market co-simulation framework,” in 2017 Workshop on Modeling and Simulation of Cyber-Physical Energy Systems (MSCPES), pp. 1–6, IEEE, 2017.
  * [26] D. P. Chassin, K. P. Schneider, and C. Gerkensmeyer, “GridLAB-D: An open-source power systems modeling and simulation environment,” in IEEE/PES Transmission and Distribution Conference and Exposition, pp. 1–5, IEEE, May 2008.
  * [27] NS-3 Network Simulator.  https://www.nsnam.org/.


# Collision-free Formation Control of Multiple Nano-quadrotors

Anh Tung Nguyen1, Ji-Won Lee2, Thanh Binh Nguyen3, and Sung Kyung Hong1 *This
work was supported by the MSIT (Ministry of Science and ICT), Korea, under the
ITRC (Information Technology Research Center) support program
(IITP-2021-2018-0-01424) supervised by the IITP (Institute for Information &
communications Technology Promotion). (Corresponding author: Sung Kyung
Hong.)1Anh Tung Nguyen and Sung Kyung Hong are with Faculty of Mechanical and
Aerospace Engineering, Sejong University, Seoul 143–747(05006), Korea
<EMAIL_ADDRESS><EMAIL_ADDRESS>Lee is with Human-Robot
Interaction Research Center, Korea Institute of Robotics and Technology
Convergence, Pohang 37553, Korea<EMAIL_ADDRESS>Binh Nguyen is with
Department of Control Engineering and Automation, Thuyloi University, 175 Tay
Son, Dong Da, Hanoi, Vietnam<EMAIL_ADDRESS>
###### Abstract

The utilisation of unmanned aerial vehicles has witnessed significant growth
in real-world applications including surveillance tasks, military missions,
and transportation deliveries. This letter investigates practical problems of
formation control for multiple nano-quadrotor systems. To be more specific,
the first aim of this work is to develop a theoretical framework for the time-
varying formation flight of the multi-quadrotor system regarding anti-
collisions. In order to achieve this goal, the finite cut-off potential
function is devoted to avoiding collisions among vehicles in the group as well
as between vehicles and an obstacle. The control algorithm navigates the group
of nano-quadrotors to asymptotically reach an anticipated time-varying
formation. The second aim is to implement the proposed algorithm on Crazyflies
nano-quadrotors, one of the most ubiquitous indoor experimentation platforms.
Several practical scenarios are conducted to tendentiously expose anti-
collision abilities among group members as well as between vehicles and an
obstacle. The experimental outcomes validate the effectiveness of the proposed
method in the formation tracking and the collision avoidance of multiple nano-
quadrotors.

###### Index Terms:

Multi-Robot Systems, Formation Control, Collision Avoidance, Obstacle
Avoidance, Nano-quadrotors.

##  I Introduction

Over the past few decades, multi-robot systems have gained massive popularity
in industrial societies due to their ability to describe large-scale
interconnected systems in a variety of real-world applications such as
transportation systems [1, 2] and power systems [3]. In the multi-robot
systems, the formation control problem aims at steering multiple robots in a
network to achieve and maintain their predefined geometric patterns in their
states, posing an immense challenge to the scientific community. Since a
colossal number of versatile robotic applications are developed, there have
been a great deal of the works devoted to unmanned aerial vehicles (UAVs) [4],
ground mobile robots [5], and especially to formation control of mobile robots
[6].

Quadrotors, one of the most ubiquitous classes of UAVs, have been intensively
developed in many broad applications for assisting humans in difficult
missions or hazard environments [7, 8, 9, 10], e.g., in agriculture [11],
industry [12], and military [13]. Among numerous commercial products of
quadrotors, this letter mainly focuses on studying nano-quadrotors (a
miniature dimension of quadrotors). The use of nano-quadrotors tolerates a
convenience and simply installed experiments to verify control algorithms as
well as conducting new research [14, 15]. In addition, actual flight tests on
the nano-quadrotors also can be deployed in incommodious applications, and
suitable for most laboratory setups. After successfully validating control
algorithms on the nano-quadrotors, researchers possibly scale up to other
sizeable quadrotors in dealing with a colossal number of civilian
applications.

Recently, collisions among robots have become a vital issue when autonomously
operating multi-robot systems, i.e., collisions among group members, and
collisions between robots and obstacles. Based on actual applications [16, 17,
18, 19, 20, 21], the desired spots of robots are generally marked with a view
to guaranteeing a safe distance among group members. Because of this
arrangement, collisions among vehicles possibly occur when they move from
their initial positions to the desired spots to result in an anticipated
formation. Additionally, after completing the given formation, the group of
robots may track a predefined trajectory. This task is probably unsuccessful
when obstacles appear in the trajectory and cause collisions with the robots.
With the purpose of dealing with the collision avoidance problem, path
planning based on optimisation problems was presented in [16, 22, 17]. Due to
the computational cost of solving the optimisation problems, outcomes of such
existing studies were limited by numerical results or operating a vehicle. A
coverage planning for ground mobile robots was introduced in [23] to solve a
problem of obstacle-cluttered environments in the two-dimensional space. The
potential field approach was developed to tackle a collision issue in [24, 18,
21, 25]. The authors in [18, 19] showed outdoor experimental results of
formation control of multiple UAVs with consideration to inter-vehicle
collisions. Due to unclear collision points, the demonstrations of their
algorithms would be vague.

Motivated by the above observations, this letter is concerned with the
formation control problem of multiple nano-quadrotor systems. With
consideration to the aforementioned collision issues, a collision-free
formation control algorithm is proposed based on potential functions.
Accordingly, our approach develops a new control Lyapunov function by which
the multiple nano-quadrotor systems are stable and the formation tracking
errors exponentially converge. In summary, our main contributions can be
highlighted as follows.

  1. i)

This paper addresses a collision-free problem in a group of nano-quadrotors as
well as between nano-quadrotors and an obstacle. A novel potential function is
proposed to avoid local minima phenomenons. In addition, this function also
guarantees the smoothness of the control input that increases practical
capabilities.

  2. ii)

Based on the control Lyapunov function, a collision-free formation protocol is
presented. As a result, the proposed protocol ensures that the group of nano-
quadrotors asymptotically tracks the anticipated formation trajectory with no
collisions.

  3. iii)

It is worth nothing that the scope of our study mainly focuses on indoor
applications. To the best of our knowledge, due to the complexity of anti-
collision algorithms, comparatively little experimental results of controlling
multi-aerial vehicles have been published. The efficacy of the proposed method
is validated by an indoor experimental scenario. Crazyflies, one of the most
ubiquitous indoor experimentation platforms, are employed to perform actual
experiments. Experimental outcomes give us a powerful demonstration of the
presented method.

Notation:  The notations $X\geq Y$ and $X>Y$ mean that $X-Y$ is positive semi-
definite and positive definite, respectively. $\otimes$ stands for Kronecker
product of two arbitrary-size matrices. ${\mathbb{R}}_{+}$ denotes the set of
real positive numbers; ${\mathbb{R}}^{n}$ and ${\mathbb{R}}^{n\times m}$ stand
for sets of real $n$-dimensional vectors and $n$-row $m$-column matrices,
respectively; $I_{n}$ is the $n\times n$ identity matrix; and
$\lambda_{\min}\left\\{W\right\\}$
$\left(\lambda_{\max}\left\\{W\right\\}\right)$ represents the minimum
(maximum) eigenvalue of matrix $W$. Next, for a vector $x\in{\mathbb{R}}^{n}$,
$\|x\|_{2}=\sqrt{x^{T}x}$ stands for Euclidean norm in ${\mathbb{R}}^{n}$.

##  II Preliminaries

Figure 1: Crazyflie Nano-quadrotor model.

###  II-A Algebraic graph theory

Let ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathcal{A}})$ be a weighted
digraph with the set of vertices ${\mathcal{V}}=\\{1,2,...,N\\}$, the set of
edges ${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$, and the
weighted adjacency matrix ${\mathcal{A}}=[a_{ij}]_{{i,j\in{\mathcal{E}}}}$.
For any $(i,j)\in{\mathcal{E}},~{}i\neq j$, the element of the weight
adjacency matrix $a_{ij}$ is positive if vertices $i$-th and $j$-th can
communicate with each other, while $a_{ij}=0$ in the cases of
$(i,j)\notin{\mathcal{E}}$ or $i=j$. The degree of a vertex $i$-th is denoted
as $deg_{i}^{in}=\sum_{j=1}^{n}a_{ij}$, and the degree matrix of the graph
${\mathcal{G}}$ is defined as
${\mathcal{D}}=\boldsymbol{diag}\big{(}deg_{1}^{in},deg_{2}^{in},\dots,deg_{N}^{in}\big{)}$.
The Laplacian matrix is defined as
${\mathcal{L}}=[\ell_{ij}]_{{i,j\in{\mathcal{E}}}}={\mathcal{D}}-{\mathcal{A}}$.
Further, ${\mathcal{G}}$ is called an undirected graph if and only if
${\mathcal{A}}$ is a symmetric matrix. An edge of the undirected graph
${\mathcal{G}}$ is denoted by an unordered pair $(i,j)\in{\mathcal{E}}$. The
undirected graph is strongly connected if for any pair of vertices, there
exists at least a path between two vertices. The set of all neighbors of the
vertex $i$-th is denoted as
${\mathcal{N}}_{i}=\\{j\in{\mathcal{V}}:(i,j)\in{\mathcal{E}}\\}$.

Figure 2: Control diagram of a nano-quadrotor

###  II-B Nano-quadrotor model description

Let us take into account a nano-quadrotor system in Fig. 1 including: body
frame $\\{Oxyz\\}_{B}$, position
$p_{i}(t)=\big{[}p_{x,i}(t),~{}p_{y,i}(t),~{}p_{z,i}(t)\big{]}^{T}\in{\mathbb{R}}^{3}$
, velocity
$v_{i}(t)=\big{[}v_{x,i}(t),~{}v_{y,i}(t),~{}v_{z,i}(t)\big{]}^{T}\in{\mathbb{R}}^{3}$,
and roll/pitch/yaw angles $\phi_{i}(t)$/${\theta}_{i}(t)$/$\psi_{i}(t)$ in the
Earth-fixed frame $\\{Oxyz\\}_{E}$. In addition, the rotors $\Omega_{1}(t)$
and ${\Omega}_{3}(t)$ ($\Omega_{2}(t)$ and ${\Omega}_{4}(t)$) rotates
clockwise (anticlockwise) in order to generate a thrust force and moments.
While the thrust force lifts the vehicle along $z$-axis, the moments rotate
the system following $x$-, $y$-, and $z$-axes. Thus, the vehicle is able to
move to any positions in the three-dimensional space by adjusting the thrust
force and the moments. Inspired by [19], this paper successfully applied the
cascade control strategy (see Fig. 2) that includes inner and outer loops in
each controlled quadrotors. Based on the setups, the dynamics of the nano-
quadrotor in the view of the outer loop can be described as the double
integrator:

$\displaystyle\begin{cases}\dot{p}_{i}(t)=v_{i}(t),\\\
\dot{v}_{i}(t)=u_{i}(t),\\\ \end{cases}$ (1)

where
$u_{i}(t)=\big{[}u_{x,i}(t),u_{y,i}(t),u_{z,i}(t)\big{]}^{T}\in{\mathbb{R}}^{3}$
is a control input of the vehicle, representing the accelerations along $x$-,
$y$-, and $z$-axes. It should be remarked that the nano-quadrotor systems
enable us to independently design the outer-loop controller from which the
control input $u_{i}(t)$ is used to calculated the thrust force $T_{i}^{*}(t)$
and the attitude reference $\phi_{i}^{*}(t),~{}{\theta}_{i}^{*}(t)$ (refer Eq.
(2)) toward the inner-loop controller. Since the fast dynamics of the inner
loop, it can be assumed that the attitude immediately tracks its desired value
(see [8, 18, 5] and references therein). Particularly, the relationship among
$u_{i}(t),T_{i}^{*}(t),\phi_{i}^{*}(t)$ and ${\theta}_{i}^{*}(t)$ is given by

$\displaystyle\begin{cases}T_{i}^{*}(t)=m_{i}\sqrt{u_{x,i}^{2}(t)+u_{y,i}^{2}(t)+\big{(}u_{z,i}(t)+g\big{)}^{2}},\\\
\phi_{i}^{*}(t)=\arcsin\frac{m_{i}u_{x,i}^{2}\\!(t)\sin(\psi_{i}\\!(t))-m_{i}u_{y,i}^{2}\\!(t)\cos(\psi_{i}\\!(t))}{T_{i}^{*}(t)},\\!\\!\\\
{\theta}_{i}^{*}(t)=\arctan\frac{u_{x,i}(t)\cos(\psi_{i}(t))+u_{y,i}(t)\sin(\psi_{i}(t))}{u_{z,i}(t)+g},\end{cases}$
(2)

where $m_{i}$ and $g$ denote the mass and the gravitational acceleration of
the $i$-th nano-quadrotor, respectively.

###  II-C Formation description

Figure 3: Formation shape of four nano-quadrotors.

In this letter, the group of nano-quadrotor systems is considered as the
undirected graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathcal{A}})$,
in which each nano-quadrotor is a vertex in ${\mathcal{V}}$. Additionally, a
pair $(i,j)\in{\mathcal{E}}$ implies that the $i$-th nano-quadrotor can invoke
states of the $j$-th nano-quadrotor and vice versa. In an attempt to our
control objectives, the group of $N$ nano-quadrotors is navigated by a
proposed control algorithm to reach an anticipated formation. Generally, a
formation of the multi-quadrotor is considered as a geometric shape in the
three-dimensional space that satisfies some prescribed constraints achieved
and preserved by the group of nano-quadrotors.

In what follows, let us define a reference trajectory of the position-based
formation $r(t):{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}^{3}$ and a formation
shape vector of the $i$-th nano-quadrotor
$f_{p,i}(t):{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}^{3}$, as seen in Fig. 3.
For more details of this formation shape, $f_{p,i}(t)$ is a continuously twice
differentiable function, and
$\|\dot{f}_{p,i}(t)\|=\|{f}_{v,i}(t)\|<\varpi_{fv}$,
$\|\ddot{f}_{p,i}(t)\|<\varpi_{fa}$, where $\varpi_{fv}$ and $\varpi_{fa}$ are
positive constants. Furthermore, the $i$-th nano-quadrotor only knows its
position in the formation via $f_{p,i}(t)$. In the group, at least a nano-
quadrotor knows the reference trajectory of the formation $r(t)$. Let
$\delta_{i}=1$ if the $i$-th nano-quadrotor knows $r(t)$ and $\delta_{i}=0$
otherwise. Let us consider that $r(t)$ satisfies the following

$\displaystyle\dot{r}(t)=v_{0},$ (3)

where $v_{0}\in{\mathbb{R}}^{3}$ is the constant reference velocity of the
formation.

###### Definition 1

The multiple nano-quadrotors (1) are said to achieve the state formation
specified by the vectors $f_{p,i}(t)$, $i\in{\mathcal{V}}=\\{1,2,\dots,N\\}$
for any given bounded initial states if

$\displaystyle\lim\limits_{t\rightarrow\infty}\left(p_{i}(t)-f_{p,i}(t)-r(t)\right)=0,~{}i\in{\mathcal{V}}.$
(4)

######  Lemma 1 ([10])

By letting
$\Delta=\boldsymbol{diag}\big{(}\delta_{1},~{}\delta_{2},\ldots,\delta_{N}\big{)}$
$\in{\mathbb{R}}^{N\times N}$, if the graph ${\mathcal{G}}$ is undirected and
strongly connect, the matrices ${\mathcal{L}}+\Delta$ and
${\mathcal{M}}=({\mathcal{L}}+\Delta)\otimes I_{3}$ are symmetric positive-
definite where ${\mathcal{L}}$ is the Laplace matrix of the graph
${\mathcal{G}}$.

For the convenience, the tracking errors
$e_{p,i}(t),~{}e_{v,i}(t):{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}^{3}$ of the
$i$-th nano-quadrotor can be defined as follow:

$\displaystyle e_{p,i}(t)$ $\displaystyle=p_{i}(t)-f_{p,i}(t)-r(t),$
$\displaystyle e_{v,i}(t)$ $\displaystyle=v_{i}(t)-f_{v,i}(t)-v_{0}.$ (5)

###### Remark 1

When it comes to the complicated trajectories, they are generally divided into
a sequence of desired points that describe the desired position of the
formation. Hence, the reference trajectory of the group of nano-quadrotors can
be established by combining many straight lines connected two consecutive
points in the sequence, i.e., each straight line is considered as a constant
velocity represented in (3).

###### Remark 2

Let us consider an obstacle as another agent freely moving in the experimental
space. Further, the position of this agent is available in the other agents.
This letter mainly focus on the scenario in which the obstacle only appears
after the group of nano-quadrotors completes the given formation shape.

###  II-D Collision-free condition

Figure 4: Two virtual zones created surrounding a nano-quadrotor

With the aim of achieving the collision-free, this letter proposes a finite
cut-off potential function. Each nano-quadrotor system possesses its own
potential function besides a formation control law. While this law manages the
anticipated formation flight of the group of multiple nano-quadrotor systems,
the control input generated from the proposed potential function guarantees
that there are no collisions during the flight time. How to manipulate the two
control inputs can be described by creating three zones, i.e., a risky zone, a
cautionary zone, and a non-jeopardising zone in Fig. 4. Let us define the
radii of the risky and the cautionary zones of the $i$-th nano-quadrotor as
$d_{i}^{(r)}$ and $d_{i}^{(c)}$, respectively. There are three possible
circumstances around the $i$-th nano-quadrotor when an object appears in:

  * •

Non-jeopardising zone: the $i$-th nano-quadrotor solely tracks its desired
spot in the anticipated formation shape.

  * •

Cautionary zone: the $i$-th nano-quadrotor prepares for jeopardy the detected
object possibly causes.

  * •

Risky zone: repulsive forces are generated to steer the $i$-th nano-quadrotor
away from the jeopardising object.

##  III Collision-free formation control

The consideration of collisions to controlling multiple vehicles is one of the
most challenges for automated driving. Autonomous systems can be interrupted
by some collisions among group members and between vehicles and an obstacle.
The aim of this section is to propose an algorithm to manoeuvre nano-quadrotor
systems tracking their given trajectories as well as avoiding collisions.

###  III-A Finite cut-off potential function

This part proposes a novel finite cut-off potential function
$\Phi_{ij}(d_{ij}):{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}_{+}$ (see Fig. 5)
that describes the impact of the $j$-th nano-quadrotor on the $i$-th nano-
quadrotor ($i\neq j$) as follows:

$\displaystyle\Phi_{ij}(d_{ij})=f_{ij}(d_{ij}|\mu_{ij})+\lambda_{ij}g_{ij}(d_{ij}),$
(6)

where $d_{ij}$ denotes the Euclidean distance between the $i$-th and the
$j$-th nano-quadrotors; two scalars $\lambda_{ij}$ and $\mu_{ij}$ are positive
constants such that

$\displaystyle{\lambda_{ij}<\mu_{ij},~{}d_{i}^{(r)}<d_{i}^{(c)}<\infty}.$ (7)

Next, let $f_{ij}(d_{ij}|\mu_{ij}),~{}\forall
d_{ij}\in\big{[}0,~{}\infty\big{)}$ as:

$\displaystyle
f_{ij}(d_{ij}|\mu_{ij})=\begin{cases}\displaystyle\frac{\big{(}d_{i}^{(r)}-d_{ij}\big{)}^{3}}{d_{ij}\\!+\\!d_{i}^{(r)3}\mu_{ij}^{-1}}{,}&\\!\\!\text{if}~{}d_{ij}\\!\in\\!\big{[}0,~{}d_{i}^{(r)}\big{]},\\\
\displaystyle
0{,}&\\!\\!\text{if}~{}d_{ij}\\!\in\\!\left(d_{i}^{(r)},~{}\infty\right).\end{cases}$
(8)

Figure 5: The proposed finite cut-off potential function (6).

Further, a smooth differentiable step function
$g_{ij}(d_{ij}):{\mathbb{R}}_{+}\rightarrow\big{[}0,~{}1\big{]}$ is defined
(9). Based on the definitions of the differentiable step function
$g_{ij}(d_{ij})$ and the function $f_{ij}(d_{ij}|\mu_{ij})$, one has the
derivative of the proposed potential function (10). Furthermore, all the
properties of the above-defined functions are provided in Appendix A. The
smoothness and the finiteness of the proposed potential function (6) depicted
in Fig. 5 afford us a great opportunity to implement our method on
experimentation platforms.

$\displaystyle
g_{ij}(d_{ij})=\begin{cases}0{,}&\text{if}~{}~{}d_{ij}\in\big{(}0,~{}d_{i}^{(r)}\big{]},\\\
\displaystyle\bigg{(}\frac{d_{ij}-d_{i}^{(r)}}{d_{i}^{(c)}-d_{i}^{(r)}}\bigg{)}^{3}\sum_{k=0}^{2}\binom{k+2}{k}\binom{5}{2-k}\bigg{(}\frac{d_{i}^{(r)}-d_{ij}}{d_{i}^{(c)}-d_{i}^{(r)}}\bigg{)}^{k}{,}&\text{if}~{}~{}d_{ij}\in\big{(}d_{i}^{(r)},~{}d_{i}^{(c)}\big{)},\\\
1{,}&\text{if}~{}~{}d_{ij}\in\big{[}d_{i}^{(c)},~{}\infty\big{)}.\end{cases}$
(9)

$\displaystyle\frac{\partial\Phi_{ij}(d_{ij})}{\partial
d_{ij}}=\begin{cases}\displaystyle\frac{-\big{(}d_{i}^{(r)}-d_{ij}\big{)}^{2}\big{(}2d_{ij}+3d^{(r)3}_{i}\mu_{ij}^{-1}+d_{i}^{(r)}\big{)}}{\big{(}d_{ij}+d^{(r)3}_{i}\mu_{ij}^{-1}\big{)}^{2}}{,}&\text{if}~{}d_{ij}\in\big{[}0,~{}d_{i}^{(r)}\big{)},\\\
\displaystyle\lambda_{ij}\frac{\partial g_{ij}(d_{ij})}{\partial
d_{ij}}{,}&\text{if}~{}d_{ij}\in\big{[}d_{i}^{(r)},~{}d_{i}^{(c)}\big{)},\\\
0{,}&\text{if}~{}d_{ij}\in\big{[}d_{i}^{(c)},~{}\infty\big{)}.\end{cases}$
(10)

###  III-B Control synthesis

This section presents a control algorithm for multiple nano-quadrotors with
the purpose of steering the group of nano-quadrotors from initial positions to
form the anticipated formation shape (see Fig. 3). Moreover, the multiple
nano-quadrotors also follow the given formation trajectory (3), and there is
no collision between nano-quadrotors in the group and obstacles as well.
Because of such main goals, the designed control algorithm is constructed from
two parts, i.e., a collision-free input $u_{i}^{c}(t)$ and a formation control
input $u_{i}^{f}(t)$ as follows:

$\displaystyle u_{i}(t)$ $\displaystyle=u_{i}^{c}(t)+u_{i}^{f}(t),$ (11)
$\displaystyle u_{i}^{c}(t)$
$\displaystyle=\sum_{j\in{\mathcal{N}}_{i}}\frac{\partial\Phi_{ij}}{\partial
d_{ij}}\dot{d}_{ij},$ $\displaystyle u_{i}^{f}(t)$
$\displaystyle=\Gamma_{i}\Big{(}-\gamma_{p}\delta_{i}e_{p,i}(t)-\gamma_{v}\delta_{i}e_{v,i}(t)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\gamma_{p}\sum_{j\in{\mathcal{N}}_{i}}\ell_{ij}\big{(}e_{p,i}(t)-e_{p,j}(t)\big{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\gamma_{v}\sum_{j\in{\mathcal{N}}_{i}}\ell_{ij}\big{(}e_{v,i}(t)-e_{v,j}(t)\big{)}\Big{)},$

where
$\Gamma_{i}=\prod_{j\in{\mathcal{N}}_{i}}\big{(}1-g_{ij}(d_{ij})\big{)},~{}\big{(}0\leq\Gamma_{i}\leq
1\big{)}$, $\gamma_{p}$ and $\gamma_{v}$ are the positive scalars.

Let us consider all the circumstances in which there exists an obstacle or
another nano-quadrotor inside the detection zone of the $i$-th nano-quadrotor.
First, if the $i$-th nano-quadrotor detects an object (e.g., the $j$-th
quadrotor) in its risky zone, i.e., $d_{ij}\leq d_{i}^{(r)}$, the part
$u_{i}^{c}(t)$ is non-zero to the $i$-th nano-quadrotor. The control action
(11) prefers avoiding collisions with the detected object to forming the
formation shape. In light of A.9 in Appendix A, the purpose of this part is to
decrease the value of the function $\Phi_{ij}(d_{ij})$, leading to an increase
of the distance from the $i$-th nano-quadrotor to the detected object. Next,
in the second circumstance, the $j$-th nano-quadrotor is detected in the
cautionary zone of the $i$-th nano-quadrotor, i.e., $d_{i}^{(r)}<d_{ij}\leq
d_{i}^{(c)}$. Both two terms of the dedicated control input (11) manage the
system. In which, the part $u_{i}^{c}(t)$ is capable of keeping the detected
object inside the cautionary zone and reducing the repulsive force generated
when this object is in the risky zone. Meanwhile, the part $\Gamma_{i}$
regulates the impact of the consensus formation control (the part
$u_{i}^{f}(t)$). In the last circumstance, there is no object detected in the
range of the cautionary zone. Thanks to the property of the function
$\Phi_{ij}(d_{ij})$ (A.9 in Appendix A), only the part $u_{i}^{f}(t)$ (11)
manoeuvres the $i$-th nano-quadrotor. The aim of this part is to drive the
vehicle such that the group of nano-quadrotors achieves the anticipated
formation shape (see Fig. 3).
Next, let us investigate the last circumstance by constructing the closed-loop
error dynamic model of the group of multiple nano-quadrotor systems. In this
circumstance, there is no object in the cautionary zone of the $i$-th nano-
quadrotor, i.e., $\Gamma_{i}=1$. From the $i$-th nano-quadrotor dynamic model
(1) and the tracking errors (5), one obtains the following closed-loop error
dynamics:

$\displaystyle e_{p}(t)=$ $\displaystyle e_{v}(t),$ $\displaystyle e_{v}(t)=$
$\displaystyle-\gamma_{p}\delta_{i}e_{p,i}(t)-\gamma_{v}\delta_{i}e_{v,i}(t)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\gamma_{p}\sum_{j\in{\mathcal{N}}_{i}}\ell_{ij}\big{(}e_{p,i}(t)-e_{p,j}(t)\big{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\gamma_{v}\sum_{j\in{\mathcal{N}}_{i}}\ell_{ij}\big{(}e_{v,i}(t)-e_{v,j}(t)\big{)}.$
(12)

Then, the closed-loop error dynamic of $N$ nano-quadrotors follows:

$\displaystyle\boldsymbol{\dot{e}}_{p}(t)$
$\displaystyle=\boldsymbol{\dot{e}}_{v}(t),$
$\displaystyle\boldsymbol{\dot{e}}_{v}(t)$
$\displaystyle=\gamma_{p}{\mathcal{M}}\boldsymbol{e}_{p}(t)-\gamma_{v}{\mathcal{M}}\boldsymbol{e}_{v}(t),$
(13)

where
$\boldsymbol{e}_{p}(t)=\big{[}e_{p,1}^{T}(t),~{}e_{p,2}^{T}(t),\ldots,e_{p,N}^{T}(t)\big{]}^{T}$
and $\boldsymbol{e}_{v}(t)=$
$\big{[}e_{v,1}^{T}(t),~{}e_{v,2}^{T}(t),\ldots,e_{v,N}^{T}(t)\big{]}^{T}$.
The following theorem provides a formation control algorithm for the multiple
nano-quadrotor systems, by which there are no collisions among the group
members and obstacles.

###### Theorem 1

(Proof in Appendix A) Let us consider the multiple nano-quadrotor systems (1).
Suppose that the graph ${\mathcal{G}}$ is undirected and strongly connected,
and there exist positive scalar coefficients $\lambda_{ij},\mu_{ij},\gamma$,
$\gamma_{p}$, $\gamma_{v}$, ${\theta}_{p}$, and ${\theta}_{v}$ such that:

$\displaystyle\mu_{ij}>\Phi_{ij}(d_{ij}(0)),~{}\forall
i\in{\mathcal{V}},~{}j\in{\mathcal{N}}_{i},$ $\displaystyle
0<\gamma_{p}-{\theta}_{p},$ $\displaystyle
0<\left(\gamma_{v}-{\theta}_{v}\right){\mathcal{M}}^{2}-\gamma{\mathcal{M}},$
$\displaystyle
0<\left(\gamma_{p}+\gamma_{v}\gamma\right){\mathcal{M}}^{2}-\gamma^{2}{\mathcal{M}},$
(14)

where ${\mathcal{M}}$ is defined in Lemma 1. Then, under the control input
(11), the following statements hold:

  1. i)

The group of nano-quadrotors is collision-free, and

  2. ii)

State formation in Definition 1 is achieved.

##  IV Experimental Results

Figure 6: Experimental scenario of five Crazyflies.

(a)

(b)

(c)

(d)

Figure 7: Actual tests of five Crazyflies: (a) taking off at ground locations,
(b) moving to desired spots of an anticipated formation shape, (c) achieving
the formation shape, (d) tracking a straight trajectory. Figure 8: Distances
between couples of agent-Crazyflies. Figure 9: Distances between an agent-
Crazyflies and an obstacle-Crazyflies. Figure 10: Position snapshots of agent-
Crazyflies.
Figure 11: Position snapshots of agent-Crazyflies moving through an obstacle.

By carrying out flight experiments, this section shows a powerful
demonstration of the proposed algorithm introduced in Section III. The
multiple nano-quadrotors utilized in our flight experiments were Crazyflies
(see Fig. 1) developed by Bitcraze [26]. To be specific, Four Crazyflies were
employed as agents to form a given rectangular shape in Fig. 3 and another
Crayzlife was used as an obstacle. In addition, each Crazyflie was equipped
with an indoor positioning receiver to perceive its position in an
experimental room. This system has afforded us a Python library to send outer-
loop control commands via a wireless communication network.

In order to validate our method, the scenario was constructed by three main
stages in Fig. 7 (total of 28-second flight time). At the first stage, the
four agent-Crazyflies were placed at initial positions (see Fig. 6(a)).
Trajectories of the four agent-Crazyflies from such initial positions to their
desired spots were designed such that they crossed each other (see Fig. 6(b)).
This setup possibly resulted in collision points and offered us a golden
opportunity to verify the proposed algorithm. The aim of the four agent-
Crazyflies at the next stage was to form the anticipated formation shape in
Fig. 3 without any collisions among the group members (see Fig. 6(c)). At the
final stage, the four agent-Crazyflies tracked a straight trajectory in Fig.
6(d). There was another Crazyflie described as an obstacle at a middle point
of the trajectory. For more details, the initial positions of the four agent-
Crazyflies and the obstacle-Crazyflie were placed at different locations in
the experimental room, while their desired formation shapes of all the agent-
Crazyflies were selected as: $f_{p,1}=[0.4,0.45,~{}0.0]^{T}$ ,
$p_{p,2}=[-0.4,~{}0.45,~{}0.0]^{T}$, $p_{p,3}=[0.4,~{}-0.45,~{}0.0]^{T}$, and
$p_{p,4}=[-0.4,~{}-0.45,~{}0.0]^{T}$, $p_{ob}=[0.2,~{}0.2,~{}0.4]^{T}$. Based
on Theorem 1, all the coefficients were chosen as:
$\gamma=0.2,\gamma_{v}=3,\gamma_{p}=2,{\theta}_{p}=1.3,{\theta}_{v}=1.3,\mu_{ij}=0.5,\lambda_{ij}=10^{-3},d_{i}^{(r)}=0.4$m,
and $d_{i}^{(c)}=0.7$m. Laplacian matrix was selected such that each nano-
quadrotor is able to communicate with the other nano-quadrotors. The reference
velocity was set $v_{0}=0.4$ m/s.

As mentioned above, Figs. 10 and 11 illustrated position snapshots of the four
agent-Crazyflies in the 28-second flight time. The four agent-Crazyflies
smoothly moved to their desired spots until the $6^{\text{th}}$ second when
they sensed near objects that might occur collisions. In the following six
seconds, collision avoidance was witnessed in Fig. 10. Distances between
couples of the agent-Crazyflies were shown in Fig. 8. When the agent-
Crazyflies moved near the collision points, these distances gradually
decreased and surpassed $d_{i}^{(r)}=0.4$m at the $6^{\text{th}}$ second.
Thanks to the advances of the proposed potential function, the part $u^{c}(t)$
of the control input (11) was fairly activated to create repulsive forces
among the agent-Crazyflies. The activated repulsive forces smoothly pushed the
agent-Crazyflies far away from the collision points. These forces guaranteed a
safe distance among the agent-Crazyflies. Further, the lines in Fig. 8 from
the $6^{\text{th}}$ to the $14^{\text{th}}$ seconds never hit the zero-line,
clearly illustrating no collisions among the agent-Crazyflies. After the
14-second fight time, all the agent-Crazyflies reached their desired spots in
the formation shape (see Fig. 10). The consensus formation flight was
achieved. At the next stage, the formation trajectory was built as a straight
line in Fig. 6(d). Due to the presence of the obstacle, the four agent-
Crazyflies changed their movements to avoid this obstacle-Crazyflie. This
execution was shown from the $19^{\text{th}}$ to the $28^{\text{th}}$ seconds
in Fig. 11. Fig. 9 illustrated that distances between the four agent-
Crazyflies and the obstacle-Crazyflie never hit the zero-line. This highly
confirmed that no obstacle collision occurs. Please refer to a video at
https://bit.ly/37nno4t for further observations of the actual experiments.

##  V Conclusion

This letter has addressed the formation control of multiple nano-quadrotor
systems regarding collisions among group members and between vehicles and an
obstacle. The collision-free formation protocol has been presented to deal
with such a problem. The effectiveness and validity of the proposed method
have been vividly demonstrated by the actual experiments via Crazyflies nano-
quadrotors. In future work, we plan to study the control of multiple nano-
quadrotor systems under cyber-attacks in which some agents lose control and
are managed by attackers.

## APPENDIX A
Properties of the proposed potential functions

$\bullet$ Properties of the function $f_{ij}(d_{ij}|\mu_{ij})$

  * A.1

$0\leq f_{ij}(d_{ij}|\mu_{ij})\leq\mu_{ij}$

  * A.2

Derivative $\displaystyle\frac{\partial f_{ij}(d_{ij}|\mu_{ij})}{\partial
d_{ij}}$ exists and is continuous.

$\bullet$ Properties of the function $g_{ij}(d_{ij})$

  * A.3

$g_{ij}(d_{ij})$ is continuous and differentiable $\forall
d_{ij}\in\big{(}0,~{}\infty\big{)}$, $0<g_{ij}(d_{ij})<1,~{}\forall
d_{ij}\in\big{(}d_{i}^{(r)},~{}d_{i}^{(c)}\big{)}$,

  * A.4

$\frac{\partial g_{ij}(d_{ij})}{\partial d_{ij}}$ is continuous $\forall
d_{ij}\in\big{(}0,~{}\infty\big{)}$,

  * A.5

$\frac{\partial g_{ij}(d_{ij})}{\partial d_{ij}}>0,~{}\forall
d_{ij}\in\big{(}d_{i}^{(r)},~{}d_{i}^{(c)}\big{)}$,
$\frac{\partial^{k}g_{ij}(d_{ij})}{\partial d_{ij}^{k}}=0$, $\forall
d_{ij}\in\big{(}0,~{}d_{i}^{(r)}\big{]}\cup\big{[}d_{i}^{(c)},~{}\infty\big{)}$.

$\bullet$ Properties of the function $\Phi_{ij}(d_{ij})$

  * A.6

$0\leq\Phi_{ij}(d_{ij})<\mu_{ij},~{}\forall
d_{ij}\in\big{(}0,~{}\infty\big{)}$,

  * A.7

$\max_{d_{ij}\geq 0}\Phi_{ij}(d_{ij})=\Phi_{ij}(0)=\mu_{ij}$,

  * A.8

$\frac{\partial\Phi_{ij}(d_{ij})}{\partial d_{ij}}$ is bounded and continuous
$\forall d_{ij}\in\big{(}0,~{}\infty\big{)}$,

  * A.9

$\frac{\partial\Phi_{ij}(d_{ij})}{\partial d_{ij}}<0,\forall
d_{ij}\in\big{(}0,d_{i}^{(r)}\big{)}$,
$\frac{\partial\Phi_{ij}(d_{ij})}{\partial d_{ij}}>0,~{}\forall
d_{ij}\in\big{(}d_{i}^{(r)},~{}d_{i}^{(c)}\big{)}$,
$\frac{\partial\Phi_{ij}(d_{ij})}{\partial d_{ij}}=0,~{}\forall
d_{ij}\in\big{[}d_{i}^{(c)},~{}\infty\big{)}$.

## APPENDIX B
Proof of Theorem 1

In the scope of this study, collision points are mainly addressed when nano-
quadrotors move from their initial positions to their desired spots in the
formation shape (Fig. 3). Collisions between vehicles and an obstacle are also
considered after the group members reach their anticipated formation. The part
$u_{i}^{c}(t)$ of the proposed control law (11) guarantees that nano-
quadrotors prefer avoiding jeopardizing objects inside their cautionary zones
to forming the given formation shape. This means that a nano-quadrotor tends
to automatically change its direction with a view to placing all the detected
objects outside its cautionary zone, which achieves i).
By utilizing Schur’s complement, the condition (14) and Lemma 1 hold that:

$\displaystyle{\mathbb{P}}=\left[\begin{array}[]{cc}(\gamma_{p}+\gamma_{v}\gamma){\mathcal{M}}^{2}&\gamma{\mathcal{M}}\\\
\gamma{\mathcal{M}}&{\mathcal{M}}\\\ \end{array}\right]>0.$ (17)

Next, let us take a Lyapunov function candidate as follows:

$\displaystyle
V(t)=\big{[}\boldsymbol{e}_{p}^{T}(t),~{}\boldsymbol{e}_{v}^{T}(t)\big{]}{\mathbb{P}}\big{[}\boldsymbol{e}_{p}^{T}(t),~{}\boldsymbol{e}_{v}^{T}(t)\big{]}^{T}$
(18)

The time-derivative of (18) along with the solution of (13) is represented by:

$\displaystyle\dot{V}(t)=$
$\displaystyle(\gamma_{p}+\gamma_{v}\gamma)\boldsymbol{e}_{p}^{T}(t){\mathcal{M}}^{2}\boldsymbol{e}_{v}(t)+\gamma\boldsymbol{e}_{v}^{T}(t){\mathcal{M}}\boldsymbol{e}_{v}(t)$
$\displaystyle+\gamma\boldsymbol{e}_{p}^{T}(t){\mathcal{M}}\boldsymbol{\dot{e}}_{v}(t)+\boldsymbol{e}_{v}^{T}(t){\mathcal{M}}\boldsymbol{\dot{e}}_{v}(t)$
$\displaystyle=$
$\displaystyle\\!-\\!\gamma_{p}\gamma\boldsymbol{e}_{p}^{T}{\mathcal{M}}^{2}\boldsymbol{e}_{p}(t)\\!-\\!\boldsymbol{e}_{v}^{T}(t)(\gamma_{v}{\mathcal{M}}^{2}\\!-\\!\gamma{\mathcal{M}})\boldsymbol{e}_{v}(t)$
$\displaystyle=$
$\displaystyle-\big{[}\boldsymbol{e}^{T}_{p}(t),~{}\boldsymbol{e}_{v}^{T}(t)\big{]}{\mathbb{M}}\big{[}\boldsymbol{e}_{p}^{T}(t),~{}\boldsymbol{e}_{v}^{T}(t)\big{]}^{T},$
(19)

where

$\displaystyle{\mathbb{M}}=\left[\begin{array}[]{cc}\\!\\!\\!\gamma(\gamma_{p}\\!-\\!{\theta}_{p}){\mathcal{M}}^{2}&0\\\
0&(\gamma_{v}\\!-\\!{\theta}_{v}){\mathcal{M}}^{2}\\!-\\!\gamma{\mathcal{M}}\\\
\end{array}\right].$ (22)

By invoking the conditions (14), the matrix ${\mathbb{M}}$ is positive
definite. From the fact that
${\mathbb{M}}\geq\frac{\lambda_{min}({\mathbb{M}})}{\lambda_{max}({\mathbb{P}})}$,
where ${\mathbb{P}}$ is also a positive finite matrix, let us choose:

$\displaystyle\zeta=\frac{2\lambda_{min}({\mathbb{M}})}{\lambda_{max}}.$

Then, (19) gives:

$\displaystyle\dot{V}(t)\leq-\zeta V(t),$ (23)

for $\forall t\geq 0$. Therefore, in light of the comparison lemma [27, pp.
102], one has:

$\displaystyle 0\leq V(t)\leq V(0)e^{-\zeta t}.$

Thanks to the assumptions of a given reference trajectory, initial tracking
errors are bounded, leading to $V(0)<\infty$. For all bounded $V(0)$,
$\displaystyle\lim_{t\rightarrow\infty}V(0)e^{-\zeta t}=0$. This derives that
$\displaystyle\lim_{t\rightarrow\infty}V(t)=0$. Because of the positive
definite matrix ${\mathbb{P}}$ (17),
$\displaystyle\lim_{t\rightarrow\infty}\boldsymbol{e}_{p}(t)=0$ and
$\displaystyle\lim_{t\rightarrow\infty}\boldsymbol{e}_{v}(t)=0$, which proves
ii).

## References

  * [1] B. Chen and H. H. Cheng, “A review of the applications of agent technology in traffic and transportation systems,” _IEEE Trans. Intell. Transp. Syst._ , vol. 11, no. 2, pp. 485–497, 2010.
  * [2] C. Li, X. Yu, T. Huang, and X. He, “Distributed optimal consensus over resource allocation network and its application to dynamical economic dispatch,” _IEEE Trans. Neural Netw. Learn. Syst._ , vol. 29, no. 6, pp. 2407–2418, 2017.
  * [3] P. Zhao, S. Suryanarayanan, and M. G. Simoes, “An energy management system for building structures using a multi-agent decision-making control methodology,” _IEEE Trans. Ind. Appl._ , vol. 49, no. 1, pp. 322–330, 2012\.
  * [4] H. Li, P. Xie, and W. Yan, “Receding horizon formation tracking control of constrained underactuated autonomous underwater vehicles,” _IEEE Trans. Ind. Electron._ , vol. 64, no. 6, pp. 5004–5013, 2016.
  * [5] X. Dong, Y. Zhou, Z. Ren, and Y. Zhong, “Time-varying formation tracking for second-order multi-agent systems subjected to switching topologies with application to quadrotor formation flying,” _IEEE Trans. Ind. Electron._ , vol. 64, no. 6, pp. 5014–5024, 2016.
  * [6] I. Bayezit and B. Fidan, “Distributed cohesive motion control of flight vehicle formations,” _IEEE Trans. Ind. Electron._ , vol. 60, no. 12, pp. 5763–5772, 2012.
  * [7] U. R. Mogili and B. Deepak, “Review on application of drone systems in precision agriculture,” _Procedia Comput. Sci._ , vol. 133, pp. 502–509, 2018.
  * [8] A. T. Nguyen, N. Xuan-Mung, and S.-K. Hong, “Quadcopter adaptive trajectory tracking control: A new approach via backstepping technique,” _Appl. Sci._ , vol. 9, no. 18, p. 3873, 2019.
  * [9] N. P. Nguyen and S. K. Hong, “Active fault-tolerant control of a quadcopter against time-varying actuator faults and saturations using sliding mode backstepping approach,” _Appl. Sci._ , vol. 9, no. 19, p. 4010, 2019.
  * [10] A. T. Nguyen, T. B. Nguyen, and S. K. Hong, “Dynamic event-triggered time-varying formation control of second-order dynamic agents: Application to multiple quadcopterssystems,” _Appl. Sci._ , vol. 10, no. 8, p. 2814, 2020\.
  * [11] D. Yallappa, M. Veerangouda, D. Maski, V. Palled, and M. Bheemanna, “Development and evaluation of drone mounted sprayer for pesticide applications to crops,” in _2017 IEEE Global Humanitarian Technology Conference_. IEEE, 2017, pp. 1–7.
  * [12] J. A. Benito, G. Glez-de Rivera, J. Garrido, and R. Ponticelli, “Design considerations of a small uav platform carrying medium payloads,” in _Design of Circuits and Integrated Systems_ , 2014, pp. 1–6.
  * [13] D. Erdos, A. Erdos, and S. E. Watkins, “An experimental uav system for search and rescue challenge,” _IEEE AERO EL SYS MAG_ , vol. 28, no. 5, pp. 32–37, 2013.
  * [14] W. Hönig, J. A. Preiss, T. S. Kumar, G. S. Sukhatme, and N. Ayanian, “Trajectory planning for quadrotor swarms,” _IEEE Trans. Robot._ , vol. 34, no. 4, pp. 856–869, 2018.
  * [15] C. E. Luis, M. Vukosavljev, and A. P. Schoellig, “Online trajectory generation with distributed model predictive control for multi-robot motion planning,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 2, pp. 604–611, 2020.
  * [16] T. A. Wood, M. Khoo, E. Michael, C. Manzie, and I. Shames, “Collision avoidance based on robust lexicographic task assignment,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 5693–5700, 2020.
  * [17] S. S. Mansouri, C. Kanellakis, B. Lindqvist, F. Pourkamali-Anaraki, A.-A. Agha-Mohammadi, J. Burdick, and G. Nikolakopoulos, “A unified nmpc scheme for mavs navigation with 3d collision avoidance under position uncertainty,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 5740–5747, 2020.
  * [18] F. Liao, R. Teo, J. L. Wang, X. Dong, F. Lin, and K. Peng, “Distributed formation and reconfiguration control of vtol uavs,” _IEEE Trans. Control Syst. Technol._ , vol. 25, no. 1, pp. 270–277, 2016.
  * [19] X. Dong, B. Yu, Z. Shi, and Y. Zhong, “Time-varying formation control for unmanned aerial vehicles: Theories and applications,” _IEEE Trans. Control Syst. Technol._ , vol. 23, no. 1, pp. 340–348, 2014.
  * [20] K. Z. Ang, X. Dong, W. Liu, G. Qin, S. Lai, K. Wang, D. Wei, S. Zhang, S. K. Phang, X. Chen _et al._ , “High-precision multi-uav teaming for the first outdoor night show in singapore,” _Unmanned Systems_ , vol. 6, no. 01, pp. 39–65, 2018.
  * [21] T. B. Nguyen and S. H. Kim, “Distributed flocking bounded control of second-order dynamic multiple polygonal agents,” _IEEE Access_ , vol. 8, pp. 200 170–200 179, 2020.
  * [22] B. Lindqvist, S. S. Mansouri, A.-a. Agha-mohammadi, and G. Nikolakopoulos, “Nonlinear mpc for collision avoidance and control of uavs with dynamic obstacles,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 6001–6008, 2020\.
  * [23] X. Kan, H. Teng, and K. Karydis, “Online exploration and coverage planning in unknown obstacle-cluttered environments,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 5969–5976, 2020.
  * [24] M. Boldrer, M. Andreetto, S. Divan, L. Palopoli, and D. Fontanelli, “Socially-aware reactive obstacle avoidance strategy based on limit cycle,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 2, pp. 3251–3258, 2020.
  * [25] N. T. Binh, P. D. Dai, N. H. Quang, N. T. Ty, and N. M. Hung, “Flocking control for two-dimensional multiple agents with limited communication ranges,” _International Journal of Control_ , pp. 1–8, 2020.
  * [26] Bitcraze, https://bit.ly/2VkMMRE, 2020, [Online; accessed 11-November-2020].
  * [27] H. K. Khalil, “Nonlinear systems,” _Upper Saddle River_ , 2002.