\section*{Supplementary Material}\label{sec:supplementary-material}

\section{Introduction}\label{sec:introduction}

\IEEEPARstart{D}{rone swarms} have a large socio-economic potential and can serve in a variety of real-world applications \cite{floreano_science_2015}.
For example, drones can be leveraged to automatically monitor crops, safely inspect confined spaces, or quickly deliver medicine to inaccessible locations.
Operating these vehicles in swarms could bring increased robustness to failures, larger area coverage, and faster task completion times \cite{coppola_survey_2020}.

Despite this potential, decentralized control has been a limiting factor in the deployment of drone swarms.
For instance, large groups of quadcopters can be used to perform awe-inspiring aerial choreographies in the night sky.
These robotic light shows are a true feat of engineering but individual drones are far from autonomous: their motion is centrally controlled by a ground computer that precomputes their trajectories and continuously monitors their positions.
Hence, the ground computer represents a single point of failure.
Researchers attempted to remove the central computer by equipping drones with hardware that allows them to wirelessly communicate with each other.
Notable examples of these decentralized swarms feature bearing-only formation control \cite{schiano_bearing_2017}, exploration of unknown environments \cite{mcguire_minimal_2019}, as well as flocking with ten fixed-wing drones \cite{hauert_reynolds_2011} or thirty quadcopters \cite{vasarhelyi_optimized_2018}.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{overview.pdf}
    \caption{%
        Photo taken during outdoor experiments with detection annotations.
        The quadcopters avoid collisions with each other and remain cohesive as a group while performing a variety of navigation tasks.
        Each agent detects its neighbors in real-time from omnidirectional images.
        There is no communication of state information between agents, and relative positions and velocities are estimated onboard using local visual inputs.
    }\label{fig:overview}
\end{figure}

\newcommand{0.32 \textwidth}{0.32 \textwidth}
\begin{figure*}[t]
    \begin{subfigure}{\textwidth}
        \includegraphics[width=\textwidth]{method/method/method.pdf}
    \end{subfigure}
    \medskip
    \\
    \begin{subfigure}{0.32 \textwidth}
        \caption{Real-time drone detection}\label{fig:method_detection}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{Multi-agent state tracking}\label{fig:method_tracking}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{Potential-field-based control}\label{fig:method_control}
    \end{subfigure}
    \caption{Overview of the processing steps of the vision-based flocking algorithm for a single time step: a) \textit{detection}, b) \textit{tracking}, and c) \textit{control}.
    Step 1 (\textit{detection}): we detect neighboring agents from omnidirectional images to estimate their relative range and bearing from the camera intrinsics and the drone's known physical size.
    Step 2 (\textit{tracking}): we track the positions and velocities of the detected agents using a linearized observation model of range and bearing, as well as a state transition model that takes the focal agent's egomotion into account.
    Step 3 (\textit{control}): we control the focal agent using a Reynolds-rules-based flocking algorithm that keeps the swarm collision-free and cohesive while following a navigation goal.
    }\label{fig:method}
\end{figure*}

Scaling up decentralized drone swarms is complicated by the limitations of wireless communication.
As the number of robots increases, the communication channels may become saturated and possibly jammed since the data transfer volume scales quadratically with the robot count \cite{cieslewski_data-efficient_2018}.
Frequent retransmissions of messages can lead to delays that render the control of each drone extremely difficult.
Researchers have experimented with different sensory modalities such as sound \cite{basiri_onboard_2016}, but vision seems to be the most scalable approach to address the relative localization problem.
Indeed, there is evidence to support that vision is the primary sensory modality that enables collective motion in animal groups \cite{strandburg-peshkin_visual_2013}.

However, visual detection of outdoor flying drones is challenging because they are relatively small and can fly in environments with large amounts of background clutter and difficult lighting conditions \cite{rozantsev_detecting_2017}.
Additionally, the visual detection algorithm must run onboard the drones, whose own motion may generate motion blur and whose small dimensions and lightweight mass can restrict computational capabilities.
To simplify the relative localization problem, researchers have mounted different types of easily detectable visual markers on the drones \cite{roelofsen_reciprocal_2015,walter_uvdar_2019}.
However, this approach is less general since specialized hardware has to be mounted on each drone, which can result in increased weight and drag, thus reducing the energetic autonomy of the drones.

In recent years, markerless detection of drones has become an active research topic.
In \cite{sapkota_vision-based_2016}, the authors use a boosted cascade of classifiers in combination with visual tracking and finite set filtering to estimate the positions and velocities of markerless drones from a mostly static observer.
However, the method is applied in post-processing and not validated in a sense-and-avoid setting.
In \cite{vrba_marker-less_2020}, the authors propose a combination of stereo vision and convolutional detection to enable leader-follower flight.
However, the to-be-localized agent is always visible in front of the clear sky which simplifies the detection problem and would be impossible to guarantee in a self-organized flock.
Moreover, the limited field of view and processing delays in the proposed system are known to be problematic when flying in dense swarms \cite{vasarhelyi_optimized_2018,soria_influence_2019}.
Other notable examples of markerless detection include approaches based on template matching and morphological filtering \cite{opromolla_vision-based_2018}, as well as convolutional neural networks \cite{wyder_autonomous_2019}.

Our previous work \cite{schilling_learning_2019} proposes a fundamentally different approach to visual flocking based on imitation learning.
Rather than detecting neighboring agents, we predict flocking algorithm commands directly from omnidirectional visual inputs which allow the drones to remain collision-free and cohesive.
The approach is validated with leader-follower experiments in an indoor environment but its reliance on end-to-end learning means that the entire monolithic neural network has to be retrained each time the task and/or visual appearance of the drones change.
Adopting a more modular approach can be beneficial since the flocking algorithm may easily be exchanged for another task-dependent controller, and only the detector would have to be retrained for different drone appearances or environmental conditions.
In all of the above cases, the algorithms are validated on a single agent and not in a multi-robot control setting.

Here, we propose a modular detection and tracking algorithm that enables collision-free and cohesive navigation for drone swarms.
We automatically label images of drones using background subtraction to generate a dataset for the drone detector.
We show that the detector can localize other drones in the presence of background clutter from onboard a flying drone despite being trained on images from a static camera.
We assess the method with a dense group of three real quadcopters that flock in planar configurations in an outdoor environment with difficult lighting conditions.
The omnidirectional camera configuration of the drone is specifically designed to enable safe operation in swarms regardless of the agent configuration.
The overall proposed flocking algorithm is modular since each component, i.e. detection, tracking, and control, is self-contained and can be evaluated independently.
To the best of our knowledge, this is the first entirely vision-based flock that does not depend on visual markers to simplify mutual detections.

\section{Method}\label{sec:method}

The proposed approach to vision-based flocking can be divided into the following steps: \textit{detection}, \textit{tracking}, and \textit{control} (Fig.~\ref{fig:method}).
Firstly, the \textit{detection} module (Fig.~\ref{fig:method_detection}) takes grayscale images from an omnidirectional camera setup as inputs and outputs bounding boxes of nearby drones in real-time (Sec.~\ref{sec:control}).
Secondly, the \textit{tracking} module (Fig.~\ref{fig:method_tracking}) transforms the bounding boxes into range and bearing measurements using the known dimensions of the drones.
Their relative positions and velocities are then estimated from the noisy measurements using a multi-agent state tracker (Sec.~\ref{sec:localization}).
Finally, the \textit{control} module (Fig.~\ref{fig:method_control}) applies a flocking algorithm to the relative positions to obtain high-level control commands that keep the drones collision-free and cohesive (Sec.~\ref{sec:control}).
In the tracking and control steps, we assume that the agents are moving on a horizontal plane.
We mainly introduce this constraint to be able to attribute their mutual repulsion to the flocking algorithm and avoid the effects of physical downwash that may be caused by nearby agents.

\subsection{Real-time monocular drone detection}\label{sec:detection}

\subsubsection{Automatic drone labeling with background subtraction}\label{sec:background_subtraction}

We use background subtraction to automatically generate a labeled image dataset for the object detector (Fig.~\ref{fig:background_subtraction}).
We record images from a stationary camera mounted on a tripod and manually fly a quadcopter within its field of view.
The image data is recorded under varying lighting conditions in both indoor and outdoor environments that contain large amounts of background clutter (Fig.~\ref{fig:dataset}).
For each scene, the camera location and orientation are carefully chosen such that the quadcopter is the dominant source of motion.
The final dataset consists of $9\,891$ training and $1\,931$ testing examples.

In post-processing, we apply a nearest neighbor background subtraction algorithm \cite{zivkovic_efficient_2006} to the sequence of images to learn a background model of the scene.
We extract a foreground mask of the moving parts of the image (and therefore the quadcopter) by computing the element-wise difference of the input image and the background, followed by a thresholding step.
The ground truth label is obtained by filtering out the largest contour present in the foreground mask and enclosing it with an axis-aligned rectangular bounding box.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{method/background_subtraction/background_subtraction.pdf}
    \caption{%
        Schematic overview of the automatic dataset generation process.
        The foreground mask is generated using background subtraction.
        We find that the most precise bounding box labels are obtained by dilating the foreground mask since it eliminates discontinuities that occur due to the mechanical design of the drone.
        Enclosing the dilated mask with a bounding box (red rectangle) overestimates the size of the drone.
        We therefore scale the bounding box down (green rectangle) to obtain a precise label.
        We record six flights of roughly one-minute duration in the above target environment.
    }\label{fig:background_subtraction}
\end{figure}

\newcommand{0.49 \columnwidth}{0.49 \columnwidth}

\begin{figure}[t]
    \begin{subfigure}{0.49 \columnwidth}
        \includegraphics[width=\textwidth]{method/dataset/dronedome.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.49 \columnwidth}
        \includegraphics[width=\textwidth]{method/dataset/workshop.png}
    \end{subfigure

    \medskip

    \begin{subfigure}{0.49 \columnwidth}
        \includegraphics[width=\textwidth]{method/dataset/parking.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.49 \columnwidth}
        \includegraphics[width=\textwidth]{method/dataset/passage.png}
    \end{subfigure}
    \caption{%
        Example images with bounding box annotations from the dataset generated using background subtraction.
        We record image data from a static camera in both indoor (top row) and outdoor (bottom row) environments.
        The environments are selected to maximize the variety of background clutter and lighting conditions.
        We record three flights of roughly one-minute duration in each of the above environments.
    }\label{fig:dataset}
\end{figure}

\subsubsection{Training the real-time drone detector}\label{sec:training_detector}

We train a single-stage convolutional object detector on the automatically annotated image dataset to obtain drone detections.
We opt for the YOLOv3-tiny architecture \cite{redmon_yolov3_2018} due to its favorable tradeoff between detection accuracy and inference speed on embedded devices.
The network architecture is comprised of a total of $13$ convolutional layers that are interspersed with max-pooling operations and leaky rectified linear units (ReLUs).
The final detections are computed independently at two different scales and subsequently filtered by non-maximum suppression to account for bounding box overlaps.

We make a few notable modifications to the training procedure of the original publication \cite{redmon_yolov3_2018}.
Firstly, we replace the mean-squared error localization loss with an objective that is based on the generalized intersection over union (GIoU) \cite{rezatofighi_generalized_2019} since it provides scale invariance.
Secondly, we employ two recently popularized data augmentation techniques to improve detection accuracy: 1) multi-scale training, and 2) mosaic augmentation \cite{bochkovskiy_yolov4_2020}.
In multi-scale training, each training batch is scaled at random by up to $\pm 50\%$ of its side length to make the detector invariant to object scales.
Mosaic training refers to concatenating four random training samples along their spatial dimensions to obtain an image collage.
The resulting four-image mosaic is subsequently cropped randomly in the center to obtain a new training sample.
We use holdout cross-validation to find suitable values for the most critical hyperparameters.

