Mathieu Chalvidal
email: mathieuchalvid@brown.edu

Matthew Ricci

Rufin Vanrullen
email: rufin.vanrullen@cnrs.fr

Thomas Serre
email: thomasserre@brown.edu

GO WITH THE FLOW: ADAPTIVE CONTROL FOR NEURAL ODES

Despite their elegant formulation and lightweight memory cost, neural ordinary differential equations (NODEs) suffer from known representational limitations. In particular, the single flow learned by NODEs cannot express all homeomorphisms from a given data space to itself, and their static weight parameterization restricts the type of functions they can learn compared to discrete architectures with layer-dependent weights. Here, we describe a new module called neurallycontrolled ODE (N-CODE) designed to improve the expressivity of NODEs. The parameters of N-CODE modules are dynamic variables governed by a trainable map from initial or current activation state, resulting in forms of open-loop and closed-loop control, respectively. A single module is sufficient for learning a distribution on non-autonomous flows that adaptively drive neural representations. We provide theoretical and empirical evidence that N-CODE circumvents limitations of previous NODEs models and show how increased model expressivity manifests in several supervised and unsupervised learning problems. These favorable empirical results indicate the potential of using data-and activity-dependent plasticity in neural networks across numerous domains.

INTRODUCTION

The interpretation of artificial neural networks as continuous-time dynamical systems has led to both theoretical and practical advances in representation learning. According to this interpretation, the separate layers of a deep neural network are understood to be a discretization of a continuous-time operator so that, in effect, the net is infinitely deep. One important class of continuous-time models, neural ordinary differential equations (NODEs) [START_REF] Ricky | Neural ordinary differential equations[END_REF], have found natural applications in generative variational inference [START_REF] Grathwohl | Scalable reversible generative models with free-form continuous dynamics[END_REF] and physical modeling [START_REF] Köhler | Equivariant flows: sampling configurations for multibody systems with symmetric energies[END_REF][START_REF] Ruthotto | A machine learning framework for solving high-dimensional mean field game and mean field control problems[END_REF] because of their ability to take advantage of black-box differential equation solvers and correspondence to dynamical systems in nature.

Nevertheless, NODEs suffer from known representational limitations, which researchers have tried to alleviate either by lifting the NODE activation space to higher dimensions or by allowing the transition operator to change in time, making the system non-autonomous [START_REF] Dupont | Augmented neural odes[END_REF]. For example, [START_REF] Zhang | Approximation capabilities of neural odes and invertible residual networks[END_REF] showed that NODEs can arbitrarily approximate maps from R d to R if NODE dynamics operate with an additional time dimension in R d+1 and the system is affixed with an additional linear layer. The same authors showed that NODEs could approximate homeomorphisms from R d to itself if the dynamics were lifted to R 2d . Yet, the set of homeomorphisms from R d to itself is in fact quite a conservative function space from the perspective of representation learning, since these mappings preserve topological invariants of the data space, preventing them from "disentangling" data classes like those of the annulus data in Fig. 1 (lower left panel). In general, much remains to be understood about the continuous-time framework and its expressive capabilities. In this paper, we propose a new approach that we call neurally-controlled ODEs (N-CODE) designed to increase the expressivity of continuous-time neural nets by using tools from control theory. Whereas previous continuous-time methods learn a single, time-varying vector field for the whole input space, our system learns a family of vector fields parameterized by data. We do so by mapping the input space to a collection of control weights which interact with neural activity to optimally steer model dynamics. The implications of this new formulation are critical for model expressivity.

In particular, the transformation of the input space is no longer constrained to be a homeomorphism, since the flows associated with each datum are specifically adapted to that point. Consequently, our system can easily "tear" apart the two annulus classes in Fig. 1 (lower right panel) without directly lifting the data space to a higher dimension. Moreover, when control weights are allowed to vary in time, they can play the role of fast, plastic synapses which can adapt to dynamic model states and inputs.

The rest of the paper proceeds as follows. First, we will lay out the background for N-CODE and its technical formulation. Then, we will demonstrate its efficacy for supervised and unsupervised learning. In the supervised case, we show how N-CODE can classify data by learning to bifurcate its dynamics along class boundaries as well as memorize high-dimensional patterns in real-time using fast synapses. Then, we show how the flows learned by N-CODE can be used as latent representations in an unsupervised autoencoder, improving image generation over a base model.

BACKGROUND

Neural ODEs (NODEs) [START_REF] Ricky | Neural ordinary differential equations[END_REF] are dynamical systems of the form

dx dt = f (x, θ, t), (1)
where X is a space of features, θ ∈ Θ is a collection of learnable parameters, and f : X ×Θ×R → X is an equation of motion which we take to be differentiable on its whole domain. f defines a flow, i.e. a triple (X , R, Φ θ ) with Φ θ : X × R → X defined by

Φ θ (x(0), T ) = x(0) + T 0 f (x(t), θ, t)dt (2)
which relates an initial point x(0) to an orbit of points {x(t) = Φ θ (x(0), t), t ∈ R}. For a fixed T , the map x → Φ θ (x, T ) is a homeomorphism from X to itelf parametrized by θ.

Several properties of such flows make them appealing for machine learning. For example, the ODEs that govern such such flows can be solved with off-the shelf solvers and they can potentially model data irregularly sampled in time. Moreover, such flows are reversible maps by construction whose inverse is just the system integrated backward in time, Φ θ (., t) -1 = Φ θ (., -t). This property enables depth-constant memory cost of training thanks to the adjoint sensitivity method [START_REF] Gamkrelidze R V Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Gamkrelidze R V Pontryagin | The mathematical theory of optimal processes[END_REF] and the modeling of continuous-time generative normalizing flow algorithms (?).

Interestingly, discretizing Eq. 2 yields the recursive formulation of a residual network [START_REF] He | Deep residual learning for image recognition[END_REF] with a single residual operator f θ :

Φ ResNet (x(0), T ) = x(0) + T t=1 f θ (x t-1 ) (3)
In this sense, NODEs with a time-independent (autonomous) equation of motion, f , are the infinitely-deep limit of weight-tied residual networks. Relying on the fact that every nonautonomous dynamical system with state x ∈ R d is equivalent to an autonomous system on the extended state (x, t) ∈ R d+1 , Eq. 2 can also be used to model general, weight-untied residual networks. However it remains unclear how dependence of f in time should be modeled in practice and how their dynamics relate to their discrete counterparts with weights evolving freely across blocks through gradient descent.

3 N-CODE: LEARNING TO CONTROL DATA-DEPENDENT FLOWS General formulation -The main idea of N-CODE is to consider the parameters, θ(t), in Eq. 1 as control variables for the dynamical state, x(t). Model dynamics are then governed by a coupled system of equations on the extended state z(t) = (x(t), θ(t)). The initial value of the control weights, θ(0), is given by a mapping γ : X → Θ. Throughout, we assume that the initial time point is t = 0. The full trajectory of control weights, θ(t), is then output by a controller, g, given by another differentiable equation of motion g : Θ × X × R → Θ with initial condition γ(x 0 ). Given an initial point, x(0), we can solve the initial value problem (IVP) Given a initial state x(0), the module consists of an augmented dynamical system that couples activity state x and weights θ over time (red arrows). A mapping γ infers initial control weights θ0 defining an initial flow (open-loop control). This flow can potentially evolve in time as θ might be driven by a feedback signal from x (closed-loop, dotted line). This meta-parameterization of f can be trained with gradient descent by solving an augmented adjoint sensitivity system (blue arrows).

   dz dt = h(z, t) z(0) = z 0 =    dx dt , dθ dt = (f (x, θ, t), g(θ, x, t)) (x(0), θ(0)) = (x 0 , γ(x 0 )) , (4)
where h = (f, g). We may think of g and γ as a controller of the dynamical system with equation of motion, f . We model g and γ as neural networks parameterized by µ ∈ R nµ and use gradient descent techniques where the gradient can be computed by solving an adjoint problem (Pontryagin Lev Semyonovich ; Boltyanskii V G & F, 1962) that we describe in the next section. Our goal here is to use the metaparameterization in the space Θ to create richer dynamic behavior for a given f than directly optimizing fixed weights θ.

Well-posedness -If f and g are continuously differentiable with respect to x and θ and continuous with respect to t, then, for all initial conditions (x(0), θ(0)), there exists a unique solution z for Eq. 4 by the Cauchy-Lipschitz theorem. This result leads to the existence and uniqueness of the augmented flow

(X × Θ, R, Φ µ ) with Φ µ : (X × Θ) × R → X × Θ.
Moreover, considering the restriction of such a flow on X , we are now endowed with a universal approximator for at least the set of homeomorphisms on X given that this restriction constitutes a non-autonomous system. We discuss now how g and γ affect the evolution of the variable x(t), exhibiting two forms of control and noting how they relate to previous extensions of NODEs.

OPEN AND CLOSED-LOOP CONTROLLERS

If the controller outputs control weights as a function of the current state, x(t), then we say it is a closed-loop controller. Otherwise, it is an open-loop controller.

Open-loop control: First, we consider the effect of using only the mapping γ in Eq. 4 as a controller. Here, γ maps the input space X to Θ so that f is conditioned on x(0) but not necessarily on x(t) for t > 0. In other words, each initial value x(0) evolves according to its own learned flow (X , R, Φ γ(x(0)) ). This allows for trajectories to evolve more freely than within a single flow that must account for the whole data distribution and resolves the problem of non-intersecting orbits (see Figure 4). Recently, (Massaroli et al., 2020b) proposed a similar form of data-conditioned openloop control with extended state (x(t), x(0))). This is a version of our method in which γ is of the form γ(x) = θ(0) = [C : id] with C a constant vector, id is the identity function, and : denotes concatenation. Our open-loop formulation makes an architectural distinction between controller and dynamics and is consequently generalizable to the following closed-loop formulation.

d 2 x d 2 t = µ(1 -x 2 ) dx dt -x ⇐⇒ ẋ = µ(x -1 3 x 3 -θ) θ = 1 µ x
Figure 3: Example of dynamical control augmentation: The Van der Pol oscillator: As a second order differential equation, the dynamics cannot be approximated by a single 1-dimensional NODE with a constant control θ (degenerate solution in red). However, if the dynamics are decomposed into a planar system with a dynamic control variable θ(t), then the parameter µ can be adapted to fit a particular oscillatory regime (in black). This is a particular form of augmentation discussed in [START_REF] Dupont | Augmented neural odes[END_REF]; Norcliffe et al. ( 2020)) that showcases the benefit of using additional variables as evolving parameters of the dynamical system.