The parameters of the network are initialized with weights that are pre-trained on the COCO \cite{lin_microsoft_2014} dataset.
The detector is then fine-tuned on a single \tit{drone} class with stochastic gradient descent and Nesterov momentum of $\mu = 0.937$.
We modulate the learning rate using a cosine annealing schedule with an initial learning rate $\eta = 0.01$ and a final learning rate $\eta_f = 0.0005$.
We also employ a weight decay term of $\lambda = 0.0005$.
The total loss is computed as

\begin{equation}
    \mcl{L} = k^\txt{bal} \mcl{L}_{\text{loc}} + \mcl{L}_{\text{obj}}
\end{equation}
where $\mcl{L}_{\text{loc}}$ denotes the GIoU-based localization loss and $\mcl{L}_{\text{obj}}$ refers to the binary cross-entropy objectness loss.
We set the hyperparameter $k^\txt{bal} = 0.055$ to balance the losses and to account for their different magnitudes during training.
Note that unlike the original paper, we omit the classification loss term since we are training on a single class.

Our highest-scoring model achieves an average precision (AP@0.5) of $98.9\%$ \cite{everingham_pascal_2010} at a confidence threshold of $p^\txt{conf} = 0.001\%$ on the hold-out test set after $77$ epochs of training.
The model can be trained in less than $\SI{1}{\hour}$ using a batch size of $64$ on a recent GPU such as the Nvidia GeForce RTX 2080Ti.
We also evaluate the model at different image scales on the onboard computer (Sec.~\ref{sec:hardware}) to determine a reasonable speed/accuracy tradeoff of the detector.
We find that performing inference at a resolution of $512 \times 384 $ pixels (in batches of four, one image per camera) provides accurate detections at a frequency of around $\SI{5}{\hertz}$.
For the experiments, we set the confidence threshold to $p^\text{conf} = 50\%$ and use a non-maximum suppression threshold of $p^\text{nms} = 60\%$.

\subsection{Multi-agent localization and tracking}\label{sec:localization}

\subsubsection{Relative localization based on known physical size}\label{sec:relative_localization}

We compute the relative location of the drone detections using their apparent size in the field of view and the camera parameters.
The camera parameters are obtained using the Kalibr visual-inertial calibration toolbox \cite{furgale_unified_2013}.
We use the equidistant camera model for its compatibility with fisheye lenses and its resulting sub-pixel reprojection errors.

To obtain the relative position estimate from the detection, we first compute the unit-norm bearing vector to the center $\bm{\beta}^\txt{ctr}$ and one of the extreme points $\bm{\beta}^\txt{ext}$ of the image bounding box from the intrinsic camera parameters.
We assume the drone can be enclosed by a bounding cube with side length $l$ which is reasonable given its mechanical design (Fig.~\ref{fig:drone_hardware}).
We can then compute the approximate distance to the three-dimensional center of the detected object as

\begin{equation}
    d = \frac{l / 2}{\tan({\alpha})} + l / 2
\end{equation}
where $\alpha = \cos^{-1}(\bm{\beta}^\txt{ctr} \cdot \bm{\beta}^\txt{ext})$ denotes the angle between the unit-norm bearing vectors.
Note that the second term in the above equation accounts for the depth of the object.

\subsubsection{Multi-agent state estimation using random finite sets}\label{sec:multi_agent_tracking}

We use the Gaussian mixture probability hypothesis density (GM-PHD) filter \cite{vo_gaussian_2006} to filter out spurious false-positive detections and to estimate the positions and velocities of nearby agents over time.
We briefly describe the workings of the filter but refer the reader to \cite{vo_gaussian_2006} for more details.
We omit the subscript $i$ to denote the dependence on the focal agent for notational brevity.
The following steps are computed independently for each agent in a decentralized fashion.

\paragraph{Theory}

The GM-PHD filter takes a set of relative localization measurements $\mathcal{Z}_k = \{ \mbf{z}_{k,1}, \dots, \mbf{z}_{k,M_k} \}$ as input and computes an output set of agent states $\mathcal{X}_k = \{ \mbf{x}_{k,1}, \dots, \mbf{x}_{k,J_k} \}$ for each discrete time step $k$.
The states can be described as a single intensity that consists of a weighted sum of Gaussian components of the form
\begin{equation}
    v_k(\mbf{x}_k) = \sum_{i=1}^{J_k} w_k^{(i)} \mcl{N}(\mbf{x}_k; \mbf{m}_k^{(i)}, \mbf{P}_k^{(i)})
\end{equation}
where $w_k^{(i)}$ denotes the weight associated with each of the $J_k$ Gaussian components which are described by their mean $\mbf{m}_k^{(i)}$ and covariance $\mbf{P}_k^{(i)}$.
Each of the Gaussian components is then propagated with a \tit{prediction} and \tit{update} step, similar to the Kalman filter.

The \tit{prediction} step can be formalized as
\begin{equation}
    v_{k|k-1}(\mbf{x}_k) = \sum_{i=1}^{J_k} w_{k|k-1} \mcl{N}(\mbf{x}_k; \mbf{m}_{k|k-1}^{(i)}, \mbf{P}_{k|k-1}^{(i)}) + \gamma(\mbf{z}_k)
\end{equation}
with respective weight, mean, and covariance
\begin{align}
    w_{k|k-1}^{(i)} &= p_{s,k} w_{k-1}^{(i)} \\
    \mbf{m}_{k|k-1}^{(i)} &= \mbf{F}_{k-1} \mbf{m}_{k-1}^{(i)} + \mbf{B}_{k-1} \mbf{u}_{k-1} \\
    \mbf{P}_{k|k-1}^{(i)} &= \mbf{F}_{k-1} \mbf{P}_{k-1}^{(i)} \mbf{F}_{k-1}^\top + \mbf{Q}_{k-1}
\end{align}
where $\mbf{F}_{k-1}$ is the state transition matrix and $\mbf{Q}_{k-1}$ the process noise covariance.
To account for the egomotion of the observing drone, we include a control input matrix $\mbf{B}_{k-1}$ and its control input $\mbf{u}_{k-1}$.
We let $p_{s,k}$ denote the probability that a Gaussian component survives the prediction step.
We assume an adaptive agent birth model $\gamma_k(\mbf{z}_k)$ in which each observation generates a new Gaussian component with weight $w_{\gamma}$, mean $\mbf{m}_{\gamma}$, and covariance $\mbf{P}_{\gamma}$.
We further assume that new agents cannot be spawned from existing ones and that there are no spontaneous births without associated measurement.

The \tit{update} step can be formalized as
\begin{align}
    v_{k}(\mbf{x}_k) &= (1 - p_{d,k}) v_{k|k-1}(\mbf{x}_k) \\
    &+ \sum_{\mbf{z}_k \in Z_k} \sum_{i=1}^{J_k} w_k^{(i)} (\mbf{z}_k) \mcl{N}(\mbf{x}_k; \mbf{m}_{k|k}^{(i)} (\mbf{z}_k), \mbf{P}_{k|k}^{(i)})
\end{align}
with respective weight, mean, and covariance
\begin{align}
    w_k^{(i)}(\mbf{z}_k) &=
        \frac   {p_{d,k} w_{k|k-1}^{(i)} q_k^{(i)} (\mbf{z}_k)}
                {\kappa_k(\mbf{z}_k) + \sum_{j=1}^{J_{k|k-1}} p_{d,k}  w_{k|k-1}^{(j)} q_k^{(j)}(\mbf{z}_k)} \\
    \mbf{m}_{k|k}^{(i)} (\mbf{z}_k) &= \mbf{m}_{k|k-1}^{(i)} + \mbf{K}_k^{(i)} (\mbf{z}_k - \mbf{H}_k \mbf{m}_{k|k-1}^{(i)}) \\
    \mbf{P}_{k|k}^{(i)} &= (\mbf{I} - \mbf{K}_k^{(i)} \mbf{H}_k) \mbf{P}_{k|k-1}^{(i)}
\end{align}
and
\begin{align}
    q_k^{(i)} (\mbf{z}_k) &= \mcl{N}(\mbf{z}_k; \mbf{H}_k \mbf{m}_{k|k-1}^{(i)},\mbf{H}_k \mbf{P}_{k|k-1}^{(i)} \mbf{H}_k^\top + \mbf{R}_k) \\
    \mbf{K}_k^{(i)} &= \mbf{P}_{k|k-1}^{(i)} \mbf{H}_k^\top (\mbf{H}_k \mbf{P}_{k|k-1}^{(i)} \mbf{H}_k^\top + \mbf{R}_k)^{-1}
\end{align}
where $\mbf{H}_k$ is the measurement matrix and $\mbf{R}_k$ the measurement noise covariance.
We let $p_{d,k}$ denote the probability that an agent is detected during the update step.
We model false positive detections as clutter $\kappa_k$.

Since the number of Gaussian components increases at each filter iteration, the intensity quickly becomes computationally intractable.
Therefore, we prune the components according to the following three conditions to guarantee fast tracking performance.
Firstly, we discard components with a weight of less than the truncation threshold of $T$.
Secondly, we merge components with Mahalanobis distance less than the merging threshold of $U$.
Finally, we retain only the $J_\txt{max}$ components with the largest weights.

\paragraph{Implementation}

We model the state of each agent as $\mbf{x}_k = \left[ p_{x,k}, p_{y,k}, v_{x,k}, v_{y,k} \right]^\top $ which consists of relative position $(p_{x,k}, p_{y,k})$ and velocity $(v_{x,k}, v_{y,k})$.
The position and velocity components of the state are two-dimensional since the agents are assumed to fly in a planar configuration at approximately the same altitude.
The control input $\mbf{u}_k$ is defined as the linear velocity of the drone in the body frame which we obtain from the internal state estimate of the autopilot.
Finally, the measurements $\mbf{z}_k = [ d_k, \beta_k ]^\top$ consist of range and bearing, respectively.

The process and measurement noise covariances are modeled as

\begin{align} \label{eq:process_obs_noise}
    \mbf{Q}_k &= \sigma_v^2
        \begin{bmatrix}
            \frac{\Delta_k^4}{4} \mbf{I}_2 & \frac{\Delta_k^3}{2} \mbf{I}_2 \\
            \frac{\Delta_k^3}{2} \mbf{I}_2 & \Delta_k^2 \mbf{I}_2
        \end{bmatrix} &
        &\txt{and} &
    \mbf{R}_k &=
        \begin{bmatrix}
            \sigma_d^2 & 0 \\
            0 & \sigma_\beta^2
        \end{bmatrix}
\end{align}
where $\Delta_k$ denotes the time elapsed since the last measurement and is computed as the difference between consecutive timestamps $\Delta_k = t_k - t_{k-1}$.
We further let $\sigma_v$, $\sigma_d$ and $\sigma_\beta$ denote the standard deviation of the process, range, and bearing noise, respectively.

The state transition function follows a linear Gaussian model and is defined as

\begin{equation}
    f(\mbf{x}_{k-1}, \mbf{u}_k) = \mbf{F}_k \mbf{x}_{k-1} + \mbf{B}_k \mbf{u}_k
\end{equation}
where
\begin{align}
    \mbf{F}_k &=
        \begin{bmatrix}
            \mbf{I}_2 & \Delta_k \mbf{I}_2 \\
            \mbf{0}_2 & \mbf{I}_2
        \end{bmatrix} &
        &\txt{and} &
    \mbf{B}_k &= \Delta_k \mbf{I}_2.
\end{align}
The observation model is nonlinear and consists of measurements of range and bearing

\begin{equation}
    h(\mbf{x}_k) =
        \begin{bmatrix}
            d_k \\
            \beta_k
        \end{bmatrix} =
        \begin{bmatrix}
            \sqrt{p_{x,k}^2 + p_{y,k}^2} \\
            \atantwo(p_{y,k}, p_{x,k})
        \end{bmatrix}
\end{equation}
where we use the two-argument function $\atantwo$ to avoid ambiguities in the conversion from cartesian to polar coordinates.
We linearize the measurement model by computing the Jacobian with respect to the state as

\begin{equation}
    \mbf{H}_k = \frac{\partial h(\mbf{x}_k)}{\partial \mbf{x}_k} =
        \begin{bmatrix}
            \begin{matrix}
                \frac{p_{x,k}}{\sqrt{p_{x,k}^2 + p_{y,k}^2}} & \frac{p_{y,k}}{\sqrt{p_{x,k}^2 + p_{y,k}^2}} \\
                -\frac{p_{x,k}}{p_{x,k}^2 + p_{y,k}^2} & \frac{p_{y,k}}{p_{x,k}^2 + p_{y,k}^2}
            \end{matrix}
            & \mbf{0}_2
        \end{bmatrix}.
\end{equation}

We set the probability of detection to $p_d = 90\%$ to provide a slightly more conservative estimate of the detection performance during real-time inference than the results on our test dataset suggest (Sec.~\ref{sec:training_detector}).
We assume a probability of survival of $p_s = 100\%$ since agents that are detected once should not disappear.
New Gaussian components are initialized directly from the measurements using an adaptive birth model with weight, mean, and covariance
\begin{align}
    w_\gamma &= 10^{-5} \\
    \mbf{m}_\gamma &= \left[ d_k \cos(\beta_k), d_k \sin(\beta_k), 0, 0 \right]^\top \\
    \mbf{P}_\gamma &= \diag\left( \left[ \sigma_p^2, \sigma_p^2, \sigma_v^2, \sigma_v^2 \right] \right)
\end{align}
where $\sigma_p$ and $\sigma_v$ are the standard deviation of the position and velocity which we set to $\SI{1}{\meter}$ and $\SI{10}{\meter\per\second}$, respectively.
The mean is computed by converting the raw measurement from polar to cartesian coordinates, assuming zero initial velocity.
We model false positive detections by assuming that we observe one clutter return per time step and therefore set $\kappa_k = 1 / a^2$ where $a = \SI{10}{\meter}$ is the side length of the virtual arena.
Although the bearing measurements are precise to a single degree, we set its standard deviation to a slightly larger value of $\sigma_\beta = 3^\circ$ to account for factors such as calibration inaccuracies or imprecisions in the detections due to background clutter.
Since we experimentally determined that the range noise varies with the distance to the observer, we model its standard deviation as a function of the measured distance itself (Sec.~\ref{sec:relative_localization_results}).
Finally, we model the truncation threshold as $T = 10^{-5}$, the merging threshold as $U = 0.5$, and the maximum number of components to $J_\txt{max} = 100$.

\subsection{Flocking algorithm}\label{sec:control}

The method described so far could be leveraged by any flocking algorithm.
In this work, we use a control algorithm based on the Reynolds flocking rules \cite{reynolds_flocks_1987} to compute high-level velocity commands from the relative position estimates of nearby drones \cite{schilling_learning_2019}.
The weighted velocity commands are: 1) a repulsive \tit{separation} term to steer nearby drones away from each other, 2) a \tit{cohesion} term to keep the drones close to each other, and 3) a \tit{migration} term that provides a navigation goal to the swarm.
The respective velocity terms for separation, cohesion, and migration can be formalized as
\begin{align}
    \mbf{\vel}_i^\txt{sep} &= k^\txt{sep} \frac{1}{|\set{A}_i|} \sum_{j \in \set{A}_i} \frac{-\mbf{\rel}_{ij}}{\norm{\mbf{\rel}_{ij}}} \\
    \mbf{\vel}_i^\txt{coh} &= k^\txt{coh} \frac{1}{|\set{A}_i|} \sum_{j \in \set{A}_i} {\mbf{\rel}_{ij}} \\
    \mbf{\vel}_i^\txt{mig} &= k^\txt{mig} \frac{\mbf{\rel}^\txt{mig}}{\norm{\mbf{\rel}^\txt{mig}}}
\end{align}
where $\mbf{\rel}_{ij}$ denotes the relative position of the $j$-th agent with respect to agent $i$, and $\set{A}_i$ the set of neighbors of the $i$-th agent.
The migration point is denoted by $\mbf{\rel}^\txt{mig}$ and expressed relative to the body frame.
We further let $k^\txt{sep}$, $k^\txt{coh}$, and $k^\txt{mig}$ denote the gains that modulate the strength of the separation, cohesion, and migration terms, respectively.
The final velocity command is then given by $\mbf{\vel}_i = \tilde{\mbf{\vel}}_i / \norm{\tilde{\mbf{\vel}}_i} \min(\norm{\tilde{\mbf{\vel}}_i}, v^\txt{max})$ where $\tilde{\mbf{\vel}}_i = \mbf{\vel}_i^\txt{sep} + \mbf{\vel}_i^\txt{coh} + \mbf{\vel}_i^\txt{mig}$ is the sum of all velocity terms and $v^\txt{max}$ is the desired maximum speed cutoff.

During the experiments, we set the maximum speed to $v^\txt{max} = \SI{0.5}{\meter \per \second}$, and the separation, cohesion, and migration gains to $k^\txt{sep} = 7$, $k^\txt{coh} = 1$, and $k^\txt{mig} = 1$, respectively.
The gains are chosen such that the agents converge to an equilibrium distance of approximately $\SI{2}{\meter}$ during migration.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{drone_hardware/drone_hardware.pdf}
    \caption{The LeQuad drone features omnidirectional vision via four cameras, a high performance onboard computer with embedded GPU for real-time inference, and a RTK-enabled GNSS to obtain centimeter-accurate ground-truth positions.}\label{fig:drone_hardware}
\end{figure}

\section{Experimental setup}\label{sec:experimental_setup}

\subsection{Hardware}\label{sec:hardware}

We use a custom-built quadcopter named LeQuad for all experiments (Fig.~\ref{fig:drone_hardware}).
Each quadcopter features four FLIR Firefly S global-shutter cameras mounted at a right angle from each other to obtain omnidirectional visual inputs.
Each camera is equipped with an OpenMV ultra-wide angle lens which provides a horizontal and vertical field of view of $ \SI{166}{\degree} $ and $\SI{116}{\degree}$, respectively.
We operate the cameras over a powered USB 2 hub at a binned resolution of $720 \times 540$ to obtain grayscale images at a frequency of $\SI{10}{\hertz}$.
We refrain from using USB 3 to avoid electromagnetic interference with the Drotek F9P RTK-GNSS receiver, which provides centimeter-accurate absolute positions.
We use the Nvidia Jetson TX2 mounted on a ConnectTech Orbitty carrier board as an onboard computer and the Holybro Pixhawk 4 as an autopilot.

\subsection{Software}\label{sec:software}

The onboard computer runs Ubuntu 18.04 bundled with the Linux4Tegra (L4T) distribution, and we use ROS Melodic as a robotics middleware.
The autopilot runs PX4 and is responsible for hardware-triggering the cameras to provide the resulting images with IMU-synchronized timestamps.
The neural networks are trained and evaluated using PyTorch.

\section{Results}\label{sec:results}

We first report the accuracy of the visual relative localization of drones in a controlled indoor environment (Sec.~\ref{sec:relative_localization_results}) and then show vision-based flocking with three real quadcopters performing several navigation tasks outdoors (Sec.~\ref{sec:outdoor_results}).

The proposed approach has also been extensively tested using the Gazebo simulator with up to ten vision-based agents.
For instance, we simulated noisy detections with different false positive and negative rates, as well as processing delays, to find suitable parameter values and operational requirements for the detection, tracking, and control modules.
However, we find that it is difficult to realistically model misdetections in simulation since their distribution highly depends on environmental conditions such as visual clutter.

\subsection{Visual relative localization}\label{sec:relative_localization_results}

We show results on the theoretical performance of the visual relative localization system to test its operational bounds and to find suitable values for the range and bearing noise parameters $\sigma_{d}$ and $\sigma_{\beta}$ (Eq.~\ref{eq:process_obs_noise}).
To this end, we employ a setup similar to the one used for automatic labeling (Sec.~\ref{sec:background_subtraction}) except that we additionally obtain ground truth poses of the observing camera and drone from a motion capture system.
After transforming the visual detections into the frame of the motion capture system, we can directly compare the drone's true position with its estimate obtained using vision in metric space.
We find that the relative localization error varies considerably as a function of the distance to the drone, whereas the error caused by bearing variations is negligible (Fig.~\ref{sec:relative_localization}).

\newcommand{0.48 \columnwidth}{0.48 \columnwidth}
\begin{figure}[t]
    \begin{subfigure}{0.48 \columnwidth}
        \includegraphics[width=\textwidth]{results/localization/range_boxplot.pdf}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.48 \columnwidth}
        \includegraphics[width=\textwidth]{results/localization/bearing_boxplot.pdf}
    \end{subfigure}
    \caption{%
        Comparison of vision-based relative localization errors by \tit{range} (left) and \tit{bearing} (right).
        The bearing estimates are near-constant over the field of view, whereas the range errors increase with distance from the observer.
        Millimeter-accurate ground-truth positions are obtained at $\SI{100}{\hertz}$ using a motion capture system.
        The counts above the boxes indicate the number of measurements used to calculate their statitics.
    }\label{fig:localization}