Closed-loop control: Defining a differentiable mapping, g, which outputs the time-dependent control weights θ(t) given the state of the variable x(t) yields a non-autonomous system on X (see Fig. 3). This can be seen as a specific dimension augmentation technique where additional variables correspond to the parameters θ.

However, contrary to appending extra dimensions which does not change the autonomous property of the system, this augmentation results in a module describing a time-varying transformation Φ θ(t) (x 0 , t). Note that this formulation generalizes the functional parameterization proposed in recent non-autonomous NODES systems [START_REF] Choromanski | An ode to an ode[END_REF]Massaroli et al., 2020b;?), since the evolution of θ(t) depends on x(t). Much like in the case of classical control theory, we hypothesized that the use of dynamic control weights would be of particular use in reacting to a nonstationary stimulus. We evaluate this hypothesis in section. 5.3.

The expressivity of N-CODE compared to other continuous-time neural networks is encapsulated in the following proposition. The result, proven in appendix, shows that both openloop and closed-loop control systems overcome NODEs' expressivity constraint with two distinct strategies, data-conditioning and statespace augmentation.

Proposition 1 -There exists a transformation φ : R d → R d which can be expressed by N-CODE but not by NODEs. In particular, φ is not a homeomorphism.

TRAINING

Loss function: Dynamics are evaluated according to a generalized loss function that integrates a cost over some interval [0, T ]: The loss in Eq. 5 is more general than in Chen et al. ( 2018) since can be any Lebesgue-measurable function of both states, x, and control parameters, θ. In particular, this includes penalties at discrete time points or over the entire trajectory (Massaroli et al., 2020b) but also regularizations on the weights or activations over the entire trajectory rather than the final state z(T ) of the system.

l(z) := T 0 (z(t), t)dt = T 0 (x(t), θ(t), t)dt (5)
In order to estimate a control function θ(t) that is optimal with respect to Eq. 5, we invoke Pontryagin's maximum principle (PMP) (Pontryagin Lev Semyonovich ; Boltyanskii V G & F, 1962), which only requires mild assumptions on the functional control space Θ and applies to functions f that are non-smooth in θ. The PMP gives necessary conditions on θ(t) at optimality via the augmented adjoint variable a(t). This quantity is the Jacobian of with respect to both x(t) and θ(t).

In the case of θ(t) being differentiable with respect to the meta-parameters µ, solving for the augmented adjoint state a(t) as in Chen et al. ( 2018) allows us to compute the gradient of the loss with respect to µ thanks to Theorem 1.

Theorem 1 -Augmented adjoint method: Given the IVP of equation 4 and for defined in equation 5, we have: In practice, we compute the Jacobian for this augmented dynamics with open source automatic differentiation libraries using Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF], enabling seamless integration of N-CODE modules in bigger architectures. We show an example of such modules in section B of Appendix.

∂l ∂µ = T 0 a(t) T ∂h ∂µ dt, such that a satisfies    da dt = -a T . ∂h ∂z - ∂ ∂z a(T ) = 0 (6)

EXPERIMENTS

REFLECTION AND CONCENTRIC ANNULI

We introduce our approach with the 1-and 2-dimensional problems discussed in [START_REF] Dupont | Augmented neural odes[END_REF]Massaroli et al., 2020b), consisting of learning either the reflection map ϕ(x) = -x or a linear classification boundary on a data-set of concentric annuli. Earlier results have shown that vanilla NODEs cannot learn these functions since NODEs preserve the topology of the data space, notably its linking number, leading to unstable and complex flows for entangled classes. However, both the open and closed loop formulations of N-CODE easily fit these functions, as shown in Figs 4 and5. The adaptive parameterization allows N-CODE to learn separate vector fields for each input, allowing for simpler paths that ease model convergence. Informally, we may think of the classification decisions as being mediated by a bifurcation parameterized by the data space.

Transiting across the true classification boundary switches the vector field from one pushing points to one side of the classification boundary to the other. Moreover, the construction of γ and g as a neural network guarantees a smooth variation of the flow with respect to the system state, which can potentially provide interesting regularization properties on the family learned by the model.

SUPERVISED IMAGE CLASSIFICATION

We now explore different data-dependent controls for continuous-time models on an image classification task. We tested the open-loop version of N-CODE against augmented NODEs [START_REF] Dupont | Augmented neural odes[END_REF] and the data-control model of (Massaroli et al., 2020b). We perform this experiment on the MNIST and CIFAR-10 image datasets. Our model consists of the dynamical system with equation of motion f expressed as a single convolutional block

(1 × 1 → 3 × 3 → 1 × 1) with 50 channels.
The control weights θ ∈ R K are a real-valued vector of K concatenated parameters instantiating the convolution kernels. The 10-class prediction vector is obtained as a linear transformation of the final state x(T ) of the dynamics with a softmax activation. We aimed at simplicity by defining the mapping γ as a two-layer perceptron with 10 hidden units. We varied how γ affects θ by defining different linear heads for each convolution in the block. (see C.2 of Appendix). We did not perform any fine-tuning, although we did augment the state dimension, akin to ANODE, by increasing the channel dimension by 10, which helped stabilize training. We show that the open-loop control outperforms competing models. Results are shown in Table . 1 with further details in Section C.2 of Appendix.

REAL-TIME PATTERN MEMORIZATION

To validate our intuition that closed-loop control is particularly valuable on problems involving flexible adaptation to non-stationary stimuli and to distinguish our system from related data-conditioned methods, we trained N-CODE on a dynamic, few-shot memorization task [START_REF] Miconi | Differentiable plasticity: training plastic neural networks with backpropagation[END_REF].

Here, the model is cued to quickly memorize sets of sequentially presented n-dimensional binary patterns and to reconstruct one of these patterns when exposed to a degraded version. For each presentation episode, we exposed the model sequentially to m = 2, 5, or 10 randomly generated n-length (n = 100, 1000) bit patterns with {-1, 1}-valued bits. Then, we presented a degraded version of one of the previous patterns in which half of the bits were zeroed, and we tasked the system with reconstructing this degraded pattern via the dynamic state, x(t). This is a difficult task since it requires the memorization of several high-dimensional data in real-time. [START_REF] Miconi | Differentiable plasticity: training plastic neural networks with backpropagation[END_REF] have shown that, although recurrent networks are theoretically able to solve this task, they struggle during learning, while models storing transient information in dynamic synapses tend to perform better. Following this result, we evaluated the performance of a closed-loop N-CODE model in which the output of the control function takes the form of dynamic synapses (Eq. 7). Here, the control parameters are a matrix of dynamic weights, θ(t) ∈ R m×m , which are governed by the outer product

µ xx T ,    dx dt , dθ dt = ρ (θx) , µ xx T (x(0), θ(0)) = (x stim , θ 0 ) ( 7
)
where µ is matrix of learned parameters, θ 0 is a learned initial condition, ρ is an element-wise sigmoidal non-linearity, is element-wise multiplication, and x stim is the presented stimulus. We have suppressed t for concision. We ran this experiment on the closed-loop version of N-CODE as well as an LSTM, continuous-time ODE-RNN [START_REF] Rubanova | Latent odes for irregularly-sampled time series[END_REF] and the data-controlled model of Massaroli et al. (2020b). All systems had the same number of learnable parameters, except for the LSTM which was much larger. We found that all considered models, except closed-loop N-CODE, struggled to learn this task (Fig. 6). Note that chance is .25 since half the bits are degraded.

A discrete LSTM endowed with a much larger hidden state of 5000 neurons learns marginally better than chance level. ODE-RNN and data control versions plateau at 17.5% error rate. For both n tested, N-CODE learned the problem not only much better, but strikingly faster.

UNSUPERVISED LEARNING : IMAGE AUTOENCODING WITH CONTROLLED FLOW

Finally, we use our control approach in an autoencoding model for low-dimensional representation learning in order to see if merely adding N-CODE to a simpler model can improve performance. Our idea is to use the solution to a system of linear ODEs to encode the model's latent representations. First, images from U ∈ R n are mapped by an encoder to a space of couplings, θ, by

γ µ : u → γ µ (u) = θ ∈ R m×m . (8)
The latent code is taken as the vector x(T ) acquired by solving the autonomous differential system presented in Eq. 9 from 0 to T starting from the initial condition (x 0 , γ µ (u)) with x 0 ∼ N (0, I):

   dx dt , dθ dt = (θx, 0) (x(0), θ(0)) = (x 0 , γ µ (u)) (9)
This model, which we call autoN-CODE, can be interpreted as a hybrid version of autoencoders and generative normalizing flows, where the latent generative flow is data-dependent and parameterized by the encoder output.

We found that a simple linear system was already more expressive than a discrete linear layer in shaping the latent representation. Since the output dimension of the encoder grows quadratically with the desired latent code dimension of the bottleneck, we adopted a sparse prediction strategy, where as few as two elements of each row of θ were non-zero, making our model match exactly the number of parameters of a VAE with same architecture. We believe that the increase in representational power over a vanilla autoencoder, as shown by improvements in Fréchet Inception Distance (FID; [START_REF] Heusel | Gans trained by a two time-scale update rule converge to a local nash equilibrium[END_REF]) and more disentangled latent representations (see A.7), despite these simplifying assumptions, shows the potential of this approach for representation learning (Fig. 2).

Training details can be found in Sec. C.4.

Here we measure the marginal improvement acquired by adding latent flows to a baseline VAE. In order to endow our deterministic model with the ability to generate image samples, we perform, similarly to [START_REF] Ghosh | From variational to deterministic autoencoders[END_REF], a post hoc density estimation of the distribution of the latent trajectories' final state. Namely, we employ a gaussian mixture model fit by expectation-maximization and explore the effect of increasing the number of components in the mixture. We report in Table 2, contrary to [START_REF] Ghosh | From variational to deterministic autoencoders[END_REF], a substantial decrease in FID of our sampled images to the image test set with an increasing number of components in the mixture (see Fig. 7 for latent evolution and Fig. 8 Considering the identical architecture of the encoder and decoder to a vanilla model, these results suggest that the optimal encoding control formulation produced a structural change in the latent manifold organization. (See section C.7 in Appendix). Finally, we note that this model offers several original possibilities to explore, such as sampling from the control coefficients, or adding a regularization on the dynamic evolution as in [START_REF] Finlay | How to train your neural ode: the world of jacobian and kinetic regularization[END_REF], which we leave for future work.

RELATED WORK

Neural ODEs -This work builds on the recent development of Neural ODEs, which have shown great promise for building generative normalizing flows [START_REF] Grathwohl | Scalable reversible generative models with free-form continuous dynamics[END_REF] used for modeling molecular interactions [START_REF] Köhler | Equivariant flows: sampling configurations for multibody systems with symmetric energies[END_REF] or in mean field games [START_REF] Ruthotto | A machine learning framework for solving high-dimensional mean field game and mean field control problems[END_REF]. More recent work has focused on the topological properties of these models [START_REF] Dupont | Augmented neural odes[END_REF][START_REF] Zhang | Approximation capabilities of neural odes and invertible residual networks[END_REF], introducing time-dependent parameterization (Zhang et al., 2019b;Massaroli et al., 2020b;[START_REF] Choromanski | An ode to an ode[END_REF], developing novel regularization methods using optimal transport or stochastic perturbations [START_REF] Finlay | How to train your neural ode: the world of jacobian and kinetic regularization[END_REF][START_REF] Viktor Oganesyan | Stochasticity in neural odes: An empirical study[END_REF] and adapting them to stochastic differential equations [START_REF] Tzen | Neural stochastic differential equations: Deep latent gaussian models in the diffusion limit[END_REF].

Optimal Control -Several authors have interpreted deep learning optimization as an optimal control problem [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF][START_REF] Benning | Deep learning as optimal control problems: models and numerical methods[END_REF][START_REF] Liu | Deep learning theory review: An optimal control and dynamical systems perspective[END_REF]Zhang et al., 2019a;[START_REF] Jacob | Robust deep learning as optimal control: Insights and convergence guarantees[END_REF] providing strong error control and yielding convergence results for alternative algorithms to standard gradient-based learning methods. Of particular interest is the mean-field optimal control formulation of (E et al., 2018), which notes the dependence of a unique control variable governing the network dynamics on the whole data population.

Hypernetworks -Finally, our system is related to network architectures with adaptive weights such as hypernetworks [START_REF] Ha | [END_REF], dynamic filters [START_REF] Jia | Dynamic filter networks[END_REF] and spatial transformers [START_REF] Jaderberg | Spatial transformer networks[END_REF]. Though not explicitly formulated as neural optimal control, these approaches effectively implement a form of modulation of neural network activity as a function of input data and activation state, resulting in compact and expressive models. These methods demonstrate, in our opinion, the significant untapped potential value in developing dynamically controlled modules for deep learning.

CONCLUSION

In this work, we have presented an original control formulation for continuous-time neural feature transformations. We have shown that it is possible to dynamically shape the trajectories of the transformation module applied to the data by augmenting the network with a trained control mechanism, and further demonstrated that this can be applied in the context of supervised and unsupervised representation learning. In future work, we would like to investigate the robustness and generalization properties of such controlled models as well as their similarities with fast-synaptic modulation systems observed in neuroscience, and test this on natural applications such as recurrent neural networks and robotics. An additional avenue for further research is the connection between our system and the theory of bifurcations in dynamical systems and neuroscience.

APPENDIX A THEORETICAL RESULTS

A.1 PROOF OF PROPOSITION 1

In order to show that the set of functions defined by N-CODE are not necessarily homeomorphisms on X , we show that they are not generally injective.

Open-loop: Consider the one-dimensional N-CODE system Φ : (X , R) → X with equation of motion f (x, θ) = -θ where θ = γ(x(0)) = x(0). This system has solution

Φ(x(0), T ) = x(0) + f (x, θ, t) dt = x(0) -x(0)T. (10)
However, this implies Φ(0, 1) = Φ(1, 1) = 0, (11) giving the result.

Closed-loop: Similarly, we can show that the following 2D oscillatory system:

ẋ(t) = θ(t) θ(t) = -x(t) =⇒ x(t) = α cos(t) + β sin(t) =⇒ x(t) = x(0) cos(t) if θ(0) = 0 (12) But Φ(x, π/2) = 0, so x → Φ(x, t) = x(t) is not injective.
This shows that the two control forms, open and closed, are not restricted to learn homemorphisms on the space X .

A.2 PROOF OF THEOREM 1

Proof. The proof is inspired from Massaroli et al. (2020a). We place our analysis in the Banach space R n . Since and h are continuously differentiable, let us form the Lagrangian functional L where a is an element of the dual space of z = x θ :

L(z, a, t) := T 0 [ -a, ∂z ∂t -h ]dt (13)
The constraint a(t), ∂z ∂t -h(z, t) is always active on the admissible set by construction of z in equation 4, such that we have ∀z, a, ∂L ∂µ (z, a, t) = ∂l ∂µ (z, a, t). Moreover, integrating the left part of the integral in equation 13 gives:

T 0 a, ∂z ∂t dt = a, z T 0 - T 0 ∂a ∂t , z dt (14)
Now, given that -da dt , z is differentiable in µ for any t ∈ [0, T ], Leibniz integral rule allows to write

∂l ∂µ = T 0 ∂ ∂µ [ + da dt , z + a, h ]dt + a(T ), ∂z(T ) ∂µ -a(0), ∂z (0)
The last equation can be reordered as:

∂l ∂µ = T 0 da dt -a ∂h ∂z - ∂ ∂z , ∂z ∂µ dt + T 0 a, ∂h ∂µ dt + a(T ), ∂z(T ) ∂µ (18)
Posing that a(T

) = O |X |+|Θ| and da dt -a ∂h ∂z - ∂ ∂z , ∂z ∂µ = 0, the result follows.
Remark 1: This result allows to compute gradient estimate assuming that µ directly parametrizes h such that the gradient ∂h ∂µ are straightforward to compute. However, in the case of optimizing the initial mapping γ that infer the variable θ , this result can be combined with usual chain rule to estimate the gradient:

∂l ∂µ γ = ∂θ ∂γ ∂l ∂θ (19)
B IMPLEMENTATION with k the number of filters, c the number of channels (1 for MNIST and 3 for CIFAR10) and a the augmentation channel [START_REF] Dupont | Augmented neural odes[END_REF]. The open-loop controller is a very simple linear transformation of the data into 10 hidden units followed by multiple heads outputting the vector of weights for each convolution kernel. In order to examine the effect of the size of dimensionality of the equation of motion on learnability, we tried different parameterizations for f :

• NCODE : The three kernels are conditioned by the controller γ.

• NCODE 3 × 3: The inner 3x3 kernels are dynamically conditioned.

• NODE DC : This corresponds to the form of data-control proposed by (Massaroli et al., 2020b) where the image data is concatenated with state x(t) at each evaluation of the first 1x1 convolution kernel.

All parameters not adaptively controlled are simply learned as fixed parameters. For all models we used the adaptative solver Dormund-Prince with tolerance of 1e-3. The results show that N-CODE models achieve higher accuracy in fewer epochs and overall lower loss (Fig. A.2, top panels). However, the simplicity of the mapping γ in this experiment and the absence of regularization of the induced map Φ γ [START_REF] Finlay | How to train your neural ode: the world of jacobian and kinetic regularization[END_REF] yields a complex flow requiring a higher total number of function evaluations (Fig. A.2, lower panels). This potentially prevents NCODE from achieving even better accuracy. Interestingly, the results for the partially controlled systems suggest that gradient descent accommodates the learning of hybrid modules f with both static and dynamic convolution parameters. These results advocate for the exploration of specific closed-loop formulations and regularization of control for convolution layers which we leave for future work.

C.3 PATTERN MEMORIZATION

This experiment is an update of [START_REF] Miconi | Differentiable plasticity: training plastic neural networks with backpropagation[END_REF], which showed that using a simple Hebbian update rule between weights endowed recurrent networks with fast memory properties. The base model is a single-cell recurrent network with all-to-all coupling, θ(t). The considered models are

Open-loop control: We adopt the implementation of (Massaroli et al., 2020b) by appending the presented pattern x stim to the dynamic variable x(t) such that θ(t) = θ is linear layer L(R 2×N , R N ).

Closed-loop control: The weights are fully dynamic and the influence of the plastic evolution is tuned by learning the components µ that apply an element-wise gain for every weight velocity, according to Eq. 7.

ODE-RNN: Here, a static-weight NODE models the evolution of a continuous-time single-cell recurrent network whose hidden state is carried over from one presentation to the next. A supplementary linear read-out L(R N , R N ) outputs the model's guess and balances the number of parameters with other models.

LSTM: A vanilla LSTM cell [START_REF] Hochreiter | Long short-term memory[END_REF]. In coherence with [START_REF] Miconi | Differentiable plasticity: training plastic neural networks with backpropagation[END_REF] results, the hidden state dimension need to be greatly increased for the module to start memorizing. We tested a 5000-dimensional hidden state for the reported results.

For each model, we learn the model weights with gradient descent using an Adam optimizer with a learning rate λ = 3e -4 and use the L 2 distance between the reconstruction and objective pattern as our objective function. For continuous models, each episode consists of a sequential presentation of each pattern for 0.5 sec followed by a query time of 0.5 sec. For the LSTM, we adopted the setting of [START_REF] Miconi | Differentiable plasticity: training plastic neural networks with backpropagation[END_REF] where the sequence of presentation is discretized into 5 time-steps of presentations. We tested for episodes of 3,5 and 10 patterns with 2 presentations in random order and different proportions of degradation (0.5,0.7 and 0.9). Static-weight model performance rapidly deteriorated in more challenging settings, whereas N-CODE reconstruction converged to a residual error below 1% in all cases. MNIST CIFAR-10 CelebA

x ∈ R 3×28×28 x ∈ R 3×32×32 x ∈

C.4 AUTON-CODE ARCHITECTURES

Adapting the models used by [START_REF] Tolstikhin | Wasserstein autoencoders[END_REF] and [START_REF] Ghosh | From variational to deterministic autoencoders[END_REF], we use a latent space dimension of 25 for MNIST, 128 for CIFAR-10 and 64 for CelebA. All convolutions and transposed convolutions have a filter size of 4×4 for MNIST and CIFAR-10 and 5×5 for CELEBA. We apply batch normalization to all layers. They all have a stride of size 2 except for the last convolutional layer in the decoder. We use Relu non-linear activation and batch normalisation at the end of every convolution filter. Official train and test splits are used for the three datasets. For training, we use a mini-batch size of 64 in MNIST and CIFAR and 16 for CelebA in AutoN-CODE. (64 for control models.) All models are trained for a maximum of 50 epochs on MNIST and CIFAR and 40 epochs on CelebA. We make no additional change to the decoder. We train the parameters of the encoder and decoder module for minimizing the mean-squared error (MSE) on CIFAR-10 and CelebA [START_REF] Liu | Deep learning face attributes in the wild[END_REF] or alternatively the Kullback-Leibler divergence between the data distribution and the output of the decoder for MNIST (formulas in Appendix.) Gradient descent is performed for 50 epochs with the Adam optimizer (Kingma & Ba, 2014) with learning rate λ = 1e -3 reduced by half every time the loss plateaus. All experiments are run on a single GPU GeForce Titan X with 12 GB of memory.

C.5 VISUALIZATION OF LATENT CODE DYNAMICAL EVOLUTION

C.6 EXPLORING SAMPLING DISTRIBUTIONS

For random sampling, we train the VAE with a N (0, I) prior. For the deterministic models, samples are drawn from a mixture of multivariate gaussian distributions fit using the testing set embeddings. The distribution is obtained through expectation-maximization [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] with one sin- gle k-means initialization, tolerance 1e -3 and run for at most 100 iterations. We compute the FID using 10k generated samples evaluated against the whole test set for all FID evaluations, using the standard 2048-dimensional final feature vector of an Inception V3 following [START_REF] Heusel | Gans trained by a two time-scale update rule converge to a local nash equilibrium[END_REF] implementation.

The results show that 10K components (matching the number of data points in the test sets) naturally overfits on the data distribution. Generated images display marginal changes compared to test images. However, 1000 components does not, showing that our AutoN-CODE sampling strategy mediates a trade-off between sample quality and generalization of images. We alternatively tested non-parametric kernel density estimation with varying kernel variance to replace our initial sampling strategy. We report similar results to the gaussian mixture experiment with an overall lower FID of AutoN-CODE for small variance kernels. As the fitting distribution becomes very rough (σ ≈ 5), the generated image quality is highly deteriorated.(see Figure A6).

Figure 1 :

 1 Figure 1: Vector fields for continuous-time neural networks. Integral curves with arrows show the trajectories of data points under the influence of the network. Top left: Standard NODEs learn a single time-independent flow (in black) that must account for the whole data space. Top right: N-CODE learns a family of vector fields (red vs yellow vs blue), enabling the system to flexibly adjust the trajectory for every data point. Bottom left: Trained NODE trajectories of initial values in a data set of concentric annuli, colored green and yellow. The NODE transformation is a homeomorphism on the data space and cannot separate the classes as a result. Colored points are the initial state of the dynamics, and black points are the final state. Bottom right: Corresponding flows for N-CODE which easily separate the classes. Transiting from the inner to outer annulus effects a bifurcation which linearly separates the data.

Figure 2 :

 2 Figure 2: Diagram of a general N-CODE module:

Figure 4 :

 4 Figure 4: Trajectories over time for three types of continuous-time neural networks learning the 1-dimensional reflection map ϕ(x) = -x. Left: (NODE) Irrespective of the form of f , a NODE module cannot learn such a function as trajectories cannot intersect[START_REF] Dupont | Augmented neural odes[END_REF]. Middle: (Open-loop N-CODE). The model is able to learn a family of vector fields by controlling a single parameter of f conditioned on x(0). Right: (Closed-loop N-CODE) With a fixed initialization, x(0), the controller model still learns a deformation of the vector field to learn ϕ. Note that the vector field is time-varying in this case, contrary to the two others versions.

  where ∂h ∂z is the Jacobian of h with respect to z: ∂h ∂z =

Figure 5 :

 5 Figure 5: The concentric annuli problem. A classifier must separate two classes in R 2 , a central disk and an annulus that encircles it. Left: Soft decision boundaries (blue to white) for NODE (First), N-CODE open-loop (Second) and closed-loop (First) models. Right: Training curve for the three models.

Figure 6 :

 6 Figure 6: Left: Diagram of an episode of the few-shot memorization task Miconi et al. (2018). Three, 5-length -1, 1-valued bit patterns (n = 5 here for visualization) are presented to the system at regular intervals. At query time, tq, a degraded version of one of the patterns is presented and the task of the model is to complete this pattern. Right: Performance of a closed-loop N-CODE, an LSTM, a continuous-time ODE-RNN Rubanova et al. (2019) and the data-conditioned NODE of Massaroli et al. (2020b) averaged over 3 runs. All models besides N-CODE plateau at high values (chance is .25). N-CODE, on the other hand, learns immediately and with extremely high accuracy across all n. (see Appendix)

Figure 7 :

 7 Figure 7: Left to right: Temporal evolutions of an image reconstruction along controlled orbits of AutoN-CODE for representative CIFAR-10 and CelabA test images. Last column: ground truth image.

Figure 8 :

 8 Figure 8: Left Reconstructions of random test images from CelebA for different autoencoding models. Right: Random samples from the latent space. Deterministic models are fit with 100components gaussian mixture model.

  Figure A.1: Time-varying flows for two points in different classes (inner: top row; outer: bottom row) of the annuli dataset. Time advances from left to right.

Figure A. 2 :

 2 Figure A.2: Accuracy, Train and test losses over training interations and Number of function evaluations (NFE) vs loss for several versions of NCODE on CIFAR-10.

Figure A. 3 :

 3 Figure A.3: Error rate of N-CODE on the 1000-bit reconstruction task as a function of episodes for (Left) different proportions of degraded bits and (Right) different number of patterns presented in one episode. In all cases, the model learns in several hundreds episodes.

Figure A. 4 :

 4 Figure A.4: Reconstructions of the image along the controlled orbits of AutoN-CODE for MNIST. The temporal dimension reads from left to right. Last column: ground truth image

Figure A. 8 :

 8 Figure A.8: Nearest neighbors search of random samples from the gaussian mixture with (Left) 1000 components (Left) and (Right) 10K components in the testing set of CIFAR-10. The sampled latents show overfiting on the test set in the highest case.

Figure A. 9 :

 9 Figure A.9: (Left) Evolution of FID between test set and generated samples as a function of the gaussian kernel variance used to perform kernel estimation of the the latent code distribution. (Right) Nearest neighbor search of random samples from the latent distribution fit by gaussian kernel estimation with variance σ = 2 in the testing set.

Table 2 :

 2 Fréchet

	for sampled examples).
	4.5 IMAGE GENERATION
	FID (↓)	CIFAR-10 CelebA-(64)
	VAE † (Kingma & Welling, 2013) 105.45	68.07
	Auto-encoder † (100 comp.)	73.24	63.11
	Auto-encoder † (1000 comp.)	55.55	55.60
	AutoNCODE (100 comp.)	34.90	55.27
	AutoNCODE (1000 comp.)	24.19	52.47