\end{figure}

\subsection{Collective outdoor navigation}\label{sec:outdoor_results}

\newcommand{0.32 \textwidth}{0.32 \textwidth}
\begin{figure*}[t]
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_linear_trajectories}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_rectangular_trajectories}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_circular_trajectories}
    \end{subfigure}
    \\
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/linear_trajectories.pdf}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/rectangular_trajectories.pdf}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/circular_trajectories.pdf}
    \end{subfigure}
    \\
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/linear_distances.pdf}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/rectangular_distances.pdf}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32 \textwidth}
        \includegraphics[width=\textwidth]{results/migration/circular_distances.pdf}
    \end{subfigure}
    \\
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_linear_distances}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_rectangular_distances}
    \end{subfigure}
    \begin{subfigure}{0.32 \textwidth}
        \caption{}\label{fig:migration_circular_distances}
    \end{subfigure}
    \caption{%
        Results of vision-based flocking in outdoor environments.
        We show \tit{trajectories} (top) and inter-agent \tit{distances} (bottom) for three different migration scenarios: \tit{linear} (left), \tit{rectangular} (middle) and \tit{circular} (right) migration.
        The drones remain collision-free and cohesive using only local visual information, which is processed onboard in real-time.
        Centimeter-accurate ground-truth positions are obtained at $\SI{10}{\hertz}$ using RTK-enabled GNSS receivers mounted on the drones.
        These positions are not shared between the drones during the experiments and only serve to analyze the inter-agent distances.
    }\label{fig:migration}
\end{figure*}

\newcommand{0.24 \textwidth}{0.24 \textwidth}
\begin{figure*}[t]
    \begin{subfigure}{0.24 \textwidth}
        \includegraphics[width=\textwidth]{results/detections/easy.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.24 \textwidth}
        \includegraphics[width=\textwidth]{results/detections/difficult.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.24 \textwidth}
        \includegraphics[width=\textwidth]{results/detections/false_negative.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.24 \textwidth}
        \includegraphics[width=\textwidth]{results/detections/false_positive.png}
    \end{subfigure}
    \\
    \begin{subfigure}{0.24 \textwidth}
        \caption{High confidence}\label{fig:detection_easy}
    \end{subfigure}
    \begin{subfigure}{0.24 \textwidth}
        \caption{Low confidence}\label{fig:detection_difficult}
    \end{subfigure}
    \begin{subfigure}{0.24 \textwidth}
        \caption{False negative}\label{fig:detection_false_negative}
    \end{subfigure}
    \begin{subfigure}{0.24 \textwidth}
        \caption{False positive}\label{fig:detection_false_positive}
    \end{subfigure}
    \caption{%
        Qualitative examples of detection results categorized by confidence score (high and low) and type of error (false negative and false positive).
        High confidence (\ref{fig:detection_easy}): the drone is easily detected in front of the grass texture despite direct sunlight.
        Low confidence (\ref{fig:detection_difficult}): the drone blends in with the background due to the line features and lighting conditions.
        False negative (\ref{fig:detection_false_negative}): the right drone (yellow, manually labeled) can not be reliably distinguished from the background clutter.
        False positive (\ref{fig:detection_false_positive}): this type of spurious detection occurs relatively frequently but is filtered out by the multi-agent state tracker.
    }\label{fig:detections}
\end{figure*}

We report results for three different navigation scenarios: \tit{linear}, \tit{rectangular}, and \tit{circular} migration.
Before each flight, we place the drones at roughly $\SI{2.5}{\meter}$ distance from each other and wait for their RTK-GNSS receivers to converge to a fixed solution which provides centimeter-accurate absolute positions at $\SI{10}{\hertz}$.
These measurements are only used to provide a reliable ground-truth for the evaluation of the experiments.

After the RTK fix is obtained, we let all agents take off simultaneously and reach an altitude of $\SI{2}{\meter}$ before we let the vision-based flocking algorithm take over the control of their motion.
The agents are given the same list of migration points depending on the type of navigation scenario.
We switch from one migration point to the next as soon as an agent enters an acceptance radius of $r^\txt{acc} = \SI{3}{\meter}$.
As the list of migration points is exhausted, we repeat the procedure from the first waypoint.
We stop the experiment as soon as the battery level of one of the agents reaches a critical capacity of $15\%$.

The agents' altitude is individually regulated using a proportional controller to constrain their motion to a horizontal plane.
However, the planar constraint may be lifted by mounting additional cameras that point to the top and bottom, or by equipping the existing cameras with lenses that provide a larger field of view.

During the linear migration experiment, the agents fly between two waypoints that are located $\SI{10}{\meter}$ apart from each other (Fig.~\ref{fig:migration_linear_trajectories}).
Over a total flight duration of around $\SI{2.5}{\minute}$, the minimum inter-agent distance the agents reach is $\SI{1.82}{\meter}$ and the overall mean $\SI{2.42}{\meter}$ (Fig.~\ref{fig:migration_linear_distances}).

The rectangular migration experiment defines four waypoints that are located at the corners of a square with side length $\SI{10}{\meter}$ (Fig.~\ref{fig:migration_rectangular_trajectories}).
The total flight time is around $\SI{3.3}{\minute}$ and the overall minimum and mean inter-agent distances are $\SI{1.45}{\meter}$ and $\SI{2.36}{\meter}$, respectively (Fig.~\ref{fig:migration_rectangular_distances}).

Finally, the circular migration experiment leads the agents through a series of twelve waypoints that are linearly spaced around a circle with $\SI{10}{\meter}$ diameter (Fig.~\ref{fig:migration_circular_trajectories}).
During an overall flight time of around ${\SI{5}{\minute}}$, the agents remain collision-free while reaching a minimum inter-agent distance of $\SI{1.37}{\meter}$ and an overall mean distance of $\SI{2.32}{\meter}$ (Fig.~\ref{fig:migration_circular_distances}).

The above experiments are three representative flights taken from a total of $\SI{30}{\min}$ of collision-free experimental recordings.
The progressively lower inter-agent distances across linear, rectangular, and circular migration experiments can be explained by examining the distribution and/or density of waypoints.
During the rectangular migration experiment, the lowest inter-agent distances are reached close to the corners where directional changes occur.
In the case of the circular migration experiment, the larger number and density of waypoints have a cohesive effect since the agents simultaneously approach points that are more densely spaced.

During the experiments, a variety of objects that can potentially confuse the detector are present in the background (Fig.~\ref{fig:detections}).
In descending order of frequency, the most common objects aside from drones are trees, buildings, cars, people, fences, traffic signs, tables, and even dogs.
The other drones are generally detected despite background clutter and adverse lighting conditions (Fig.~\ref{fig:detection_easy} and \ref{fig:detection_difficult}).
False-negative detections occur most frequently when another drone is flying far away and/or in front of the ground control station, i.e. a camping table with the experimental equipment covered by a sun umbrella (Fig.~\ref{fig:detection_false_negative}).
False-positive detections appear most often in image regions that contain many line-like features since they resemble the mechanical design of the drone (Fig.~\ref{fig:detection_false_positive}).
They are mostly present for the duration of a single frame and the tracker can reliably reject them.
False tracks are occasionally created if false-positive detections are present for more than one frame.
However, the false tracks do not cause instabilities that lead to collisions.

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

We presented a vision-based detection and tracking algorithm that enables dense groups of drones to fly cohesively and without mutual collisions.
The proposed approach does not depend on visual markers or inter-agent communication and is thus suitable for flocking operation in GNSS-denied environments or in situations where wireless links are unreliable.
The approach is fully decentralized since each agent relies exclusively on onboard processing of local visual information to estimate the positions and velocities of neighboring drones.
The outdoor navigation experiments show that the system is robust to background clutter and enables collision-free flight even in demanding lighting conditions.
These experiments should be considered as minimal validation conditions of vision-based flocking algorithms without explicit communication or localization infrastructure.
Further experiments are needed to ensure the scalability of the proposed system to larger numbers of vision-based agents.

\section*{Acknowledgments}\label{sec:acknowledgments}

We thank Olexandr Gudozhnik for the contributions to the drone hardware, Vivek Ramachandran and Enrico Ajanic for the help with conducting outdoor experiments, as well as Enrica Soria, Valentin Wüest, and Davide Zambrano for the helpful discussions.

\balance

\bibliographystyle{IEEEtran}


\section{Construction of the Multi-LED VLC Model}\label{sec:multielement-architecture}

\section{Summary and Future Work}
\label{sec:summary}

In this paper, we have explored a novel mirror employment approach for performance improvement for an MEMD VLC system. In order to maximize the average throughput and illumination, we have formulated an optimization problem to obtain optimum mirror placement, power allocation and LED-user association. As the problem is NP-complete, we have proposed a two-stage solution of the optimization problem. As a part of the solution of the first stage, we have introduced several mirror placement approaches and analyzed their performance. To deal with heavy computation complexity in the second stage with the communication problem, we have presented multiple heuristic approaches to solve it with a detailed analysis on their performance. We have shown that up to threefold increase in average illumination and fourfold increase in average throughput can be achieved using a mirror placement approach.

In light of our proposed mirror placement approach, many promising future works are possible in the realm of indoor MEMD VLC networks. For example, it is worth exploring different mirror sizes and shapes which might yield further improvement in throughput or illumination level. For the communication problem, other heuristic algorithms may be explored to attain results that are closer to the optimal solution. It would also be worth observing how the proposed algorithms and mirror placement approaches perform with different system parameters (e.g., room shape, shape of the bulb) and whether there are any relationship between them. In our work, we considered a single bulb in a standard size room. Designing a MEMD VLC system for a much larger indoor environment with multiple such bulb is definitely a challenging and promising direction.

\section{Heuristic Solutions}

\section{Introduction}
Visible Light Communications (VLC) is a rapidly unfolding technology with noteworthy potential to meet the necessity of wireless access speeds. With the constantly increasing number of Internet-of-Things (IoT) devices and the humongous demand for wireless bandwidth, VLC solutions are of lofty worth.
Most of the related work in VLC have concentrated on diffuse optics \cite{2011-wu-power} and diversity combining \cite{2014-Tsiatmas-Optimum} for downloading a data stream to devices in a room. At the modulation level, OFDM \cite{2014-ijaz-Optical, aziz2019high} and Multi-Input-Multi-Output (MIMO) \cite{2014-Nuwanpriya-Angle, 2013-azhar-Gigabit} techniques were explored to increase the VLC link capacities. Unlike diffuse optics, Multi-element VLC (where multiple transmitters are involved in a VLC system design) allows several LED transmitters with narrow divergence angles to transmit different datastreams  simultaneously. Because of this capability, optical wireless communications researchers gained interest in multi-element VLC architectures ~\cite{mushfique2018software,mossaad2015visible,sifat} in the recent times. The multi-element VLC networks can offer increased aggregate throughput via simultaneous wireless links and attain higher spatial reuse. The downlink data transmission efficiency may be considerably enhanced by using multi-element VLC modules due to its light beam directionality where each transmitter, e.g., a Light Emitting Diode (LED), can be modulated with different data streams. Although, with this, the problem of maintaining a balanced illumination also arises because of the potential creation of dark spots in between the directional beams of the transmitters.