Inception Distance (FID) for several recent architectures (lower is better). "Components" are the numbers of components used in the gaussian mixture estimation of the latent distribution.

  R 3×64×64 Encoder: Flatten Conv 256 → BN → Relu Conv 256 → BN → Relu FC 400 → Relu Conv 512 → BN → Relu Conv 512 → BN → Relu FC 25 → Norm Conv 1024 → BN → Relu Conv 1024 → BN → Relu FC 128 * 2 → Norm FC 64 * 10 → Norm Norm FC 784 → Sigmoid ConvT 512 → BN → Relu ConvT 512 → BN → Relu ConvT 256 → BN → Relu ConvT 256 → BN → Relu ConvT 3 ConvT 3 Table A.1: Model architectures for the different data sets tested. FC n and Conv n represent, respectively, fully connected and convolutional layers with n output/filters. We apply a component-wise normalisation of the control components which proved crucial for good performance of the model. The dynamic is run on the time segment [0,10] which empirically yields good results.

	Latent dynamics:	ODE[0,10]	ODE[0,10]	ODE[0,10]
	Decoder:	FC 400 → Relu	FC 1024 * 8 * 8 → Norm	FC 1024 * 8 * 8 →

ACKNOWLEDGMENTS AND DISCLOSURE OF FUNDING

This work was funded by the ANR-3IA Artificial and Natural Intelligence Toulouse Institute (ANR-19-PI3A-0004).

Additional support provided by ONR (N00014-19-1-2029), NSF (IIS-1912280, OSCI-DEEP ANR grant (ANR-19-NEUC-0004) and the Center for Computation and Visualization (CCV). We acknowledge the Cloud TPU hardware resources that Google made available via the TensorFlow Research Cloud (TFRC) program as well as computing hardware supported by NIH Office of the Director (S10OD025181).

to fit the latent distribution. We also tested a vanilla autoencoder with twice the capacity (236-dimensional), exactly matching the number of parameters of our model. Interestingly, this model performs worst in a low component regime. We interpret this a as a manifestation of the "curse of dimensionality", as the latent example population becomes less dense in this augmented space, making the naive gaussian mixture fitting less accurate for fitting the latent distribution. Right: Samples of generated images for different number of mixture components. The shape and details becomes clearer with increasing components.


Mahmoud Wafik Eltokhey

Mohammad-Ali Khalighi

Zabih Ghassemlooy

UAV Location Optimization in MISO ZF Pre-coded VLC Networks

Keywords: Visible-light communications, unmanned-aerial vehicles, ZF pre-coding, particle swarm optimization

come

I. INTRODUCTION

Visible-light communication (VLC) systems exploit the light-emitting diode (LED)-based luminaires to provide simultaneous lighting and data communications [START_REF] Ghassemlooy | Visible Light Communications: Theory and Applications[END_REF]. The growing use of LEDs for illuminations in unmanned aerial vehicles (UAVs) (i.e., drones) has opened up opportunities for VLC applications such as surveillance and monitoring, data collection from Internet-of-things devices, off-loading traffic data from base stations (BSs), and improving quality-of-service for the users. Recently, the concept of Twinkle was proposed in [START_REF] Deng | Twinkle: A flying lighting companion for urban safety[END_REF], where LEDs mounted on UAVs were used for illumination in urban areas. In addition, UAVs with LED-based lights have been reported for operations at night time [START_REF] Yang | Power efficient visible light communication with unmanned aerial vehicles[END_REF].

Use of UAVs as flying BSs to provide on-demand wireless communications offers a cost-effective and flexible solution in a number of applications [START_REF] Wang | Deep learning for optimal deployment of UAVs with visible light communications[END_REF]. Here, RF-based transmission may result in interference with the ground networks, and may also decrease the applicability of high-speed communications due to the limited energy [START_REF] Wang | Deep learning for optimal deployment of UAVs with visible light communications[END_REF]. VLC-based transmission is an interesting alternative, since it maximizes the resource utilization by using the existing LED-based lights. This combined feature of illumination and data communication is much favored in M. W. Eltokhey and M. A. Khalighi are with Aix-Marseille University, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France (e-mail: Mahmoud.Eltokhey@centrale-marseille.fr, Ali.Khalighi@fresnel.fr) Z. Ghassemlooy is with Optical Communications Research Group, Faculty of Engineering and Environment, Northumbria University, Newcastle upon Tyne, NE1 8ST, UK (e-mail: z.ghassemlooy@northumbria.ac.uk)

This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764461 (VisIoN), and was based upon work from European Union's Horizon 2020 COST Action CA19111 (NEWFOCUS).

search and rescue, delivery by drone at night, path illumination by UAVs in areas with no lighting, etc. [START_REF] Yang | Power efficient visible light communication with unmanned aerial vehicles[END_REF], as well as in large space indoor industrial scenarios. The use of VLC in UAV-based networks has received particular interest in the recent literature. For instance, optimizing UAV location and user association was considered in [START_REF] Yang | Power efficient visible light communication with unmanned aerial vehicles[END_REF] to minimize UAV power consumption, which was further extended in [START_REF] Wang | Deep learning for optimal deployment of UAVs with visible light communications[END_REF] using deep learning. Also, [START_REF] Pham | Sumrate maximization for UAV-assisted visible light communications using NOMA: Swarm intelligence meets machine learning[END_REF] considered network sum-rate maximization for non-orthogonal multiple access UAV-based VLC networks, by optimizing the power allocation and UAVs' placement.

In [START_REF] Yang | Power efficient visible light communication with unmanned aerial vehicles[END_REF]- [START_REF] Pham | Sumrate maximization for UAV-assisted visible light communications using NOMA: Swarm intelligence meets machine learning[END_REF], each user was associated with a single cell, which needs dealing with handover and cell-based interuser interference. To avoid handovers while mitigating interuser and inter-cell interferences [START_REF] Eltokhey | Multiple access techniques for large space indoor scenarios: A comparative study[END_REF], multiple-input multipleoutput techniques can be used [START_REF] Gesbert | From theory to practice: an overview of MIMO space-time coded wireless systems[END_REF] but the constraints on the receiver (Rx) array size to ensure a full-rank channel matrix makes multiple-input single-output (MISO) technique a more suitable solution in practice [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]. In multi-user (MU) MISO VLC systems with single PD-based Rxs [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF], [START_REF] Ma | Coordinated broadcasting for multiuser indoor visible light communication systems[END_REF], linear zero forcing (ZF) pre-coding is widely considered, as it offers simplicity and improved performance at high signal-to-noise ratios (SNRs) [START_REF] Eltokhey | Hybrid NOMA and ZF pre-coding transmission for multi-cell VLC networks[END_REF], [START_REF] Pham | Multi-cell VLC: Multi-user downlink capacity with coordinated precoding[END_REF]. However, the performance of linear ZF pre-coding is affected by the user locations, which impacts the correlation between the users' channel gains. In fact, a higher correlation in the network channel matrix results in a degraded ZF pre-coding performance, i.e., a lower network sum-rate. In this paper, we propose to decrease this correlation by optimizing the locations of UAVs using particle swarm optimization (PSO). Note that, PSO has been used in optical wireless communications, e.g., for optimizing the diffuse spot locations [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF] and resource allocation [START_REF] Demir | Unified resource allocation and mobility management technique using particle swarm optimization for VLC networks[END_REF]. Our main contributions include: (i) proposing the use of ZF pre-coding in UAV-based VLC networks, which to the best of our knowledge has not been addressed before; (ii) proposing to decrease the correlation in the network channel matrix by PSO-based UAV location optimization; and (iii) studying the effect of varying the optimization parameters on the network performance.

II. SYSTEM MODEL AND MATHEMATICAL FORMULATION

A. VLC System Model

Consider a UAV-based VLC system as in Fig. 1. Each UAV acts as an access point (AP), where the LED-based APs mounted on UAVs provide illumination and VLC downlink to the Rxs for data communication. A central control unit, which is mounted on one of the UAVs, or located in a separate central UAV (not used for direct communication with the Rxs), exchanges information with the APs and is, in turn, connected to a cellular network via an RF or a free-space optical link. Infrared uplinks are used for sending data and estimated channel-state information (CSI) to the APs; this latter being used in the central control unit for determining the ZF pre-coding matrix and for UAV location optimization.

Let us denote by h ij the channel coefficient between the i th AP and the j th Rx. We account only for the line-of-sight (LOS) path in h ij , given the heights of the APs and the dominance over the non-LOS paths. For complying with the constraints of ZF pre-coding, we consider that the number of Rxs (N r ) is less than or equal to the number of APs (N t ) [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]. The transmission by the considered LED luminaires can be modeled by Lambertian pattern of order m, where [START_REF] Ghassemlooy | Visible Light Communications: Theory and Applications[END_REF]:

h ij = ρ j S (m + 1) 2π l 2 ij cos m (φ ij ) cos(θ ij )A PDj . (1)
Here, φ ij and θ ij refer to the angles of emission and incident with respect to AP i and Rx j , respectively, and l ij denotes the related link distance. Also, ρ j and A PDj are the responsivity and the surface area of the PD, respectively, and S is the LED conversion efficiency. Note, for ij larger than the fieldof-view (FOV) of Rx j , we have h ij = 0.

B. ZF Pre-coding

In the considered MU-MISO system, to eliminate MU interference, ZF pre-coding is applied to the signals of the Rxs. At the i th AP, the transmit signal is [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]:

x i = Nr j=1 w ij d j + I DCi , (2)
where d j is the desired signal of the j th Rx, w ij is the precoding weight, and I DCi is the DC bias. At Rx j , after removing the DC offset, the received signal is given by:

r j = Nt i=1 x i h ij + z j = h T j w T j d j + k =j h T j w T k d k + z j , (3)
where h j = [h 1,j , ..., h Nt,j ] T is the vector of the channel gains associated with Rx j ; (•) T stands for transposition;

d = [d 1 , ..., d Nr ]
is the vector of the desired signals; w j = [w 1,j , ..., w Nt,j ] is the pre-coding vector associated with Rx j ; and z j is Gaussian noise with variance σ 2 n , representing the ambient light induced shot noise and the thermal noise. The 1 st and 2 nd terms on the right hand side in (3) represent the desired signal and the interference, respectively. Let r = [r 1 , ..., r Nr ] denote the vector of the received signals. We have:

r = H W d + z, (4)
where

H = [h T 1 , ..., h T Nr ] T , W = [w T 1 , ..., w T Nr ]
, and z = [z 1 , ..., z Nr ] denote the network channel matrix, the pre-coding weights matrix, and vector of noise component of all Rxs, respectively. For linear ZF pre-coding, W is given by [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]:

W = H T (H H T ) -1 diag(γ), ( 5
)
where γ denotes a diagonal matrix, composed of diagonal entries γ j > 0, such that γ j values can be regarded as the coefficients of parallel sub-channels determined based on a design criterion. For the sake of simplicity, we consider the criterion of maximizing the minimum achievable throughput. Defining µ = [µ 1 , ..., µ Nr ] = γ1 σ1 , ..., γ Nr σ Nr , this gives [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]:

µ = µ * 1; µ * = √ P e max(A 1) , (6)
where A = abs H T (HH T ) -1 diag(σ); σ is the Rxs' noise variance vector; abs(.) applies the absolute value operation for each element; 1 is a vector with all entries equal to one; and P e is the transmit electrical power per AP. For Rx j , the SNR and the maximum achievable throughput are given as, respectively:

SNR j = h T j w T j 2 σ 2 n = γ 2 j σ 2 n , (7)
R j = B 2 log 2 1 + SNR j ) (bps), (8)
where [START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF] assumes DC-biased optical-orthogonal frequencydivision multiplexing with a penalty of factor 2 due to Hermitian symmetry, and B denotes the system bandwidth. Note that, the achievable data rate depends on the UAVs' locations by affecting the path loss, as well as the correlation in H.

III. UAV LOCATION OPTIMIZATION

To optimize the UAV locations, we consider the wellestablished algorithm of PSO, which is a metaheuristic optimization technique that mimics the swarms movement [START_REF] Kennedy | Particle swarm optimization[END_REF]. Given the non-linear nature of the considered problem and the available degrees of freedom by change in the UAVs' locations, PSO is an efficient approach, which can converge to high quality solutions at low complexity. This is an important point, allowing real-time adaptivity to the mobility of the users, given the possible limitations on the computational resources.

A. Particle Swarm Optimization

In PSO, each point in the D-dimension search space represents a solution; particles travel between the points to find the optimal solution. For each particle's movements, both the speed and the directions are influenced by its best experience (i.e., its best solution so far, P best ), and by the best position found by the ensemble of particles (i.e., the global best position, G best ). The quality of the solution is evaluated based on a fitness function. Assuming the number of iterations N maxit and the number of particles N maxp , applied to an optimization problem with the number of variables D, for the particle q at iteration , the algorithm updates the vectors of particle's position p q = [p q1 , p q2 , ..., p qD ]; velocity v q = [v q1 , v q2 , ..., v qD ]; and personal best position P best,q = [p best,q1 , p best,q2 , ..., p best,qD ]. In addition, if a q th particle offers a superior fitness than that achieved by the position in the global best position vector G best = [g best,1 , g best,2 ..., g best,D ], G best is updated by the position of the q th particle. To update the velocity of particle q at iteration + 1 for the variable d, the following equation is used [START_REF] Kennedy | Particle swarm optimization[END_REF]:

v +1 qd = w v qd + c P rand P (p best,qd -p l qd ) + c G rand G (g best,d -p qd ), (9)
where the 1 st , 2 nd , and 3 rd terms correspond to the contributions from the old velocity of the particle, P best , and G best vectors, respectively; w refers to the inertia weight, which controls the impact of old velocity on the new value; c P and c G denote the weights for controlling P best and G best contributions to the new velocity; and rand P and rand G are random numbers taking values between 0 and 1, respectively. Defining ∆t as the time step (which is set to 1), the new position of the particle q at iteration + 1 for variable d is given by:

p +1 qd = p qd + v +1 qd ∆t. (10)
In our simulations we take 0.4 ≤ w < 0.9 as the appropriate interval allowing a good performance for the algorithm [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF]. Also, the new positions of the particles that exceed the solution space boundaries are clipped to the boundaries.

B. UAV Location Tuning Using PSO

The flow chart of the proposed PSO-based UAV location optimization is depicted in Fig. 2. Firstly, random positions and velocities for each particle are generated prior to evaluation of their associated UAV locations in the ZF pre-coding. The considered optimization problem targets decreasing the correlation in H for maximization of the sum-rate, Nr j=1 R j . To calculate the correlation in H, the 2-norm condition number, denoted cond(H), is considered, which is defined as the ratio of the largest singular value of H to the smallest one. To decrease this correlation, we consider minimizing cond(H).

The variables used for controlling cond(H) are the UAVs' locations that directly impact H, and are optimized using PSO. The sum-rate and cond(H) for the corresponding network configuration are used to determine a fitness function F for evaluating the solution quality. Here, we consider:

F = F 1 × Nr j=1 R j -F 2 × cond(H) -(F 3 × N 0 ), (11)
where N 0 is the number of users with no LOS link with any AP. The 1 st term in [START_REF] Pham | Multi-cell VLC: Multi-user downlink capacity with coordinated precoding[END_REF] promotes the solutions with higher sum-rates, the 2 nd term downgrades those resulting in a poor condition number of H, and the 3 rd term avoids solutions with users experiencing blocked LOS. Constants F 1 , F 2 , and F 3 control the contribution of the different terms. Note that, minimizing cond(H) may result in a less optimal sum-rate performance compared with when maximizing the sum-rate directly. However, the former approach results in a lower probability of obtaining an ill-conditioned H. Based on the value of F , both P best and G best are updated, followed by calculating the new position and velocity of the particle in the solution space. After evaluating each particle for all iterations, the G best vector will contain the optimized UAV locations.

When only discrete and non-uniformly separated UAV locations are considered (e.g., due to constraints on the localization accuracy), the possible UAV locations are converted to integers, and the problem is considered as an integer programming problem. Then, the indexes of UAV locations in the x and y coordinates form the solution space, to ensure equal spacing between successive solutions in each coordinate. The number of variables D is equal to the number of UAVs multiplied by the number of optimized variables per UAV.

To manage how particles search for the optimal solution in the solution space and benefit from their experiences, the weight c l P in ( 9) is decreased from 2.5 (at the start of optimization) to 0.5 (at the end) for a better exploration of the solution space at the beginning by maximizing the reliance on P best . On the other hand, for G best , c l G is increased from 0.5 (at the start of optimization) to 2.5 (at the end) to enhance the convergence to the optimal solution [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF].

The complexity of the proposed optimization mainly depends on calculating the particles' locations in PSO, evaluating the fitness function, and ZF pre-coding. For the considered PSO algorithm, the time and space complexities are O(N maxit × N maxp × D) and O(N maxp × D), respectively [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF].

IV. PERFORMANCE ANALYSIS

For the two cases of indoor industrial and outdoor scenarios, consider 4 LED-based APs, each mounted on a rotary-wing UAV, and 4 Rxs positioned at a height of 0.85 m above the floor level. The key system parameters are given in Table I. Considering a 10 MHz system bandwidth, the sum-rate R j corresponding to the 1 st term in (11) will be around tens of Mbps, whereas cond(H) can reach orders of 10 3 . Therefore, Fig. 3: Top view for the grid of possible UAV locations and the parameters considered for indoor industrial and outdoor scenarios.

we set F 1 , F 2 , and F 3 to 1, 10 5 , and 10 10 , respectively, to prioritize LOS blockage avoidance, while degrading the solutions with large condition numbers. In fact, although for the considered UAV heights it is unlikely that a Rx experiences a LOS blockage, the 3 rd term in [START_REF] Pham | Multi-cell VLC: Multi-user downlink capacity with coordinated precoding[END_REF] guarantees the applicability of the proposed algorithm to a wide range of scenarios. Figure 3 shows the top view for the considered communication environment, highlighting the grid of possible locations for the four UAVs, and the corresponding parameters for the considered indoor industrial and outdoor scenarios.

A. Effect of Parameter Optimization

To find the best combinations of the numbers of particles and iterations for a given number of evaluations, we compare in Fig. 4 the average sum-rate for different combinations of N maxp and N maxit , over 1000 random user positions, for both indoor and outdoor scenarios. As shown, the average sum-rate improves with increasing N maxp × N maxit . For the same evaluation number, generally, a higher average sum-rate is achieved for N maxp > N maxit , which is because of a better exploration of the solution space prior to converging to the optimal solution. Note that, lower sum-rates are achieved for the outdoor scenario since UAVs are considered with higher altitudes with the same transmit power as for the indoor case. Also, note that the required number of evaluations is much smaller than when evaluating all possible solutions (4 UAVs, with 25 possible locations each, see Fig. 3). Based on these results, we set N maxp = 20 and N maxit = 10 in the sequel.

B. Case Study of UAVs With Optimized Locations

To see how the sum-rate improves with the proposed optimization through iterations, we consider in Fig. 5 an indoor industrial scenario with (x, y) of the Rxs' positions of (0.26, 4.63), (2.93, 3.55), (6.17, 3.18), and (7.54, 9.87). The evolution of the fitness function F , from [START_REF] Pham | Multi-cell VLC: Multi-user downlink capacity with coordinated precoding[END_REF], is illustrated in the dashed blue boxes, together with the UAVs' location indexes (with respect to the grid), sum-rate, and cond(H), over every change in F (i.e., the best performance achieved). The dashed red box shows the case UAVs are positioned at the center of their corresponding grid space. Also, at the right side, we have shown the optimal positions of the UAVs with respect to the Rxs positions. We notice that PSO converges at iteration 5 to the optimal solution. Meanwhile, no improvement is achieved from iteration 2 to 3, i.e., the new positions found in iteration 3 did not result in an improvement of G best . In fact, in the first iterations, PSO focuses on the solution space exploration, rather than convergence to the optimal solution. It is worth mentioning that the obtained solution via PSO after convergence (i.e., a sum rate of 0.973 × 10 8 bps) is quite close to the global optimal solution of 1.007 × 10 8 bps, obtained via parameter sweep (with considerably higher complexity).

C. Optimizing UAV Locations in Case of Mobility

In a dynamic scenario where the Rxs move around, the UAVs' locations need to be updated accordingly to ensure the best network performance. Figure 6 shows a comparison of the average sum-rate over different Rx mobilities, for the indoor industrial and outdoor scenarios, for different optimization rates used. For each mobility case, 1000 users' positions are generated according to the random way-point mobility model, where the time interval between each two successive positions is set to 0.5 sec. The considered maximum user speeds of 0.5 and 1.4 m/s correspond to walking and fast walking users, whereas the cases of 2 and 5 m/s correspond to equipment movements in current and future industrial scenarios [START_REF] Sauer | Delay tolerant networks in industrial applications[END_REF], respectively. The considered time interval between optimizations are 0.5, 1, 2, and 5 sec. We have also considered the nonoptimized cases where UAVs are located at their grid space in order to vary the inter-UAV distances from the minimum to the maximum. For the case of indoor scenario, as UAVs get closer, the similarity in the users' channel gains increases, resulting in an increased correlation in H, and subsequently a decreased sum-rate. For outdoor scenarios, the advantages of maximizing the inter-UAV distances are reduced due to the increased path loss (due to increased LOS path length between the UAVs and the Rxs). Given the highlighted trade-off between the decrease in the channel correlation and the increase in the channel gain, the maximum distance between UAVs should be constrained in the system design by parameters such as the users' data rate requirements. On the other hand, results show an improved performance for optimized UAVs' locations for all cases. As expected, decreasing the optimization rate results in a performance degradation, which is more considerable with increased mobility. Indeed, larger variations in the users' positions over the given time window results in a faster deviation from the optimal UAVs configuration. Interestingly, even for the highest mobility case, updating optimized UAV positions at the lowest rate (i.e., every 5 sec) offers an improved performance, thus demonstrating the robustness of the proposed optimization.

Fig. 1 :

 1 Fig. 1: Illustration of the considered system with the flow of information between a cellular network, a central control unit, UAV APs, and Rxs.

Fig. 2 :

 2 Fig.2: Flow chart of the proposed algorithm for location optimization of three UAVs, where and q are the indexes of the number of iterations and particles, respectively. Dashed boxes on the left and right illustrate how UAV locations are converted to integers (Int.), and how two particles move in the solution space to find the optimal solution, respectively.

Fig. 4 :

 4 Fig. 4: Effect of Nmax p and Nmax it on the average sum-rate of 4 Rxs with 1000 random positions, in indoor industrial and outdoor scenarios.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Evolution of fitness value over iterations. Blue and red boxes correspond to optimized and non-optimized cases, respectively.