Most of the VLC literature can be categorized into four groups based on the number of transmitters and datastreams involved in the design. \textit{Single Element Single Datastream (SESD)} designs are based on the concept of a single LED transferring a single  datastream to a particular receiver. Most of the studies conducted in the early stages of VLC research were based on SESD VLC systems with a focus on improving the SINR \cite{wang2012performance,liu2014cellular}. \textit{Single Element Multi Datastream (SEMD)} is another form of VLC design where considerable amount of exploration has been done where a single transmitter is able to serve multiple receivers simultaneously or via time sharing methods like Time Division Multiple Access (TDMA) \cite{mmbaga2016performance,abdelhady2017resource} - most of these techniques involved diffuse optics. Apart from these two designs, other two can be categorized in multi-element VLC. \textit{Multi Element Single Datastream (MESD)} design allows multiple LEDs to send data towards one receiver which increases the received light intensity at it. Usage of more than one transmitter allows MESD to employ cooperative transmission~\cite{cooperative-little} to reduce Bit Error Rate (BER) and also various MIMO techniques by placing multiple photo-detectors at the receiver end where such demonstrated an improvement in data rate \cite{butala2013svd}. Despite having these benefits, MESD VLC designs cannot afford high directionality when transferring different data streams to different receiver at the same time, which is possible in \textit{Multi Element Multi Datastream (MEMD)}. These designs have surfaced in recent times, and various combinations in LED positioning are possible in these designs due to the flexibility in the number of elements and datastreams. Our proposed design is under this category where the advantage of light beam directionality is utilized in multi-element VLC system. More details on this categorization can be found in one of our earlier works \cite{mushfique2020optimization}.

Illumination uniformity is a critical element to be examined in a multi-element VLC architecture where uniform light distribution is essential to be maintained as each LED's transmit power is being tuned. Otherwise, inconsistent lighting might appear while the transmit powers of LEDs are being tuned for maximal SINR. It is very important to ensure an uniform lighting distribution in the room, while optimizing the LED assignment problem with a source power constraint on each LED. This type of multi-element architecture can improve the system performance as investigated in~\cite{chen2014improving,Combining_Receiver_2014}. Our work in this paper uses a totally different approach from the previous works in the literature, where we explore a novel approach of using mirror placement in multi-element multi-datastream VLC networks with narrow-angle LEDs in order to use the nLoS light beams for improving overall throughput and the evenness of lighting.

One of the main challenging tasks in  a multi-element VLC architecture is how to assign group of transmitters (i.e., LEDs) to each receiver in the room at the same time, i.e., serving all users in the room simultaneously and taking the illumination requirements into consideration. There have been several studies to find optimum LED arrangements in MEMD VLC systems to obtain different goals like evenly distributed lighting \cite{burton2012study}, reducing SINR fluctuation \cite{wang2012performance}, improving SINR while maintaining a certain illumination requirement \cite{liu2014cellular}, using cooperative beamforming to optimize total system throughput \cite{cen2019libeam} etc. Our idea is vastly different from these works for mainly two reasons - 1) Optimizing MEMD VLC system design with highly directional LEDs, and 2) Employing mirrors in the walls to use the nLoS component of the VLC channel looking to improve both SINR and illumination uniformity. There have not been much research on using mirrors in optimizing VLC networks.
In~\cite{park2017novel}, deployment of a double-sized mirror between the receiver photo-detectors is utilized to propose a mirror diversity receiver design which can help in reduction of the channel correlation by attaining two things - obstructing the reception of the light from a particular direction and enhancing the channel gain from another direction. The authors claim this mirror diversity receiver to be an encouraging approach to boost the performance of the VLC system, but the use and benefit from the mirror placement approach are still very limited unlike our work - where we have explored optimum mirror placement in all walls of the room, facilitating to broaden the scope of improvement for both lighting distribution and the SINR.

 In this paper, we propose a novel VLC framework  design with a hemispherical shaped multi-element bulb containing multiple LEDs and several mirrors mounted in the walls. This hemispherical shape allows high spatial reuse in addition to improving the illumination uniformity of the room. Light is emitted in different directions from the LEDs, and hence, the spherical multi-element bulb helps attaining an evenly scattered lighting. We further improve the average SINR and illumination uniformity of the system by formulating and solving an optimization problem, which optimizes the number of the mirrors as well as source power and assignment of the LEDs to the users to improve system performance. Apart from the the multi-element hemispherical bulb, this paper provides a framework using mirror placement for joint optimization of the LEDs' power and association to receivers to improve the aggregate data rate (of multiple data streams) and illumination uniformity in multi-element VLC networks for the first time to the best of authors' knowledge. The main contributions are as follows:

\begin{itemize}
	\item Investigating a hemispherical multi-element bulb architecture in downlink VLC transmission, where each LED can be assigned to a receiver to contribute towards data transmission or towards uniform illumination.
	\item Introducing a mirror placement approach that maximizes the ratio of SINR over the illumination uniformity of the system by formulating it as an optimization problem.
	\item Proposing a two-stage solution of the optimization problem to avoid the impracticality of changing mirror locations based on changes in the position of mobile users.
	\item Proposing multiple ways to place mirrors across the room walls, comparing and analyzing their performance with respect to the case of without any mirror placement.
	\item Proposing multiple low complexity heuristic algorithms to achieve efficient solutions for the second stage of the problem which focuses on maximizing the throughput of the system.
\end{itemize}

The rest of the paper is organized as follows:
Section~\ref{SystemModel} describes our spherical multi-element VLC architecture and the overall system model. The formulation of the proposed joint optimization problem along with the details of the two-stage solution is presented in Section \ref{ProblemFormulation}. Selected numerical simulation results are presented in Section \ref{sec:results}. Finally, we summerize our work in Section \ref{sec:summary}.

\section{Near Optimal Solution}

\section{Problem Formulation and Solution}\label{ProblemFormulation}
In a nutshell, we propose an optimization problem to optimize the combination of SINR and illumination uniformity based on mirror placement on the wall ($\xi_z$), LED-user assignment ($\epsilon_{mu}$) and source power for each of the $m$ LEDs ($P_m$). Although it is possible to re-optimize the objective value by updating $\epsilon_{mu}$ whenever a user changes its position, it is not feasible to do so for the mirrors. For this reason, we have divided the problem in two separate optimization problems. First, we find the optimal mirror placements and LED source powers that maximize the illumination uniformity, and then, we optimize the  LED-user association along with tuning the LED source powers to maximize the combination of SINR and illumination uniformity using the results from the first problem. This staging makes the optimal mirror placement more practical. In particular, the solution to the first problem will yield the best places to set up the grid mirrors so that high illumination uniformity is attained. Then, the second problem can be solved on-the-fly as users are moving in the room, yielding the best LED transmit powers and LED-user association depending on the user movements.

In this section we formulate two optimization problems defined as mirror design problem and communication problem.

\subsection{Design Problem}\label{design}

The main goal of this problem is to achieve the best illumination quality by optimizing not only the mirror placements but also the LEDs' transmit powers. Therefore, the mirror design optimization can be formulated as follows:
\begin{align}
	&\hspace{-0.5cm}\underset{\substack{\xi_ z, \chi_{mz} \in\{0,1\}, P_m}}{\text{maximize}} \quad  \underset{n}\min\left(\sum\limits_{m=1}^M \alpha_0 P_{m} H_{mn}\right)  \label{of_design}\\
	&\hspace{-0.5cm}\text{subject to \eqref{eqn:grid0} and:}\nonumber\\
	&\hspace{-0.5cm} P_{min} \leq P_m \leq \bar{P},  \quad \forall m,\label{power0}\\
	&\hspace{-0.5cm}
	\frac{\underset{n}\min\left(\sum\limits_{m=1}^M \alpha_0 P_{m} H_{mn}\right)}{\frac{1}{N}\sum\limits_{n=1}^N  \sum\limits_{m=1}^{M} \alpha_0 P_{m} H_{mn}} \geq \mu \label{unifor_constraint0}
	\\
	&\hspace{-0.5cm} \phi_2 \leq \sum\limits_{m=1}^M \alpha_0 P_{m} H_{mn} \leq \phi_1, \quad \forall n,  \label{minimum_illumination}
\end{align}
where $P_{min}$ and $\bar{P}$ are the minimum and maximum source power respectively that can be allotted to an LED and $\mu$ is the minimum acceptable illumination uniformity for an indoor setting~\cite{standardization2002lighting}. $\chi_{mz}$ is a dependent variable here, i.e., it is directly related to the values of $\xi_z$. $\phi_1$ and $\phi_2$ are the maximum and minimum levels of total illumination allowed at a particular sensing point. These are introduced to ensure that the room is illuminated at a minimum level and also it does not go beyond a level where it can be harmful to human eye.

We aim to solve the above design optimization problem optimally.
In order to do this, we first linearize the objective function \eqref{of_design} and constraint \eqref{minimum_illumination} by introducing a new decision variable $\phi$ as follows:
\begin{equation}\label{SNR}
\phi= \underset{n}\min\left(\sum\limits_{m=1}^M \alpha_0 P_{m} (h_{mn}^{\text{LoS}}+ \sum\limits_{z=1}^{Z} \chi_{mz} h_{mn}^{\text{NLoS}}(z))\right).
\end{equation}
Equation~\eqref{SNR} can be re-written as follows:
\begin{equation}\label{SNR2}
	\phi= \underset{n}\min\left(\sum\limits_{m=1}^M \alpha_0 P_{m} h_{mn}^{\text{LoS}}+ \sum\limits_{z=1}^{Z} \alpha_0 \rho_{mz} h_{mn}^{\text{NLoS}}(z)\right),
\end{equation}

\noindent where $\rho_{mz}=\chi_{mz} P_m$, which is introduced as another decision variable to linearize the product of binary variable $\chi_{mz}$ and real decision variables $P_m$.
Indeed maximizing the objective function given in \eqref{of_design} is equivalent to maximizing $\phi$ given that $\phi \leq \sum\limits_{m=1}^M \alpha_0 P_{m} H_{mn})$, $\forall n$.
With the introduction of $\rho_{mz}$, the following inequalities also have to be respected:
\begin{align}
&1) \quad P_m \geq \rho_{mz} \geq 0, \forall m, \forall z \notag\\
&2) \quad \rho_{mz} \geq \bar{P}_m \chi_{mz} - \bar{P}_m + P_m, \forall m, \forall z \notag\\
&3) \quad \rho_{mz} \leq \bar{P}_m \chi_{mz}, \forall m, \forall z.~\label{Linear}
\end{align}