Initialize random position and velocity vectors Iteration tℓ Convert integers in solution space to UAV locations Update position and velocity Use UAV locations in ZF pre-coding scenario Particle pq Optimized UAV locations are equal to G best vector

	𝐢 𝟓						𝑮 𝐛𝐞𝐬𝐭 𝑷 𝐛𝐞𝐬𝐭
					No	Yes	Particle position
	y-axis					ℓ > 𝑵 𝐦𝐚𝐱𝒊𝒕 ?	Original velocity
	𝐢 𝟏	𝐢 𝟏	x-axis	𝐢 𝟓	No	Yes 𝒒 > 𝑵 𝐦𝐚𝐱𝒑 ? ℓ = ℓ + 1	Resultant velocity Velocity towards 𝑮 𝐛𝐞𝐬𝐭 Velocity towards 𝑷 𝐛𝐞𝐬𝐭
	𝐔𝐀𝐕𝐢 location
						q = q + 1	Location
							of UAV3
	Int. of							1	2
	for 𝐔𝐀𝐕𝐢 x-coordinate Int. of y-coordinate for 𝐔𝐀𝐕𝐢	Calculate sum-rate and condition no.	Evaluate Fitness Update 𝑮 𝐛𝐞𝐬𝐭 and 𝑷 𝐛𝐞𝐬𝐭	Location of UAV1	1	of UAV2 Location 2

TABLE I :

 I Simulation parameters

		Parameter				Value
	LED conversion efficiency S		0.44 W/A [9]
	Transmit power per LED			15.84 W
	LED Lambertian order m			1
	Noise power spectral density (indoor/outdoor night)	10 -21 A 2 /Hz [10]
		Link bandwidth			10 MHz [10]
	PD responsivity ρ j			0.4 A/W [10]
		PD area				1 cm 2 [10]
		Optical Rx FOV			62 deg. [8]
		Rx height				0.85 m [8]
				X max
		13		XU max
			YU max
	y axis	5 9	Δxu	grid UAV	UAV	Y max
			Δyu
		1	5	9	13
			x axis
		Parameter	Indoor (m)	Outdoor (m)
		X max		10	20
		Y max		10	20
		XU max		2	4
		YU max		2	4
		Δxu		0.5	1
		Δyu		0.5	1
		UAV height	8	12

V. CONCLUSIONS

We proposed optimization of the UAV locations using PSO, for MU-MISO ZF pre-coded UAV-based VLC networks and demonstrated the achieved performance improvement. We showed the importance of carefully choosing the fitness function as well as the adequate number of particles for a sufficient exploration of the solution space. Generally, as the number of variables (i.e., UAVs) increases, a larger number of evaluations should be considered to avoid being trapped in local solutions. We also considered different mobility scenarios and showed the robustness of the proposed solution for different optimization rates due to limited computational resources and/or high user mobility. The proposed optimization could be extended to carry out the pre-coder design while accounting for CSI uncertainty due to outdated or noisy CSI, as proposed in [9], which could be addressed in a future work.


Mahmoud Wafik Eltokhey
email: mahmoud.eltokhey@centrale-marseille.fr

Ali Khalighi
email: ali.khalighi@fresnel.fr

Zabih Ghassemlooy
email: z.ghassemlooy@northumbria.ac.uk

Power Allocation Optimization in NOMA-Based Multi-Cell VLC Networks

Keywords: Visible-light communications, Multiple-access techniques, NOMA, Particle swarm optimization

The use of visible-light communication (VLC) to complement radio-frequency systems in satisfying the everincreasing user data demands has received particular interest since about two decades. One important issue in VLC networks is to address the multiple-access (MA) requirement. Non-orthogonal MA (NOMA) is a promising technique for multi-cell VLC networks, in particular, due to providing high spectral efficiency. However, its performance is highly impacted by the considered power allocation (PA). In this paper, we propose optimization of PA in NOMA-based VLC networks using particle swarm optimization. The presented results show that the optimized PA offers a higher performance in terms of the network achievable throughput and fairness, compared with non-optimized PA, while offering flexibility in controlling the desired performance criteria.

I. INTRODUCTION

Visible-light communication (VLC) has been receiving an increasing interest for supporting indoor and outdoor wireless access, due to several features including: the presence of unregulated spectrum, immunity against the radio-frequency (RF) interference, inherent physical layer security, and the reliance on using light-emitting diode (LED)-based lighting infrastructures [START_REF] Ghassemlooy | Visible Light Communications: Theory and Applications[END_REF]. To support lighting and data communication in relatively large space scenarios, multiple LED luminaires are used, each one serving as an access point (AP), providing service to the users within its coverage area. In such a multicell VLC network, according to the locations of the users with respect to the APs, users can be classified as either of cell-center users (CCUs) or cell-edge users (CEUs). The former refers to the users within the coverage area of only one AP, while the latter are the users located within the coverage areas of more than one AP. In such networks, the users performances are affected by the inter-cell interference (ICI), resulting from received signals from neighbouring cells, and inter-user interference (IUI) due to the received signals of other users in the same cell.

One efficient multiple-access (MA) technique to minimize IUI and ICI effects is the power-domain NOMA (simply referred to as NOMA here), which exploits the differences in the channel gain between users to allow multiplexing their signals [START_REF] Maraqa | A survey of rate-optimal power domain NOMA with enabling technologies of future wireless networks[END_REF]. At the receiver (Rx) side, successive interference cancellation (SIC) is done to separate the signals of different users [START_REF] Kizilirmak | Non-orthogonal multiple access (NOMA) for indoor visible light communications[END_REF]. Several works have compared the performances of orthogonal frequency-division MA (OFDMA) and NOMA in single-cell and multi-cell scenarios, where the advantage of the latter scheme was shown in terms of network throughput, especially for increased number of users [START_REF] Kizilirmak | Non-orthogonal multiple access (NOMA) for indoor visible light communications[END_REF]- [START_REF] Eltokhey | Multiple access techniques for large space indoor scenarios: A comparative study[END_REF].

Several works proposed solutions for improving the sumrate in single-cell NOMA-VLC networks. In [START_REF] Yang | Fair non-orthogonal multiple access for visible light communication downlinks[END_REF]- [START_REF] Li | Optimal power allocation scheme based on multi-factor control in indoor NOMA-VLC systems[END_REF], power allocation (PA) schemes were proposed, which outperformed the gain ratio PA. Cuckoo Search (CS) optimization was used in [START_REF] Tahira | Optimization of non-orthogonal multiple access based visible light communication systems[END_REF] to optimize the VLC system parameters, including PA, with the objective of maximizing both the received power and the achievable rate. In [START_REF] Feng | Joint power allocation for the multi-user NOMA-downlink in a power-line-fed VLC network[END_REF], a joint power-line communication (PLC)-VLC PA strategy was proposed to maximize the network sum-rate. Optimal user pairing followed by PA using majorization minimization was considered in [START_REF] Eroglu | Slow beam steering and NOMA for indoor multi-user visible light communications[END_REF] to maximize the sum-rate for VLC with steerable beams.

In contrast, fewer works have investigated IUI and ICI management in multi-cell NOMA-VLC networks. In [START_REF] Feng | Multiple access design for ultra-dense VLC networks: Orthogonal vs non-orthogonal[END_REF], overlapped clustering based on hybrid NOMA-orthogonal multiple access (OMA) were considered in layered asymmetricallyclipped optical orthogonal-frequency division multiplexing (LACO-OFDM) VLC networks, for minimizing interference. User assignment was considered in [START_REF] Zhang | User grouping and power allocation for NOMA visible light communication multi-cell networks[END_REF] for ICI mitigation, where PA was optimized using the gradient projection algorithm for sum-rate maximization. In [START_REF] Eltokhey | Hybrid NOMA and ZF pre-coding transmission for multi-cell VLC networks[END_REF], ICI mitigation was carried out by using zero-forcing pre-coding for broadcasting the user signals as the highest decoding order NOMA user, while [START_REF] Eltokhey | Dimming-aware interference mitigation for NOMA-based multi-cell VLC networks[END_REF] coordinated the scheduling of NOMA transmission in each cell to mitigate ICI, while satisfying the dimming requirements in the network.

In this paper, we propose the optimization of PA in multicell NOMA-VLC networks using particle swarm optimization (PSO), which has been used in solving different problems in optical communications [START_REF] Demir | Unified resource allocation and mobility management technique using particle swarm optimization for VLC networks[END_REF]- [START_REF] Eltokhey | Optimization of receivers field-of-views in multi-user VLC networks: A bio-inspired approach[END_REF]. By tuning the PA in each of the cells, both IUI and ICI in the network can be controlled, while by tuning the optimization parameters, the desired performance could be improved, considering the trade-off between network throughput maximization and the homogeneity of user performance (i.e., fairness). The results show improved performance compared with non-optimized PA in terms of both network throughput and fairness. The remainder of the paper is organized as follows: Section II presents the system model, Section III outlines a brief analysis of NOMA signalling. Section III presents the proposed PSO-based PA optimization, and Section V discusses the results of comparing the optimized and non-optimized PA. Finally, Section VI concludes the paper.

II. SYSTEM MODEL

We consider an intensity modulation/direct detection multicell VLC network, that has N t APs (LED luminaires) handling a total number of N r users, such that each AP i handles N i users, and the channel gain for the link between AP i and the user Rx ij is denoted by h ij . Assuming Lambertian pattern for the LED transmission of order O and that the line-of-sight component dominates the channel gain, we have [START_REF] Long | Investigating channel frequency selectivity in indoor visible light communication systems[END_REF]:

h ij = S (O + 1) R A 2π ℓ 2 ij cos O (ϕ ij ) cos(θ ij ), (1)
where S refers to the LED conversion efficiency, ϕ ij is the transmission angle at the AP i , θ ij is the angle of incidence at the Rx ij , R is the responsivity of the photo-detector (PD), and ℓ ij is the path length, see Fig. 1. Also, denoting the active area of the PD, the Rx field-of-view (FOV), and the optical concentrator's refractive index by A PD , θ c , and q, respectively, the Rx collection area A is given by:

A = q 2 sin 2 (θ c ) A PD . (2)
III. NOMA SIGNALING Considering N t APs, each one handling N i users, i = 1, ..., N t , the users' signals are multiplexed by the AP in power domain using superposition coding. For this, users' signals are assigned different power levels, where the users with the lowest channel gains are assigned the highest power levels, and have the lowest decoding order. At the users side, SIC detection is carried out, where users with higher decoding order decode the signals of the users with lower decoding order before accessing their signals. Assuming availability of CSI at all Rxs, the signal received by Rx ij handled by AP i is given by:

r ij = a ij P e h ij d ij + j-1 k=1 a ik P e h ij d ik + Ni k=j+1 a ik P e h ij d ik + z j , (3)
where P e refers to the total electrical transmit power of the AP i , a ij represents the PA weight, d ij denotes the desired data of Rx ij , and z j refers to Gaussian noise with variance σ 2 n . In (3), the first, second, and third terms represent the desired signal, the interference components cancelled by SIC, and residual interference following SIC, respectively. By assuming static PA for sake of simplicity, we have a 2 ij = α a 2 i j-1 , such that 0 < α < 1 refers to the PA coefficient, which decides the ratio between the power allocated to a user and the preceding user in the decoding order. Note that, Ni j=1 a 2 ij = 1 for ensuring normalized electrical transmit power. The signalto-interference-plus-noise ratio (SINR) for the Rx ij is then calculated as:

SINR NOMA,Rxij = h 2 ij P e a 2 ij I ICI + h 2 ij P e k>j a 2 ik + σ 2 n , (4)
where I ICI refers to the ICI component. The maximum achievable throughput of Rx ij is then given by:

R ij = B 2 log 2 1 + SINR NOMA,Uij (bps), (5)
where B denotes the system bandwidth, and the loss factor of 2 is due to the Hermitian symmetry constraint, assuming DCO-OFDM signalling.

To measure the homogeneity of users' throughput performance, we consider Jain's fairness index, defined as [START_REF] Eltokhey | Dimming-aware interference mitigation for NOMA-based multi-cell VLC networks[END_REF], [START_REF] Jain | A quantitative measure of fairness and discrimination[END_REF]:

FI = 1 N r   Nt i=1 Ni j=1 R ij   2 Nt i=1 Ni j=1 R 2 ij . (6)
FI = 1 in case that all R ij are equal; it decreases with increase in the difference between the users achievable throughputs.

IV. PA OPTIMIZATION USING PARTICLE SWARM OPTIMIZATION

For optimizing the PA coefficients to improve the network performance, we consider the use of PSO algorithm [START_REF] Kennedy | Particle swarm optimization[END_REF]. In PSO, particles move in the solution space to find their optimal solution, in a way similar to the behavior of swarms of bees searching for locations with highest densities of flowers [START_REF] Kennedy | Particle swarm optimization[END_REF], [START_REF] Shi | A modified particle swarm optimizer[END_REF].

The movements of each of the particles in PSO are influenced by the best position reached by the particle itself so far, which is referred to as personal best position P best, and by the global best position Gbest reached by all particles. For evaluating the quality of the solution according to the desired performance criteria, a fitness function is considered. For the number of iterations N it and of particles N p used in D dimensions optimization problem (i.e., D optimized variables), at iteration l, the PSO algorithm updates, for the particle m, the vectors of position p m = [p m1 , p m2 , ... , p mD ]; velocity v m = [v m1 , v m2 , ... , v mD ]; and personal best position P best m = [pbest m1 , pbest m2 , ... ,pbest mD ]. In case that the m th particle performance is better than that of the global best position vector Gbestv = [gbest 1 , gbest 2 , ... , gbest D ], the position of the particle is used as the Gbestv. The velocity of the m th particle at iteration l for the variable d is then given by [START_REF] Kennedy | Particle swarm optimization[END_REF]:

v l+1 md = w l v l md + c l P rand P (pbest md -p l md ) + c l G rand G (gbest d -p l md ), (7)
where the first, second, and third terms refer to the contributions from the components of the particle's old velocity, and the vectors of P best and Gbest, respectively. w l denotes the inertia weight that controls the contribution of the old value of the velocity on calculation of its new value, and rand P and rand G refer to random numbers bounded between 0 and 1. c ℓ P and c ℓ G correspond to the weights that control the contributions of P best and Gbest vectors on determining the new velocity. Considering the time step ∆t, the new position of the m th particle at iteration l + 1 for the variable d is:

p l+1 md = p l md + v l+1 md ∆t. (8)
For optimizing PA coefficients, in order to highlight the role of the fitness function in controlling the network performance, three different fitness functions are considered. The first fitness function focuses only on maximizing the sum-rate:

F SR = Nt i=1 Ni j=1 R ij , (9)
which promotes solutions resulting in sum-rate improvement.

The second fitness function maximizes both sum-rate and FI:

F SR+F I = Nt i=1 Ni j=1 R ij + W 1 FI, ( 10
)
where the first and second terms promote solutions improving the sum-rate and FI, respectively. The weight W 1 ensures comparable contributions of the two terms. Assuming the sumrate on the order of 10 Mbps and given the maximum FI of 1, W 1 is set to 10 8 here. The third considered fitness function focuses only on maximizing the FI:

F F I = W 2 FI, ( 11
)
where W 2 is set to 10 9 . Figure 2 shows the flow chart for the proposed optimization algorithm. It highlights generating random positions, evaluation of the positions (PA coefficients), and the update of the new positions to be explored, for finding the optimal solution (i.e., the final global best position).

V. PERFORMANCE STUDY

We compare here the performances of the non-optimized PA and the proposed optimized PA over different scenarios.

A. Main Assumptions and Considered Scenarios

We consider three different 4-cell VLC network scenarios, with different combinations of CCUs and CEUs. In Scenario 1, N r = 8 and the number of CEUs is set to 4; in Scenario 2, N r = 12 including 4 CEUs; and in Scenario 3, N r = 16 including 8 CEUs. The APs are positioned at the cell centers with the transmit optical power P o = 1.584 W [START_REF] Ma | Coordinated broadcasting for multiuser indoor visible light communication systems[END_REF]. Also, the heights of the APs and the Rxs from the floor level are set to 2.5 and 0.85 m, respectively. We consider N p = N it = 10 as appropriate parameters for PSO. For each scenario, 200 random user positions are generated to determine the average sum-rate and fairness performance. 1 Table I summarizes the simulation parameters. Note that, the PA coefficients in any cell are bounded between 0.1 and 0.5. Also, for the cases of non-optimized PA, we consider PA coefficients of 0.1, 0.2, 0.3, 0.4, and 0.5. The PSO parameters w ℓ , c l P , and c l G are varied according to [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF], [START_REF] Shi | A modified particle swarm optimizer[END_REF], [START_REF] Ratnaweera | Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients[END_REF] to ensure a quick conversion of the algorithm.

B. Numerical Results

Figure 3 contrasts the average sum-rate and fairness performances over the 200 random user positions in Scenarios 1, 2, and 3 for the cases of non-optimized and the optimized PA.

Consider first the case of non-optimized PA coefficients, where the results show that increasing the PA coefficient results in increased network sum-rate and decreased FI. This results from the more opportunistic approach adopted by the network. Indeed, increasing the PA coefficients results in allocating higher power levels for the users with higher decoding orders (having higher channel gains). This results in maximized network sum-rate but in decreasing the power allocated to the users with lower decoding order, thus penalizing them. Consequently, the margin between the achievable throughputs of users is increased, thus the fairness performance. We notice an improvement in the sum-rate performance when moving from Scenario 1 to Scenario 2, due to handling of more users. Also, we notice a degraded sum-rate performance by moving from Scenario 2 to Scenario 3, due to the increase in the IUI experienced by the users (see ( 4)).

Generally, the results show the merits of adopting optimized versus non-optimized PA coefficients. In particular, by adopting F SR and F SR+F I fitness functions, an improved sumrate performance is achieved for all considered values of nonoptimized PA. On the other hand, adopting F F I results in an improved FI performance, compared with the cases with non-optimized PA. Note that, by moving from the fitness functions F SR to F SR+F I and then F F I , the sum-rate performance degrades, while the fairness performance improves. This was expected, since F SR considers only the sum-rate, thus resulting in the best sum-rate and the worst fairness performances. On the other hand, F F I promotes only the fairness performance, thus resulting in the best FI performance and the lowest sum-rate. The F SR+F I fitness function offers a tradeoff between sum-rate and fairness, compared to the two other ones. The results show that, by tuning the fitness function, the desired performance criterion could be improved. This could be considered as a degree of freedom in adapting NOMA scheme to the changing network conditions or requirements. Note that, the space and time complexities of the considered PSO algorithm are given by O(N p ×D) and O(N it ×N p ×D), respectively [START_REF] Eltokhey | Optimization of intensities and locations of diffuse spots in indoor optical wireless communications[END_REF].

VI. CONCLUSIONS

We proposed PSO-based optimization of PA coefficients in NOMA-based multi-cell VLC networks using different fitness functions, and compared the achieved performance with the non-optimized PA NOMA case. Our results showed improved performance in terms of sum-rate and fairness. In particular, tuning the fitness function provides the flexibility to adapt to the changing network requirements or criteria. The achieved performance improvement by optimized PA comes at the cost of increased network complexity, which remains rather acceptable, given the relatively small number of evaluations needed for carrying out the optimization. Future work will consider applying the proposed optimization to other PA techniques (such as gain ratio PA), and including other system parameters in the optimization (e.g., user association for CEUs).

Fig. 1 :

 1 Fig. 1: Illustration of the link parameters.

Fig. 2 :

 2 Fig. 2: Flow chart of the proposed PA optimization algorithm. Parameters l and m denote the indexes of numbers of iterations and particles, respectively.

Fig. 3 :

 3 Fig. 3: Comparison of (a) sum-rate and (b) fairness index for the different considered optimization fitness functions and different values for non-optimized PA coefficients.

TABLE I :

 I Simulation parameters[START_REF] Ma | Coordinated broadcasting for multiuser indoor visible light communication systems[END_REF],[START_REF] Yu | Multi-user MISO broadcasting for indoor visible light communication[END_REF]

	Parameter	Value
	Room dimension	(7 × 7 × 3 m 3 )
	LED Lamertian order O	1
	Number of LED chips per LED luminaire	36
	Signal current per LED chip	100 mA
	LED conversion efficiency S	0.44 W/A
	PD area	1 cm 2
	Refractive index of the optical concentrator	1.5
	Rx's FOV	62 deg.
	PD responsivity R	0.4 A/W
	System bandwidth B	10 MHz
	Equivalent Rx noise power spectral density	10 -21 A 2 /Hz

Note that in the simulations we exclude those scenarios resulting in CEUs located in coverage areas of more than two APs for simplicity. However, the proposed optimization could be applied to such cases in a future work.

ACKNOWLEDGMENT

This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764461 (Vi-sIoN), and was also based upon work from European Union's Horizon 2020 COST Action CA19111 (NEWFOCUS).