The first two inequalities ensure that $\rho_{mz}$ value is between $\chi_{mz}$ and $P_m$. The third inequality guarantees that $\rho_{mz}=0$ if $\chi_{mz}=0$, and $\rho_{mz}=P_m$ if $\chi_{mz}=1$.
Therefore, the optimization problem can be reformulated as binary linear optimization problem as follows:
\begin{align}
	&\hspace{-0.5cm}\underset{\substack{\xi_ z, \chi_{mz}, P_m \\ \phi, \rho_{mz}}}{\text{maximize}} \quad  \phi  \label{of_design2}\\
	&\hspace{-0.5cm}\text{subject to \eqref{eqn:grid0}, \eqref{power0}, \eqref{Linear}, \eqref{minimum_illumination} and:}\nonumber\\
	&\hspace{-0.5cm} \phi \geq \frac{\mu}{N}\sum\limits_{n=1}^N  \sum\limits_{m=1}^{M} \alpha_0 \left(P_{m} h_{mn}^{\text{LoS}} + \sum\limits_{z=1}^{Z} \rho_{mz} h_{mn}^{\text{NLoS}}(z)\right) \label{unifor_constraintm2}\\
    &\hspace{-0.5cm} \phi \leq \sum\limits_{m=1}^M \alpha_0 \left(P_{m} h_{mn}^{\text{LoS}} + \sum\limits_{z=1}^{Z} \rho_{mz} h_{mn}^{\text{NLoS}}(z) \right), \forall n \label{minimum_illumination3}
\end{align}

It can be noticed that the solution for such binary linear optimization problem can be determined optimally using on-the-shelf software such as Gurobi/CVX interface~\cite{Gurobi}.

\subsection{Communication Problem}\label{communication}

After solving the mirror placement design optimization problem, we focus now on solving the communication phase aiming to maximize the SINR of the system to ensure the best possible data rate to the users. We choose to use Max-Min utility of the SINR. The approach of maximizing the total data rate which is known in the literature as Max C/I \cite{bifairness}, promotes users with favorable channel and interference conditions by allocating them most of the resources, whereas users suffering from higher propagation losses and/or interference levels will have very low data rates. Therefore, due to the unfairness of total sum data rate utility, the need for more fair utility metrics arises.  The Max-Min utilities are a family of utility functions attempting to maximize the minimum SINR in a network~\cite{Min-Max}. Our goal is to induce more fairness in the network by increasing the priority of users having lower SINR using the Max-Min utility. Thus, we formulate an optimization problem aiming to maximize the minimum SINR of all users by taking the association and illumination intensity constraints into consideration. This optimization problem can be expressed as

\begin{align}
&\hspace{-0.5cm}\underset{\substack{\epsilon_{mu} \in\{0,1\}, P_m \geq 0}}{\text{maximize}} \quad  \Gamma_{min} \label{of1}\\
&\hspace{-0.5cm}\text{subject to \eqref{eqn:eps1}, \eqref{power0} and \eqref{unifor_constraint0},}\nonumber
\end{align}

\noindent where $\Gamma_{min} = \underset{u}\min(\Gamma_u)$ is the minimum SINR among all users.
\\

\subsubsection{NP-Completeness}
As the number of LEDs $M$ being considered can be quite large, the number of ways to assign these LEDs to multiple users can be very large and it is practically infeasible to obtain an optimal solution as the users might move frequently inside the room, and the problem is needed to be solved again whenever any of the user coordinates is updated. In fact, this LED assignment problem in consideration is proven to be NP-Complete as elaborated in one of our earlier works~\cite{mushfique2020optimization}, so we propose three low complexity heuristic solutions to solve it.
\\
\subsubsection{Heuristic Approaches}
We describe each of our proposed heuristic approaches in this subsection.
\begin{figure}[t]
	\centering

	{\includegraphics[width=0.85\columnwidth]{Figs/Cone3.pdf}

	}

	\caption {A case of interference for an user-by-user approach such as \textit{NUA} or \textit{SSA-LED}: LEDs L1, L2 and L3 are assigned to user 1, 2 and 3 respectively. In this case, received power from L2 and L3 will be considered as interference at user 1 - since these LEDs are assigned to user 2 and user 3.}
	\label{fig:cone}
\end{figure}

\vspace{1mm}\noindent\emph{Nearest User Assignment (NUA)}:
We denote the resulting $P_m$ values from the solution of \eqref{of_design2} as intermediate power values $P^{prev}_{1}$, $P^{prev}_{2}$, .... $P^{prev}_{M}$ for $M$ LEDs, and  assign LED $m$ to the closest user to its beam projection if the user lies in the cone of LED $m$.
If there is one or no user in the cone, then we allocate $\min(P^{prev}_{m}(1+\tau),\bar{P})$ where $\tau$ is a measurement of how much we can deviate the allocated power to LED $m$ from its intermediate value $P^{prev}_{m}$. We use 0.1 as the value for $\tau$ in our simulations.
If there is more than one user in the cone of LED $m$, then we assign this LED to the user which is nearest to the center of the cone of the LED, but with fractional power, since there will be interference in this case as shown in Fig. \ref{fig:cone}. More details on this heuristic algorithm can be found in \cite{mushfique2020optimization} where it was introduced for the first time.

\begin{figure}[t]
	\centering

	{\includegraphics[width=0.85\columnwidth]{Figs/Cone3-SSA-Userpptx.pdf}

	}

	\caption {A case of conflict for an user-by-user approach such as \textit{SSA-User}: both User 1 and User 2 are getting strongest signal from LED L1, but since L1 is assigned to User 1 first, User 2 has to be assigned to L2, which provides the second strongest signal for User 2.}
	\label{fig:conessa}
\end{figure}

 	\vspace{1mm}\noindent\emph{Strongest Signal-based Assignment with User-first Approach (SSA-User)}:
	An LED-by-LED assignment approach is implemented in the \textit{NUA} algorithm. In contrast, in this approach, we do a user-by-user assignment. We inspect each user one by one and check the user in consideration is under how many LED beams, either via LoS or via NLoS. There are three possible scenarios -

	1) If there is no incoming LED beam towards this user, then we do not assign any LED to it.

	2) If there is only one incoming LED beam, then we follow similar approach to \textit{NUA} and assign the user with source power $\min(P^{prev}_{m}(1+\tau),\bar{P})$.

	3) If there are more than one incoming LED beams, then we take total channel model values  between each of these LEDs and this particular user into consideration. There are two possible cases in this scenario - i) There is no incoming LED  which is already assigned to another user, which means no possible interference - and ii) There is one or more LED(s) which are already assigned to other user(s), which creates interference.
	In the first case, since there is no chance of interference from any other user, we assign all the incoming LEDs to the user in consideration. However, a much more careful approach is needed for the second case, which is shown in Fig. \ref{fig:conessa}. Here we assign all the incoming LED beams, or simply `incoming LEDs', to the current user which are not already assigned to other user(s).
 We calculate the fraction, $\kappa$, as the ratio of maximum channel value of all the incoming LEDs including those also which are already assigned to other user(s) denoted as $H_{m_ju_i}$ assuming $m_j$ is the one among these LEDs with strongest channel value with the current user $u_i$ - and maximum of the channel values between all the incoming LEDs to the current user which are still unassigned, denoted as $H_{m_iu_i}$ assuming $m_i$ is the one among these LEDs with strongest channel value with $u_i$. It can be expressed as $\kappa = \frac{H_{m_ju_i}}{H_{m_iu_i}}$. Finally, we allot the LED $m_i$ with the source power of $\max ( P_m^{prev}(1-\kappa), P_m^{prev}((1-\tau))$.

	\vspace{1mm}\noindent\emph{Strongest Signal-based Assignment with LED-first Approach (SSA-LED)}: This approach is a blend of the previous two approaches - the scanning procedure is LED-by-LED like \textit{NUA} and the assignment is based on channel strength like \textit{SSA-User}. For each LED $m$ in the bulb, this approach checks whether there is any user inside it's coverage, either directly or via the mirror. We allocate $\min(P^{prev}_{m}(1+\tau),\bar{P})$ amount of source power when there is no or one user is covered by this LED. Now, if there are more than one user covered, we assign the current LED $m$ to the user $u_i$ with the strongest channel value only if the number of the users covered by this LED is less than or equal to 3 per cent of the total number of users. This is done to reduce the effect for interference when multiple users are being covered by the same LED. We give fraction of $P^{prev}_{m}$ to the current LED with a ratio $\kappa = \frac{H_{mu_j}}{H_{mu_i}}$, where $u_j$ is the user with the second strongest channel value with the current LED $m$. Thus the power allocated to LED m is $\max ( P_m^{prev}(1-\kappa), P_m^{prev}((1-\tau))$.

	\vspace{1mm}
	Each of the above described approaches has it's own advantages and disadvantages. \textit{NUA} has the simplest approach, although it is the worst-performing one which is elaborated further in Section \ref{sec:results}. \textit{SSA-User} is suitable for the scenarios where the users in the room are prioritized in some manner as it is an user-by-user approach which provides the best service to the user which is chosen first. But this is the most expensive approach considering the time complexity as many of the LEDs are needed to be examined multiple times if they cover more than one user. \textit{SSA-LED} can be considered as a balance between \textit{NUA} and \textit{SSA-User} since it is less complex than \textit{SSA-User} and obtains better performance than \textit{NUA}. It also provides fairness to all the users in the system.

\subsubsection{Computational Complexity}
We are assuming $M$ total LEDs and $U$ total users in our problem. In both \textit{NUA} and \textit{SSA-LED}, we check each of the LEDs to allocate source power and assign it to a user if needed. This effective runtime of checking this is $O(M)$. Calculating the channel values for all users will take $O(MU)$ time which is the most dominant component, so the overall runtime for both \textit{NUA} and \textit{SSA-LED} is $O(MU)$.

For \textit{SSA-User}, we take an user-by-user approach, and there may be cases when there is a conflict while assigning an LED to a particular user. In the worst case, runtime of this will be $O(M^2)$ when but on average this will be $O(MlogM)$ as this situation is similar as the Quick Sort algorithm. Combining this with the running time $O(MU)$ to calculate the channel values, we find the overall time complexity of the \textit{SSA-User} approach to be $O(M(U+logM))$.

\section{Simulation Results} \label{sec:results}

 \begin{table}[t]
 	\centering
 	\caption{Simulation Parameters}
 	\label{tab:parameters}
 	\begin{tabular} {
 			| p{4.0in} | p{1.99in} |}
 		\hline
 		\textbf{Parameter} & \textbf{Value} \\
 		\hline
 		Room Size & $6m \times 6m \times 3m$ \\
 		Radius of the hemispherical bulb, $R$ & $40cm$ \\
 		Number of LEDs on the bulb, $M$ & $391$ \\
 		Radius of an LED transmitter, $r_t$ & $1.5cm$ \\
 		Divergence angle of the LEDs, $\theta_d$ & $30^o$ \\
 		Maximum transmit power of the LEDs, $\bar{P}$ & $0.1W$ \\
 		Area of a PD receiver$^*$, $r_r$ & $5cm^2$ \\
 		Area of a sensing point, $r_s$ & $10cm^2$ \\
 		Number of users, $U$ & $2-12$ \\
 		Number of sensing points, $N$ & $100$ \\
 		Visibility, $V$ & $0.5km$ \\
 		Optical signal wavelength, $\lambda$ & $650nm$ \\
 		AWGN spectral density, $N_0$ & $2.5 \times 10^{-20} W/Hz$ \\
 		Modulation bandwidth, $B$ & $20MHz$ \\
 		Minimum uniformity, $\mu$ & $0.7$ \\
 		Minimum illumination, $\phi_2$ & $400$ $lux$ \\
 		Maximum illumination, $\phi_1$ & $600$ $lux$ \\
 		\hline
 	\end{tabular}
 	\thanks \vspace{3mm} {$^*$We assume that an array of PDs is used to attain a large receiver area.}
 \end{table}

In this section, we provide simulation results to study the performance of our MEMD VLC system model. Our aim is to make comparison between the three heuristics we proposed in terms of the \emph{minimum throughput} and the \emph{average throughput} of the system. To compare our proposed methods (\textit{NUA}, \textit{SSA-LED}, and \textit{SSA-User}), we use the same input parameters for all of them. We randomly place the users on the room floor with their receiver's FOV normal looking towards the ceiling. We report the average of the minimum throughout and illumination uniformity results among these randomly generated cases. To gain confidence in our results, we repeated the simulation experiments 100 times for all the results.

We use white LEDs with luminous efficacy $\alpha_0=169$  lumen/watt~\cite{pattison2018led}. The default values of the remaining input parameters used in our simulations are given in Table~\ref{tab:parameters}. For placing LEDs on the bulb, we followed the method in Section~\ref{SystemModel}. In particular, for a bulb with $R$ = 40cm radius, we place $M$ = 391 LEDs on 20 layers with $k_{1..20}$ = [1, 6, 12, 15, 19, 26, 30, 37, 43, 33, 30, 28, 25, 21, 16, 13, 11, 10, 9, 6].

\begin{figure}[t]
	\centering

	{\includegraphics[width=0.95\columnwidth]{code/results/heatmap_comparison.pdf}

	}

	\caption {Mirror placement heatmaps in four walls for different bulb radiuses and divergence angles.}
	\label{fig:heatmap}
\end{figure}

\subsection{Placement of Mirrors}
We solve the design problem in \eqref{of_design2} using different mirror placement approaches for different scenarios. For $R = 40 cm$ and  $\theta_d=30^o$, we solve the design problem considering three separate approaches : Mirror placement on 1) two adjacent walls, 2) two opposite walls and 3) all four walls. We define the four walls as $XZ$, $YZ$, $XZ+rs$ (the plane parallel with  $XZ$-plane) and $YZ+rs$ (the plane parallel with $YZ$-plane) planes, where $rs$ can be denoted as the room floor size. The mirror placement heatmaps from each of these scenarios are shown in Fig. \ref{fig:heatmap}. We can see that when we let the optimizer place mirrors on all four walls, it is difficult to establish any concrete pattern from the heatmaps although there are some similarities in placing the mirrors at lower altitudes. This happens because when the optimizer considers the sensing points near the corners of the room, there are several options to choose a mirror location from either of the two adjacent walls sharing the corner, and since there is no particular rule to do that, the mirrors are chosen randomly from those possible locations which could yield similar performance in terms of illumination level and uniformity, thus breaking the symmetry in the heatmaps. This is not the case when we allow the optimizer to place mirrors in two opposite walls as they do not share any corner, so we can see the same pattern of mirror placement in those two walls. Even in the case of two adjacent walls, the mirror placement follows the same pattern as the optimizer can avoid placing the mirror on the sides of the two walls which share the same corner and place most of the mirrors of the sides which do not share the same corner, as illustrated in the figure.

\subsection{Performance Comparison among Mirror Placement Approaches}
We compare the performance in terms of minimum throughput, average throughput and average illumination for the different mirror placement approaches. We also solve the design problem with no mirror placement at all to determine how much improvement is possible with mirror placement. The key findings are shown in Fig. \ref{fig:nua-comp} and \ref{fig:ssa-user-comp}. As the solution of the design problem is focused on optimizing the illumination, we see a noticeable three-fold improvement in average illumination level when mirror placement is applied, with the approach of placement in all walls having the most improvement aligning with our expectations. Although there is no significant improvement in minimum throughput after employing mirror placement, we can see up to four-fold increase in the average throughput values when comparing the four-wall case with the no mirror case. This fact indicates that placing mirror in all four walls can create more interference/noise for some users but at the same time they also provide better signal strength to most of the users resulting better average throughput. Adjacent walls and opposite walls cases also yield a solid improvement in both of average throughput and illumination level.

 \begin{figure*}[t]
 	\centering
 	\subfigure[Minimum throughput vs. No. of users. ]{
 		\includegraphics[width=0.31\columnwidth] {code/results/min_tp_comparison_mirror_placement_nua.pdf}
 		\label{fig:mintp-nua}
 	}
 	\subfigure[Average throughput vs. No. of users.]{
 		\includegraphics[width=0.31\columnwidth] {code/results/avg_tp_comparison_mirror_placement_nua.pdf}
 		\label{fig:avgtp-nua}
 	}
 	\subfigure[Average illumination vs. No. of users.]{
 		\includegraphics[width=0.31\columnwidth] {code/results/Illumination_comparison_mirror_placement_nua.pdf}
 		\label{fig:ilum-nua}
 	}
 	\caption{Comparison between different mirror placement approaches for \textit{NUA} with $M=391$ and $\theta_d=30^o$.}
 	\label{fig:nua-comp}
 \end{figure*}

\begin{figure*}[t]
 	\centering
 	\subfigure[Minimum throughput vs. No. of users. ]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/min_tp_comparison_mirror_placement_ssa_user.pdf}
 		\label{fig:mintp-ssa-user}
 	}
 	\subfigure[Average throughput vs. No. of users.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/avg_tp_comparison_mirror_placement_ssa_user.pdf}
 		\label{fig:avgtp-ssa-user}
 	}
 	\subfigure[Average illumination vs. No. of users.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/Illumination_comparison_mirror_placement_ssa_user.pdf}
 		\label{fig:ilum-ssa-user}
 	}
 	\caption{Comparison between different mirror placement approaches for \textit{SSA-User} with $M=391$ and $\theta_d=30^o$.}
 	\label{fig:ssa-user-comp}
 \end{figure*}

Overall, these comparisons provide valuable insights with regard to VLC network design. Since a much higher level of illumination
can be obtained on average from the same number of LEDs, either of the three mirror placement approaches can be seen as beneficial. A trade-off between system performance and ease of mirror deployment can be observed among the approaches. Among the three approaches, placing mirror in all four walls is the one which can yield maximum performance improvement in terms of illumination. In terms of average throughput, all the three approaches demonstrate better performance compared to the no mirror approach. However, in terms of the ease of mirror deployment, placing mirrors in the opposite walls is arguably the preferable approach as the mirror placement pattern is such that most of the chosen mirror grids are connected to each other - in that case the group of connected mirrors can be placed as a larger mirror - making the deployment easier. In rest of the simulations, we use the four walls approach to obtain intermediate results from the design problem which are used to solve the communication problem.

 \begin{figure*}[t]
 	\centering
 	\subfigure[Minimum throughput vs. No. of users. ]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/mintp.pdf}
 		\label{fig:comp_mintp}
 	}
 	\subfigure[Average throughput vs. No. of users.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/avgtp.pdf}
 		\label{fig:comp_avgtp}
 	}
 	\subfigure[Average illumination vs. No. of users.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/avgilum.pdf}
 		\label{fig:comp_ilum}
 	}
 	\caption{Comparison between \textit{NUA}, \textit{SSA-LED} and \textit{SSA-User} with $M=391$ and $\theta_d=30^o$.}
 	\label{fig:comp}
 \end{figure*}

\subsection{Comparison between NUA, SSA-LED, and SSA-User}

We compare our proposed heuristic approaches for solving the communication problem, \textit{NUA}, \textit{SSA-LED} and \textit{SSA-User}. We compare both minimum and average throughput among the users with average illumination level across the sensing points from 2 to 12 users for these three approaches. As seen in Fig. \ref{fig:comp}, minimum  and average throughput decreases for all three approaches with increased number of users. For \textit{SSA-User}, minimum throughput has lower value compared to \textit{NUA} because this approach highly prioritizes the users chosen first for LED assignments, which largely affects the users chosen at the end. \textit{SSA-LED} also has lower values compared to \textit{NUA}, as the restriction of assignment for some of the LEDs leads some users to starvation. But the average throughput is better for \textit{SSA-LED} than \textit{NUA} as that restriction ensures a considerable decline in interference as well. \textit{SSA-User} does even better in terms of the average throughput as the high quality service provided to the users assigned first comprehensively overpowers the low quality service provided to the users in the end. In Fig. \ref{fig:comp_ilum}, we can see that the average illumination remains steady for all heuristics with increasing number of users since changes in LED source power allocation for these heuristics are not substantial enough to have any adverse effect on the overall illumination level.

 \begin{figure*}[t]
 	\centering
 	\subfigure[Minimum throughput vs. Divergence angle. ]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/mintp_divangle.pdf}
 		\label{fig:ang_mintp}
 	}
 	\subfigure[Average throughput vs. Divergence angle.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/avgtp_divangle.pdf}
 		\label{fig:ang_avgtp}
 	}
 	\subfigure[Average illumination vs. Divergence angle.]{
 		\includegraphics[height=36mm,width=0.31\columnwidth] {code/results/ilum_divangle.pdf}
 		\label{fig:ang_ilum}
 	}
 	\caption{Comparison between \textit{NUA}, \textit{SSA-LED} and \textit{SSA-User} for different divergence angle with $6$ users. $M=391$ and $R=40 cm$.}
 	\label{fig:ang-comp}
 \end{figure*}

\subsection{Effect of Divergence Angle}

For each of our proposed heuristics, we plot minimum and average throughput with average illumination in the room for different divergence angles when there are six users in total, shown in Fig. \ref{fig:ang-comp}. With larger divergence angles, the LED beams get larger, thus there are more chances of a particular user to  get coverage of an LED, but at the same time the chance of interference increases too. Changes in divergence angle does not seem have much impact on minimum throughput for \textit{SSA-LED} and \textit{SSA-User} because of the randomness in user locations. Looking at the average throughput, \textit{SSA-LED} and \textit{SSA-User} handle the interference better than \textit{NUA} for increasing divergence angles as shown in Fig. \ref{fig:ang_avgtp}. Average illumination level is not hampered much with either heuristic in this case as well.

\section{System Model}
\label{SystemModel}

We consider a VLC downlink system model in an indoor setting consisting of a single hemispherical bulb containing $M$ directional LED transmitters in $L$ layers to serve $U$ mobile users. The transmit power of each LED is $P_m$ Watt, $\forall m=1,..,M$ as shown in Fig.~\ref{Layering}.

\begin{figure}[t]
	\centerline{\includegraphics[width=0.95\columnwidth]{Figs/LED-Placement.pdf}}
	\caption{Placement of transmitters in the multi-element bulb with 4 layers.}\label{Layering}
\end{figure}

The bulb is a hemispherical structure with two goals: i) achieve uniform lighting illumination within the room, and ii) provide high speed wireless download to mobile users. We assume that the LEDs in the bulb are placed in different layers in order to cover different locations in the room.
In order to improve the illumination uniformity and download speed to mobile users, we consider placing mirrors to the walls into the room's wall aiming to utilize NLoS beams from the LEDs. We divide each wall into an $X \times Y$ grid. We assume that in each grid element, one mirror can be placed at most. We call this mirror as a 'grid mirror'.
We anticipate that these grid mirrors in the walls will reflect the NLoS signals and thus potentially enhance the performance of the multi-element VLC network. The optimal placement of the mirrors, among other factors, will depend on the direction of the incident lights from the LEDs to the walls.
We consider that a user's Photo-detector (PD) gain depends on the light intensity from three main directions: 1) the LoS beams directly coming from LEDs, 2) strong reflections of the LoS light beams from the grid mirrors, and 3) weak-reflected light from the regular walls as shown in Fig~\ref{ChannelwithMirror}. Since there are four walls and each one is divided as $X \times Y$ sized grid, we visualize the scope of mirror placement as a $4X \times Y$ grid, where four walls are  assumed to be placed side-by-side. Each cell of the grid can be defined with index $z$ where $z = 1,2,3, ..... Z,$ here $Z=4XY$.

To compute the illumination uniformity in our simulated environment, we consider $N$ fixed sensing points uniformly distributed inside the room. The light intensity received at these points determine how uniform the lighting is inside the room. Although any of them can be placed at a place of interest, we assume that they are uniformly distributed to the room floor in a lattice placement pattern.

We also consider that each mobile user is equipped with a single PD to receive the downlink light beams.
It is assumed that the location of the mobile users in the room can  be predicted. On the other hand, we assume that the mobile users are equipped with an RF transmitter such as Wi-Fi for uplink transmission. In this work, we focus on downlink transmissions only and assume that the uplink transmission speed is not a bottleneck which is also the case for typical wireless access at an indoor setting.

\begin{figure}[t]
	\centerline{\includegraphics[width=0.95\columnwidth]{Figs/ChannelwithMirror.pdf}}
	\caption{System Model.}\label{ChannelwithMirror}
\end{figure}

\subsection{VLC Channel Model}

\subsubsection{LoS Channel}
The LoS channel model between LED $m$ and receiver node $l$ ($l \in \{u \text{ for users}, \newline
n \text{ for sensing points} \}$) can be expressed as~\cite{ghassemlooy2019optical}:
	\begin{equation}
	\small
	h_{ml}^{\text{LoS}}=\hspace{-.1cm}\left\{
	\hspace{-.1cm}\begin{array}{ll}
	\frac{A_l}{d^2_{ml}} Q_0(\varphi_{ml}) \cos(\phi_{ml}) \hspace{2.25cm}, &\hspace{-.2cm} 0 \leq \phi_{ml} \leq \phi_c \\
	\hspace{-.0cm} 0 \hspace{5.1cm} , &\hspace{-.2cm} \phi_{ml} \geq \phi_c
	\end{array}
	\right.
	\normalsize
	\end{equation}
where $A_l$ is the receiver node's PD area and $d_{ml}$ is the distance between LED $m$ and node $l$. $\varphi_{ml}$ and $\phi_{ml}$ are the irradiance and incidence angles, respectively, as shown in Fig.~\ref{ChannelwithMirror}.
$\phi_c$ is the FOV angle of the PD.
We have assumed that no optical filter is used.
$Q_0(\varphi_{ml})$ is the Lambertian radiant intensity and expressed as
\begin{equation}
Q_0(\varphi_{ml})= \frac{(q+1)}{2\pi} \cos^q(\varphi_{ml}),
\end{equation}

where $q=-\ln(2)/ \ln(\cos(\varphi_{1/2}))$ is the order of Lambertian emission and $\varphi_{1/2}$ is the transmitter semi-angle at half power.

\begin{figure}[t]
	\centerline{\includegraphics[width=0.95\columnwidth]{Figs/Channel-NLOS.pdf}}
	\caption{Comparison between an NLoS channel with mirror and LoS channel by using a mirror image of user $u$ to the wall from the same LED $m$.}\label{NLOS-fig}
\end{figure}

\subsubsection{NLoS Channel}
For simplicity, we consider only the strong reflections of the LoS light beams for computing the NLoS channel model.
In a scenario where the light from LED $m$ is coming to user $u$ through a grid mirror located at $\xi_{z}$ (i.e., $\xi_{z} = 1$), the NLoS channel model between LED $m$ and node $l$ can be expressed as~\cite{park2017novel}:
\begin{equation}
h_{ml}^{\text{NLoS}}(z)=\frac{\eta A_l}{{\hat{d}_{ml}(z)}^2} Q_0(\varphi_{ml}) \cos(\phi_{ml}) \label{nlos-eqn}
\end{equation}
where $\eta$ is the reflectivity of a grid mirror and $\hat{d}_{ml}(z)$ is the distance between $m$ and $l$ via the mirror located at index $z$ of the grid $\xi$. This is explained in Fig. \ref{NLOS-fig} where we compare the NLoS channel between LED $m$ and user $u$ with a channel between $m$ and $u'$ which is a mirror image of $u$ with respect to the wall. We can see that the distance between $m$ and $u'$ is same as the total distance between $m$ and $u$ via the mirror. Therefore, we can express the NLoS channel between $m$ and $u$  as the LoS channel between $m$ and $u'$ with the following in consideration - the NLoS channel strength will be reduced by a factor as there is a mirror in its path, and this factor is dependent on $\eta$, the reflectivity of the mirror. Same expression can be derived for the sensing points as well, thus we present the general expression for node $l$ in \eqref{nlos-eqn}.

\subsubsection{Total Channel}
We can express the total channel model as the combination of LoS and NLoS channels as below:
\begin{equation}
H_{ml} = h_{ml}^{\text{LoS}}+ \sum\limits_{z=1}^{Z} \chi_{mz} h_{ml}^{\text{NLoS}}(z)
\label{totalch}
\end{equation}
\hspace{1.9cm} such that
\begin{equation}
\begin{aligned}
\hspace{.5cm}\chi_{mz} = \xi_z, & \hspace{2cm}  if \hspace{2mm} z \in {\Upsilon}_{zm}, \\
\hspace{-.1cm}\chi_{mz} = 0, & \hspace{2cm} otherwise.
\end{aligned}
\label{eqn:grid0}
\end{equation}

To explain the above channel clearly, we define \emph{reflection area} of LED $m$ as the wall area being covered by the beam coming from LED $m$.

\begin{figure}[t]
	\centering

	{\includegraphics[width=0.95\columnwidth]{Figs/chi_explained.pdf}

	}

	\caption {Figure explaining different values of $\chi_{mz}$ in a room wall. $\chi_{mz} = \xi_z$ when grid position $z$ is inside the reflection area of LED $m$ (marked in orange) and 0 otherwise. Grid positions with mirrors are marked in light blue.}
	\label{fig:chi}
\end{figure}

$\Upsilon_{zm}$ in eq. \ref{eqn:grid0} is the set of $z$ indices corresponding to the reflection area of LED m. $\xi_z = 1$ if a grid mirror is placed at the cell located at index $z$ of the grid $\xi$ and $\xi_z = 0$ otherwise.
Further, $\chi_{mz}$ is a binary variable that expresses whether or not a NLoS component from LED m via grid location $z$ will be accumulated in the channel. Note that $\chi_{mz}$ cannot be 1 if the grid location $z$ is not within the reflection area of LED $m$. However, if the grid location $z$ is within the reflection area of LED $m$, then it can be 0 or 1 depending on the existence of a grid mirror at $z$ which is expressed by $\xi_z$.
Fig. \ref{fig:chi} explains different values for $\chi_{mz}$.

Although the above explained model seems a bit complex primarily because of the inclusion of an extra dependent variable $\chi_{mz}$, a designer of a VLC system can enjoy much more flexibility because of this. By bringing  different combinations of $\chi_{mz}$ into the overall search space, probability of finding a near optimal solution becomes much higher as we can control the contribution of every small fractions of each LED beam to either illumination only, or both illumination and communication.

\subsection{Illumination Uniformity}
The illumination intensity distribution across the room is an important factor to be considered in VLC system. One of common metric to measure the illumination intensity is the illumination uniformity. The illumination uniformity, $\vartheta$, can be defined as the ratio between the minimum and the average illumination intensity among all
N sensors and is given as~\cite{standardization2002lighting}:
\begin{equation}
\label{eqn:uni-exp}
\vartheta=\frac{\underset{n}\min\left(\sum\limits_{m=1}^M \alpha_0 P_{m} H_{mn}\right)}{\frac{1}{N}\sum\limits_{n=1}^N  \sum\limits_{m=1}^{M} \alpha_0 P_{m} H_{mn}}
\end{equation}
where $\alpha_0$ is the luminous efficiency that depends on the LED color wavelength, e.g. $\alpha_0=60$ lumen/watt for white LED~\cite{efficacywiki}. $\min(.)$ is the minimum function.

\subsection{LED-User Association}
We use a binary variable $\epsilon$ that indicates the association between LED $m$ and user $u$ which is given as follows:
\begin{equation}
\epsilon_{mu}=\hspace{-.1cm}\left\{
\begin{array}{ll}
\hspace{-.2cm}1, \hspace{-.3cm}& \hbox{if LED $m$ is associated with user $u$.} \\
\hspace{-.2cm}0, \hspace{-.3cm}& \hbox{otherwise.}
\end{array}
\right.
\end{equation}
We assume that user $u$ can be associated with multiple LEDs at the same time. However, it is assumed that each LED can be associated with one user at most during the same time. Therefore, the following conditions should be respected:
\begin{equation}
\sum\limits_{u=1}^U \epsilon_{mu} \leq 1, \forall m=1,..,M.
\label {eqn:eps1}
\end{equation}

\subsection{SINR Calculation}
We assume that each LED is either associated with one user or used for lighting only. Therefore, SINR at user $u$ can be expressed as ~\cite{8113471}
\begin{equation} \label{eqn:sinr-exp}
	\Gamma_u=\frac{ \left(\sum\limits_{m=1}^{M} \epsilon_{mu} H_{mu} P_{m} \right)^2}{N_0 B+ \sum\limits_{\substack{k=1 \\ k\neq u}}^U \left(\sum\limits_{m=1}^{M} \epsilon_{mk} H_{mk} P_{m}  \right)^2 }
\end{equation}

where $B$ and $N_0$ are the communication bandwidth and the spectral density of the Additive White Gaussian Noise (AWGN), respectively.