# On the Combined Inverse-Square Effect of Multiple Point Sources in
Multidimensional Space

Keaton Coletti
College of Engineering
University of Georgia - Athens
Athens, GA 30602
e-mail<EMAIL_ADDRESS>Pawel Kalczynski
College of Business and Economics
California State University - Fullerton
Fullerton, CA 92834.
e-mail<EMAIL_ADDRESS>Zvi Drezner
College of Business and Economics
California State University - Fullerton
Fullerton, CA 92834.
e-mail<EMAIL_ADDRESS>
###### Abstract

The inverse-square law states that the effect a source has on its surroundings
is inversely proportional to the square of the Euclidean distance from that
source. Its applicability spans multiple fields including physics,
engineering, and computer science. We study the combined effect of multiple
point sources surrounding a closed region in multidimensional space. We
determine that the maximum effect in D dimensions is on the region’s boundary
if $D\leq 4$, and the minimum is on the boundary if $D\geq 4$.

Keywords: inverse-square law, dimensionality reduction, global optimization

##  1 Introduction

This paper investigates Euclidean distances in multidimensional space. Point
sources emit some quantity whose effect is inversely proportional to the
squared distance from the source, and the total effect at a point in space is
the sum of scalar contributions from each source. This scalar sum applies to
radiation, sound, and light. The cooperative cover location problem is
proposed in Berman et al., (2010). There are sources of light, cell phone
signals, noise, heat, etc. that decline by the square of the distance. Each
point in a region is affected by the sum of the signals. For example, light
posts in a parking lot are in three-dimensional space, and we are interested
in the total light in the two-dimensional parking lot. Which point in the
parking lot receives the minimum amount of light? This minimum value
guarantees that all the points in the parking lot receive at least this
minimum light.

Locating global extrema of such quantities is often of interest (Hansen et
al.,, 1981; Drezner et al.,, 2020). The results presented in this work reduce
optimization dimensionality by 1, thereby improving computation time and
reducing the risk of erroneously identifying local extrema.

As an example, consider one or more closed shapes in 3 dimensions surrounded
by point sources. The sources emit sound, heat, light, or other types of
energy in every direction. The emission from these sources may have different
levels of intensity. If the optimization objective is to arrange the sources
and/or the shapes in a manner that reduces the maximum total (cumulative)
effect of the emissions on the shapes (keeping all sources outside of the
shapes), we show that it is sufficient to consider solely the effect on the
boundaries of these shapes because, in 3 dimensional space, the maximum total
effect must be on the boundary.

##  2 Problem Statement

Let $D\in\\{1,2,3,...\\}$ be the dimension of a Euclidean space
$\mathbb{R}^{D}$. Consider $n$ point sources
$S_{i}=(s_{1i},s_{2i},...,s_{Di})\in\mathbb{R}^{D}$, $i=1,2,...,n$. Define
$w_{j}>0$ which represent the strength of the quantity emitted by source $i$.

Furthermore, let $\overline{\Omega}\in\mathbb{R}^{D}$, be a region affected by
the point sources, such that
$\overline{\Omega}\,\cap\,\\{S_{1},S_{2},...,S_{n}\\}=\varnothing$. The
combined effect of point sources on $X=(x_{1},x_{2},...,x_{D})$ is given by

$J(X)=\sum\limits_{i=1}^{n}\frac{w_{i}}{\sum\limits_{k=1}^{D}(x_{k}-s_{ki})^{2}}.$
(1)

##  3 Maximum and Minimum Values

###### Theorem 1

: Let $\overline{\Omega}\in\mathbb{R}^{D}\setminus\\{S_{i}\\}$ contain the
open set $\Omega$ and its boundary $\partial\Omega$. If $J(X)$ is non-constant
on $\Omega$,

  1. 1.

If $D\leq 4$, then $J(X)$ obtains maximum in $\overline{\Omega}$ only on
$\partial\Omega$.

  2. 2.

If $D\geq 4$, then $J(X)$ obtains minimum in $\overline{\Omega}$ only on
$\partial\Omega$.

Proof:  $J(X)$ is defined on $\overline{\Omega}$ and hence attains a maximum
and minimum in $\overline{\Omega}$. $J(X)$ is a continuous function and has
continuous derivatives of every order in $\overline{\Omega}$. To simplify
notation, define $x_{ki}=(x_{k}-s_{ki})$ for $k=1,2,...,D,i=1,2,...,n$ (note
that $dx_{ki}=dx_{k}$). Then Eq. (1) can be written as

$J(X)=\sum\limits_{i=1}^{n}\frac{w_{i}}{\sum\limits_{k=1}^{D}x_{ki}^{2}},$

and its unmixed first and second-partial derivatives are

$\frac{dJ(X)}{dx_{k}}=\sum\limits_{i=1}^{n}\frac{-2x_{ki}w_{i}}{\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{2}}$

$\frac{d^{2}J(X)}{dx_{k}^{2}}=\sum\limits_{i=1}^{n}\frac{-2w_{i}\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{2}+8x_{ki}^{2}w_{i}\sum\limits_{k=1}^{D}x_{ki}^{2}}{\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{4}}.$

Computing the Laplacian yields

$\Delta
J(X)=\sum\limits_{k=1}^{D}\frac{d^{2}J(X)}{dx_{k}^{2}}=\sum\limits_{i=1}^{n}\frac{-2Dw_{i}\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{2}+8w_{i}\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{2}}{\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{4}}=\sum\limits_{i=1}^{n}\frac{(8-2D)w_{i}}{\left[\sum\limits_{k=1}^{D}x_{ki}^{2}\right]^{2}}.$
(2)

We have two cases:

  1. 1.

$\Delta J(X)\geq 0$ for $D\leq 4$, so $J$ is sub-harmonic For all
$X\in\overline{\Omega}$, $\frac{d^{2}J(X)}{dx_{k}^{2}}>0$ for some $k$, so $J$
has no local maxima or plateaus in $\overline{\Omega}$.

  2. 2.

$\Delta J(X)\leq 0$ for $D\geq 4$, so $J$ is super-harmonic For all
$X\in\overline{\Omega}$, $\frac{d^{2}J(X)}{dx_{k}^{2}}<0$ for some $k$, so $J$
has no local minima or plateaus in $\overline{\Omega}$.

The maximum principle for sub-harmonic functions (case 1) and minimum
principle for super-harmonic functions (case 2) is applied by Stein and
Shakarchi, (2010); Axler et al., (2013); Evans, (2022); Ovall, (2016). $\Box$

###### Corollary 1

: If $D=4$, $J$ is harmonic on its domain and attains both maximum and minimum
values in $\overline{\Omega}$ on the boundary.

Note that $\overline{\Omega}$ need not be a connected domain. No statement has
yet been made about the maximum of $J(X)$ for $D>4$ or the minimum of $J(X)$
for $D<4$. Theorem 1 does not apply for other values of $D$.

###  3.1 When Maximum and Minimum Values Can Be in the Interior

We construct a counter example for which in $D\geq 5$ the maximum can be in
the interior of $\overline{\Omega}$, and in $D\leq 3$ the minimum can be in
the interior. Consider a set of $2D$ source points in D-dimensional space, as
depicted in Table 1 for $D=5$, with $w_{i}=1$, $i=1,2,..,2D$. All the points
are on the surface of a D-dimensional sphere centered at the origin (0, 0,
$\ldots$, 0) with a radius of $\frac{1}{2}$. We first show that the center is
a local maximum for $D\geq 5$, and a local minimum for $D\leq 3$. We then find
a radius $0<r_{max}<\frac{1}{2}$ of a D-dimensional sphere that defines the
region $\overline{\Omega}$.

Table 1: Counter example for $D=5$ | $k=1$ | $k=2$ | $k=3$ | $k=4$ | $k=5$
---|---|---|---|---|---
$i=1$ | $\frac{1}{2}$ | 0 | 0 | 0 | 0
$i=2$ | 0 | $\frac{1}{2}$ | 0 | 0 | 0
$i=3$ | 0 | 0 | $\frac{1}{2}$ | 0 | 0
$i=4$ | 0 | 0 | 0 | $\frac{1}{2}$ | 0
$i=5$ | 0 | 0 | 0 | 0 | $\frac{1}{2}$
$i=6$ | $-\frac{1}{2}$ | 0 | 0 | 0 | 0
$i=7$ | 0 | $-\frac{1}{2}$ | 0 | 0 | 0
$i=8$ | 0 | 0 | $-\frac{1}{2}$ | 0 | 0
$i=9$ | 0 | 0 | 0 | $-\frac{1}{2}$ | 0
$i=10$ | 0 | 0 | 0 | 0 | $-\frac{1}{2}$

Consider a location $\Delta X=\\{\Delta x_{k}\\}$ for $k=1,\ldots,D$ on a
D-dimensional sphere of radius $r$ centered at the origin so that
$\sum\limits_{k=1}^{D}\Delta x_{k}^{2}=r^{2}$. For source point $i=1$,
$\sum\limits_{k=1}^{D}(\Delta x_{k}-x_{i})^{2}=(\frac{1}{2}-\Delta
x_{1})^{2}+\sum\limits_{k=2}^{D}\Delta x_{k}^{2}=\frac{1}{4}-\Delta
x_{1}+r^{2}$. Similarly, for $i=D+1$, $\sum\limits_{k=1}^{D}(\Delta
x_{k}-x_{i})^{2}=(-\frac{1}{2}-\Delta x_{1})^{2}+\sum\limits_{k=2}^{D}\Delta
x_{k}^{2}=\frac{1}{4}+\Delta x_{1}+r^{2}$, and similar expressions for the
other points. Therefore,

$\displaystyle J(\Delta X)$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{D}\left[\frac{1}{\frac{1}{4}-\Delta
x_{i}+r^{2}}+\frac{1}{\frac{1}{4}+\Delta
x_{i}+r^{2}}\right]=\sum\limits_{i=1}^{D}\frac{\frac{1}{2}+2r^{2}}{(\frac{1}{4}+\Delta
x_{i}+r^{2})(\frac{1}{4}-\Delta x_{i}+r^{2})}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{D}\frac{\frac{1}{2}+2r^{2}}{(\frac{1}{4}+r^{2})^{2}-\Delta
x_{i}^{2}}~{}.$

Since $\Delta x_{i}\leq r$, we can get an upper bound for $J(\Delta X)$ by
replacing $\Delta x_{i}$ by $r$, and a lower bound by replacing $\Delta x_{i}$
by 0. However, the gap in this approximation is too large and has an error of
the order $r^{2}$ which fails when used to prove local optima. A better
approximation of order $r^{4}$ is obtained by multiplying both numerator and
denominator by $(\frac{1}{4}+r^{2})^{2}+\Delta x_{i}^{2}$ getting:

$J(\Delta
X)=2\sum\limits_{i=1}^{D}\frac{\left[\frac{1}{4}+r^{2}\right]\left[(\frac{1}{4}+r^{2})^{2}+\Delta
x_{i}^{2}\right]}{(\frac{1}{4}+r^{2})^{4}-\Delta x_{i}^{4}}$ (3)

and getting bounds by replacing $\Delta x_{i}$ by $r$ or 0 only in the
denominator, and applying the property that $\sum\limits_{i=1}^{D}\Delta
x_{i}^{2}=r^{2}$.

### Upper Bound:

$\displaystyle J(\Delta X)$ $\displaystyle\leq$ $\displaystyle
2\sum\limits_{i=1}^{D}\frac{\left[\frac{1}{4}+r^{2}\right]\left[(\frac{1}{4}+r^{2})^{2}+\Delta
x_{i}^{2}\right]}{(\frac{1}{4}+r^{2})^{4}-r^{4}}$ (4) $\displaystyle=$
$\displaystyle
2D\frac{\left[\frac{1}{4}+r^{2}\right](\frac{1}{4}+r^{2})^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}+2\sum\limits_{i=1}^{D}\frac{\left[\frac{1}{4}+r^{2}\right]\Delta
x_{i}^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}=2D\frac{(\frac{1}{4}+r^{2})^{3}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}+2\frac{\left[\frac{1}{4}+r^{2}\right]r^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}$
$\displaystyle=$
$\displaystyle\frac{\frac{1}{4}+r^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}\left[\frac{D}{8}+(D+2)r^{2}+2Dr^{4}\right]~{}.$

### Lower Bound:

$\displaystyle J(\Delta X)$ $\displaystyle\geq$ $\displaystyle
2\sum\limits_{i=1}^{D}\frac{\left[\frac{1}{4}+r^{2}\right]\left[(\frac{1}{4}+r^{2})^{2}+\Delta
x_{i}^{2}\right]}{(\frac{1}{4}+r^{2})^{4}}$ (5) $\displaystyle=$
$\displaystyle\frac{2D}{\frac{1}{4}+r^{2}}+2\sum\limits_{i=1}^{D}\frac{\Delta
x_{i}^{2}}{(\frac{1}{4}+r^{2})^{3}}=\frac{2D}{\frac{1}{4}+r^{2}}+\frac{2r^{2}}{(\frac{1}{4}+r^{2})^{3}}$
$\displaystyle=$
$\displaystyle\frac{1}{(\frac{1}{4}+r^{2})^{3}}\left[\frac{D}{8}+(D+2)r^{2}+2Dr^{4}\right]~{}.$

###### Theorem 2

: If $D>4$ (resp. $D<4$), $J(X)$ can attain maximum (resp. minimum) in
$\overline{\Omega}$ on its interior or boundary.

Proof:  To show that the optimum can be in the interior of
$\overline{\Omega}$, consider $\overline{\Omega}$ as a D-dimensional sphere of
a small radius $r$ centered at the origin.

Since $J(0)=8D$, by equation (4)

$\displaystyle J(\Delta X)-J(0)$ $\displaystyle\leq$
$\displaystyle\frac{\frac{1}{4}+r^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}\left[\frac{D}{8}+(D+2)r^{2}+2Dr^{4}\right]-8D$
(6) $\displaystyle=$
$\displaystyle\frac{(\frac{1}{4}+r^{2})\left[\frac{D}{8}+(D+2)r^{2}+2Dr^{4}\right]-8D((\frac{1}{4}+r^{2})^{4}-r^{4})}{(\frac{1}{4}+r^{2})^{4}-r^{4}}$
$\displaystyle=$
$\displaystyle\frac{-8Dr^{8}-6Dr^{6}+\frac{13D+4}{2}r^{4}+\frac{4-D}{8}r^{2}}{(\frac{1}{4}+r^{2})^{4}-r^{4}}~{},$

and by equation (5)

$\displaystyle J(\Delta X)-J(0)$ $\displaystyle\geq$
$\displaystyle\frac{1}{(\frac{1}{4}+r^{2})^{3}}\left[\frac{D}{8}+(D+2)r^{2}+2Dr^{4}\right]-8D$
(7) $\displaystyle=$
$\displaystyle\frac{\frac{D}{8}+(D+2)r^{2}+2Dr^{4}-8D(\frac{1}{4}+r^{2})^{3}}{(\frac{1}{4}+r^{2})^{3}}$
$\displaystyle=$
$\displaystyle\frac{-8Dr^{6}-4Dr^{4}+\frac{4-D}{2}r^{2}}{(\frac{1}{4}+r^{2})^{3}}~{}.$

To show that $J(\Delta X)<J(0)$ for $r>0$ in the interior of
$\overline{\Omega}$, it suffices to show that the polynomial at the numerator
of (6) is negative for $r>0$. The first derivative of the polynomial at $r=0$
is 0, and the second derivative is negative for $D\geq 5$. Therefore, the
polynomial is negative up to a certain value of $r$ for $D\geq 5$. To show
that $J(\Delta X)>J(0)$ for $r>0$ in the interior of $\overline{\Omega}$, it
suffices to show that the polynomial at the numerator of (6) is positive for
$r>0$. The first derivative of the polynomial at $r=0$ is 0, and the second
derivative is positive for $D\leq 3$. Therefore, the polynomial is positive up
to a certain value of $r$ for $D\leq 3$. $\Box$

###  3.2 A radius of assured decrease for $D\geq 5$, and increase for $D\leq
3$, in the counter example

Consider the counter example in Theorem 2. The proof of Theorem 2 showed that
for $D>4$, there exists a D-dimensional sphere of radius $r\leq r_{\max}$, for
some $r_{\max}$, centered at the origin on which $J(X)$ has guaranteed global
maximum in $\overline{\Omega}$ at the origin, and for $D<4$ there exist a
global minimum.

To find a lower bound for the radius of the D-dimensional sphere for $D\geq
5$, we find the smallest positive root $\bar{x}$ of the polynomial on the
numerator of equation (6) divided by $r^{2}$, where $x=r^{2}$:

$-8Dx^{3}-6Dx^{2}+\frac{13D+4}{2}x+\frac{4-D}{8}~{}=~{}0~{},$

and $r_{\max}=\sqrt{\bar{x}}$. By numerical solution we found lower bounds for
$r_{\max}$ values depicted in Table 2. A good approximation (which is a bit
lower) is obtained by ignoring $x^{2}$ and $x^{3}$ that have negative
coefficients. This yields: $r_{\max}\approx\sqrt{\frac{D-4}{52D+16}}$ which
converges to $\sqrt{\frac{1}{52}}=0.1387$ as $D\rightarrow\infty$. For
$D=1,000,000$ the root of the polynomial yields $r_{\max}=0.1400$. If we
ignore just the term $x^{3}$, we get by solving the quadratic equation
$r_{\max}=\sqrt{\frac{13D+4-\sqrt{157D^{2}+152D+16}}{24D}}$ which agrees with
the values in Table 2 to four significant digits. The limit as
$D\rightarrow\infty$ is
$\sqrt{\frac{13-\sqrt{157}}{24}}=\frac{1}{\sqrt{26+\sqrt{628}}}$.

A similar estimate for $r_{\max}$ for the minimum case can be developed by the
polynomial in the numerator of (7) yielding
$r_{\max}=\sqrt{\frac{1}{\sqrt{4D}}-\frac{1}{4}}$. For $D=2$,
$r_{\max}=0.3218$, and for $D=3$ $r_{\max}=0.1967$. Note that all these
values, for both cases, are less than 0.5 so that all the source points are
outside $\overline{\Omega}$.

Table 2: Limit on maximum $r$ $D$ | $r_{\max}$ | $D$ | $r_{\max}$
---|---|---|---
5 | 0.0603 | 10 | 0.1063
6 | 0.0783 | 15 | 0.1183
7 | 0.0892 | 20 | 0.1240
8 | 0.0966 | 50 | 0.1337
9 | 0.1021 | 100 | 0.1369

##  4 Conclusions

We considered the combined effect of multiple source points emitting some
quantity according to the inverse-square law, on an open domain that does not
include these points. We showed that the combined (total) effect of point
sources on the domain enjoys the following properties:

  1. 1.

The maximum must be on the boundary (not inside) for $D\leq 4$

  2. 2.

The minimum must be on the boundary (not inside) for $D\geq 4$

  3. 3.

The maximum may be inside the domain for $D>4$

  4. 4.

The minimum may be inside the domain for $D<4$.

Our findings enable dimensionality reduction in practical optimization
problems, resulting in significant efficiency gains and robustness.

## References

  * Axler et al., (2013) Axler, S., Bourdon, P., and Wade, R. (2013).  Harmonic function theory, volume 137.  Springer Science & Business Media.
  * Berman et al., (2010) Berman, O., Drezner, Z., and Krass, D. (2010).  Cooperative cover location problems: The planar case.  IIE Transactions, 42:232–246.
  * Drezner et al., (2020) Drezner, T., Drezner, Z., and Kalczynski, P. (2020).  Multiple obnoxious facilities location: A cooperative model.  IISE Transactions, 52:1403–1412.
  * Evans, (2022) Evans, L. C. (2022).  Partial differential equations, volume 19.  American Mathematical Society.
  * Hansen et al., (1981) Hansen, P., Peeters, D., and Thisse, J.-F. (1981).  On the location of an obnoxious facility.  Sistemi Urbani, 3:299–317.
  * Ovall, (2016) Ovall, J. S. (2016).  The Laplacian and mean and extreme values.  The American Mathematical Monthly, 123(3):287–291.
  * Stein and Shakarchi, (2010) Stein, E. M. and Shakarchi, R. (2010).  Complex analysis, volume 2.  Princeton University Press.


Energy-Efficient Precoding Designs for Multi-User Visible Light Communication Systems with Confidential Messages

Son T. Duong,
Thanh V. Pham, Member, IEEE,
Chuyen T. Nguyen,
and Anh T. Pham, Senior Member, IEEE
Son T. Duong and Chuyen T. Nguyen are with School of Electronics and Telecommunications, Hanoi University of Science
and Technology, Hanoi, Vietnam (email<EMAIL_ADDRESS>chuyen.nguyenthanh@hust.edu.vn).
Thanh V. Pham is with the Department of Electrical and Computer Engineering, McMaster University,
Hamilton L8S 4L8, Canada (e-mail: phamv12@mcmaster.ca).
Anh T. Pham is with the Department of Computer and Engineering, The University of Aizu,
Aizuwakamatsu 965-8580, Japan (e-mail: pham@u-aizu.ac.jp).
This work is supported in part by the Telecommunications Advancement Foundation (TAF) under Grant C-2020-2. This article will be presented in part at the 2021 IEEE 93rd Vehicular Technology Conference: VTC2021-Spring. (Corresponding author: Thanh V. Pham.)

This paper studies energy-efficient precoding designs for multi-user visible light communication (VLC) systems from the perspective of physical layer security where users' messages must be kept mutually confidential. For such systems, we first derive a lower bound on the achievable secrecy rate of each user. Next, the total power consumption for illumination and data transmission is thoroughly analyzed. We then tackle the problem of maximizing energy efficiency, given that each user's secrecy rate satisfies a certain threshold. The design problem is shown to be non-convex fractional programming, which renders finding the optimal solution computationally prohibitive. Our aim in this paper is, therefore, to find sub-optimal yet low complexity solutions. For this purpose, the traditional Dinkelbach algorithm is first employed to reformulate the original problem to a non-fractional parameterized one. Two different approaches based on the convex-concave procedure (CCCP) and Semidefinite Relaxation (SDR) are utilized to solve the non-convex parameterized problem. In addition, to further reduce the complexity, we investigate a design using the zero-forcing (ZF) technique.
Numerical results are conducted to show the feasibility, convergence, and performance of the proposed algorithms depending on different parameters of the system.

VLC, energy efficiency, physical layer security, Dinkelbach algorithm, convex-concave procedure, semidefinite relaxation.

§ INTRODUCTION

The exponential growth of Internet traffic over the past decade has motivated a great deal of research in new wireless technologies. In this regard, visible light communication (VLC) has been recognized as one of the most promising solutions in providing high-capacity data transmission with a license-free spectrum. Aside from this, the technology also takes advantage of the high energy-efficient light-emitting diodes (LEDs), which are widely deployed for indoor illumination. This naturally enables VLC to fit into the future ubiquitous networks.

A lot of theoretical and experimental research effort has been paid to realize practical VLC systems. Nonetheless, multiple challenges remain to make VLC more viable [1], among which security in terms of information privacy and confidentiality (especially in public areas) is one of the most critical issues. Conventional security mechanisms are mainly based on key-based cryptographic algorithms performed at upper layers of the Open Systems Interconnection (OSI) model [2]. The secrecy provided by those cryptographies comes from the assumption that potential eavesdroppers possess limited computational power, which renders obtaining the secret key in a meaningful period of time infeasible. However, with the unprecedented advances in hardware technologies and recent developments of quantum computing, traditional cryptographic techniques will no longer be secure in the foreseeable future [3]. This requires radically new approaches for communication security. In this respect, physical layer security (PLS) has emerged as a novel paradigm to replace and/or complement key-based cryptographies. The fundamental principle of PLS is to hide confidential messages from unauthorized users by exploiting the randomness of the wireless channels, noise, and interferences. The most promising implication of PLS is that perfect secured communication is achievable from the information-theoretic point of view. Moreover, this secrecy can be quantified by characterization of the secrecy rate that defines the maximum transmission rate at which unauthorized users are unable to extract any information from the received signals regardless of their computational capability.

§.§ Related Works

While PLS has been a well-investigated topic in radio frequency (RF) communications, it has only been receiving considerable attention in the past few years in the case of VLC. Compared to RF links, VLC channels exhibit several unique characteristics in terms of power constraints, making an exact expression of the secrecy rate unavailable. Considering a single-input single-output (SISO) VLC channel, which comprises a transmitter, a legitimate user, and an eavesdropper, lower and upper bounds secrecy rates were derived [4]. In the case of one transmitter and multiple receivers, including legitimate users and eavesdroppers, the non-orthogonal multiple access scheme (NOMA) was investigated to enhance the system's secrecy performance [5]. In practical VLC systems, multiple LED luminaries should usually be deployed to provide sufficient illumination. As such, multiple-input single-output (MISO) channels are more prevalent. In such scenarios, previous works mainly focused on linear precoding as a means of secrecy enhancement. Considering a system with one legitimate user and an eavesdropper, bounds on the secrecy rate of MISO VLC channels were derived in [6]. Then, optimal designs to maximize the lower bound rate were solved for the specific zero-forcing (ZF) precoding under perfect and imperfect channel state information (CSI). In the case of general precoding, the solution was investigated by transforming the non-convex design problem to a solvable quasiconvex line search problem in [7]. Regarding the scenario of MISO with one legitimate user and multiple eavesdroppers, the optimal precoding design for maximizing the minimum secrecy rate was examined in [8]. In the case of multi-user MISO with an eavesdropper, [9] studied the problem of maximizing users' secrecy sum-rate where eavesdropper's CSI can be either known (active eavesdropper) or unknown (passive eavesdropper) to the transmitters. For systems where users' messages must be kept confidential, that is, each user considers others as potential eavesdroppers, the authors in [10] provided a lower bound on the secrecy sum-rate. Solutions to the problems of max-min fairness, harmonic mean, proportional fairness, and weighted fairness were then presented. Most recently, ZF precoding strategies to simultaneously deal with active and passive eavesdroppers have been proposed and evaluated in [11]. Considering the randomness of users and eavesdroppers, the secrecy outage probability was derived in [12, 13] for the case of one transmitter and in [14] for the case of multiple ones. In addition to precoding techniques, the use of artificial noise in improving PLS has also been recently studied in [15, 16, 17, 18, 19]. While the above-mentioned works investigated the secrecy performance under the assumption that transmitters' locations are fixed, the authors in [20] studied PLS in VLC from the perspective of stochastic geometry where locations of the transmitters, legitimate users, and eavesdroppers were randomly distributed following homogeneous Poisson point processes.

In addition to the security aspect, energy efficiency (EE), which is defined as the ratio of all users' sum-rate over the total consumed power, has also been the focus of recent research in VLC. This development is propelled by the current global effort in reducing carbon footprint as it is reported that the telecommunication industry accounts for approximately 2% of global greenhouse gas emissions. In this research direction, several studies concentrated on energy-efficient designs for hybrid VLC-RF systems relevant for indoor scenarios where RF links are usually available and used for uplink transmissions [21, 22, 23, 24, 25, 26]. For the only-VLC systems, EEs of single-input single-out (SISO), MISO, and multiple-input multiple-output (MIMO) configurations were considered in [27, 28]. While these studies provided comprehensive analyses for EE in VLC, there was a lack of investigation in the context of PLS. Although there have been several studies concerning power consumption from the perspective of PLS [8, 29, 19], their objectives were to minimize the consumed power of the information-bearing signal subject to specific secrecy requirements. This design approach, however, is not necessarily optimal from the EE perspective.

§.§ Main Contributions

Considering these limitations of previous studies, in this paper, we aim at optimal precoding designs to maximize the EE of multi-user VLC systems taking into account PLS criteria. Specifically, the worst-case scenario in which each user sees others as eavesdroppers is assumed. In such a case, messages for users must be kept mutually confidential. The secrecy EE (SEE) is then defined as the ratio of the secrecy sum-rate of all users over the total consumed power. A lower bound on the users' secrecy sum-rate under this requirement was presented in [10] assuming that the noise powers at all users were the same. This assumption, nonetheless, may not be valid in practical VLC systems where receiver noise is signal-dependent.

The first contribution of this work is to generalize the result in [10] by deriving a lower bound for the case of unequal noise powers. A unique characteristic of VLC systems is that both illumination and communication are provided simultaneously.
Moreover, the allocated power for illumination does have an impact on the communication aspect.  As a result, to formulate the SEE, we present a detailed analysis for the total power consumption, which includes powers for data transmission, illumination, and circuitry operations. A precoding design problem is formulated to maximize the SEE, given that each user's secrecy rate must satisfy a predefined minimum threshold. The problem is shown to be non-convex fractional programming, which is challenging to be solved. To tackle it, we first utilize the well-known Dinkelbach algorithm to transform the original fractional objective function to a parameterized non-fractional one, which can be handled more conveniently. In the conference version of this paper [30], we proposed an approach based on the convex-concave procedure (CCCP or CCP)[31, 32] to solve the transformed problem. Specifically, a series of variable transformations and first-order Taylor approximations were employed to convexify the parameterized objective functions and the constraints. A local solution could then be efficiently found by solving the approximate problem iteratively. Since several slack variables and approximations are used to approximate the transformed problem, the involved CCCP may require a large number of iterations to converge. To reduce the complexity, we introduce another approach, which makes use of the semidefinite relaxation technique (SDR) [33, 34]. Compared to the previous approach, the use of SDR leads to fewer constraints, enabling the CCCP to converge faster. Moreover, in solving the surrogate problem in each CCCP iteration, though the second approach suffers higher worst-case complexity than that of the first one theoretically, we show that under the considered system parameters, its software implementation practically achieves lower average running times. To further reduce the complexity, we investigate a sub-optimal design using zero-forcing (ZF) precoding.

The rest of the paper is organized as follows. Section II describes the system model, including channel and signal models. In Section III, the total consumed power and a lower bound on the achievable secrecy sum-rate of the considered system are derived, which gives rise to a formulation of the SEE. Three precoding design approaches to maximize the SEE and their complexity analysis are presented in Section IV. Section V discusses representative simulation results regarding computational complexities, feasibilities, and performances of the three proposed approaches. Finally, we conclude the paper in Section VI.

Notation: The following notations are used throughout the paper. $\mathbb{R}^{m \times n}$ denotes the set of $m \times n$ real-valued matrices. The symbol $\mathbb{S}^n$ represents the set of $n \times n$ symmetric matrices. The lowercase symbols represent scalar variables. Bold uppercase symbols, e.g., $\mathbf{M}$, denote matrices, while  bold lowercase ones such as $\mathbf{v}$ represent column vectors. The $k$-th row vector of $\mathbf{M}$ is denoted as $[\mathbf{M}]_{k, :}$. The transpose of $\mathbf{M}$ is written as $\mathbf{M}^T$. The symbols $|\cdot|_1$, $\lVert\cdot\rVert$, and $|\cdot|$ are $\text{L}_1$ norm, the Euclidean norm, and absolute value operator respectively. The functions Tr($\cdot$) and det($\cdot$) represent the trace and determinant of a square matrix, respectively. Finally, $\text{diag}\left\{\mathbf{v}\right\}$ denotes the diagonal matrix whose diagonal entries are elements of vector $\mathbf{v}$.

§ SYSTEM MODEL

§.§ VLC System Description

An example of the considered system with $N_T=4$ LED luminaries and $K=3$ users.

Our considered MU-MISO VLC system, as illustrated via a simple example in Fig. <ref>, consists of $N_T$ LED luminaries, $K$ independent users where each user is equipped with a photodiode (PD). The LED luminaries are expected to transmit information to users confidentially. It is noted that the transmission is considered to be confidential if each user can decode only its own message, while it is not possible for the others to do so. The received signal at each user is from light-of-sight (LoS) and multiple non-light-of-sight (NLoS) directions. Nevertheless, thanks to the fact that LoS signal is strong enough compared with the NLoS one (which accounts for more than 95% of the total received optical power at the receiver [35]), only LoS propagation path is studied in our work to simplify the system model.

Let $h_{n,k}$ denote the channel state information (CSI) between the $n$-th LED luminary and the $k$-th user, which can be written as [35]

\begin{equation}
\label{eqn:chann_coeff}
\begin{cases}
\frac{A_r}{t_{n,k}^2}L(\phi)T_s(\psi_{n.k})g(\psi_{n,k})\cos(\psi_{n,k}) & 0 \leq \psi_{n,k} \leq \Psi, \\
0 & \psi_{n,k} > \Psi, \\
\end{cases}
\end{equation}

where $A_r$ is an active area of the PD, $t_{n,k}$ is the distance between the $n$-th LED luminary and $k$-th user. $\phi$ is the angle of irradiance, $\Psi$ is the optical field of view (FOV) of the PD, while $\psi_{n,k}$ is the angle of incidence. $T_s(\psi_{n,k})$ is the gain of optical filter. $g(\psi_{n,k})$ is the gain of the optical concentrator and is given by [35]

\begin{equation}
\label{eqn:gain_optical_concentrator}
\begin{cases}
\frac{\kappa^2}{\sin^2(\Psi)} & 0 \leq \psi_{n,k} \leq \Psi, \\
0 & \psi_{n,k} > \Psi, \\
\end{cases}
\end{equation}

where $\kappa$ is the refractive index of the concentrator. $L(\phi)$ is the emission intensity of Lambertian light source, which is calculated as [35]

\begin{equation}
\label{eqn:Lambertian}
\end{equation}

where $l$ is the order of Lambertian emission determined by

\begin{equation}
\label{eqn:order_Lambertian}
\end{equation}

where $\Theta_{0.5}$ is the LED's semi-angle for half illuminance.

§.§ Signal Model

Let $\mathbf{d} = \begin{bmatrix}d_1 & d_2 & \cdots & d_K\end{bmatrix}^T \in \mathbb{R}^{K \times 1}$ be the vector of data symbols for all users. Assume that the symbols are drawn from an $M$-ary pulse amplitude modulation ($M$-PAM) and are modeled as a random variable (RV) $d$ following a certain distribution over $[-1, 1]$ with zero-mean and variance $\sigma^2_d$.

An information-bearing signal $s_n$ for the $n$-th LED transmitter is generated from a linear combination of the data vector and a precoder $\mathbf{v}_n = \begin{bmatrix}w_{n,1} & w_{n,2} & \cdots & w_{n,K}\end{bmatrix} \in \mathbb{R}^{1 \times K}$ as

\begin{equation}
\label{eqn:information_bearing_signal}
s_n = \mathbf{v}_n \ \mathbf{d}.
\end{equation}

For illumination, an DC bias $I_{n}^{\text{DC}}$ should be added to $s_n$ to create a non-negative drive current $x_n$ for the LED. The drive current, in addition, needs to be constrained to a maximum threshold, i.e., $I_{\text{max}}$, to ensure that LEDs operate normally. Therefore,

\begin{equation}
\label{eqn:Imax}
0 \leq x_n=s_n+I_n^{\text{DC}} \leq I_{\text{max}}.
\end{equation}

The emitted optical power of each LED luminary is given by

\begin{equation}
{P}_n^s=\eta \left(s_n+I_n^{\text{DC}}\right),
\end{equation}

where $\eta$ as the LED conversion factor. Denote $\mathbf{h}_k = \begin{bmatrix}h_{1,k} & h_{2,k} & \cdots & h_{N_T,k}\end{bmatrix}^T$ is the $k$-th user's channel vector, where $h_{n, k}$ is the line-of-sight (LoS) channel coefficient between the $n$-th LED array and the $k$-user.
The electrical signal at the PD output is then given by

\begin{align}
\label{eqn:received_signal}
y_k &=  \gamma \mathbf{h}^T_k\begin{bmatrix}P_1^s & P_2^s & \cdots & P^s_{N_T,k}\end{bmatrix}  + n_k \nonumber \\
%= & \gamma H_k^T \eta u_n + n_k \\
&= \gamma \eta \left(\mathbf{h}_k^T \mathbf{w}_k d_k + \mathbf{h}_k\sum_{i=1,\ i\neq k}^{K} \mathbf{w}_i d_i + \mathbf{h}_k \mathbf{I}^{\text{DC}}\right) + n_k,
\end{align}

where $\gamma$ is the PD responsivity, $\mathbf{w_k} = \begin{bmatrix}w_{1,k} & w_{2, k} & \cdots & w_{N_T, k}\end{bmatrix}^T$ is the $k-$th user's precoder and $\mathbf{I}^{\text{DC}} = \begin{bmatrix}I_1^{\text{DC}} & I_2^{\text{DC}} & \cdots & I_{N_T}^{\text{DC}}\end{bmatrix}^T$. It is noted that since $|d_i| \leq 1$, we have

\begin{align}
\label{eqn:cond2}
-\lVert\mathbf{v}_n\rVert_1\leq s_n\leq \lVert\mathbf{v}_n\rVert_1.
\end{align}

To ensure both (<ref>) and (<ref>), the following constraint should be satisfied

\begin{equation}
\sum_{k=1}^K |w_{n,k}| \leq \min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right).
\label{eqn:linear_constraints}
%= \min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right),
\end{equation}

The receiver noise $n_k$ in (<ref>) can be modeled as a real-valued zero-mean Gaussian RV whose variance is given by [36]

\begin{equation}
\label{eqn:noise_var}
\sigma_k^2 = 2\gamma e \overline{P}_k^r B + 4\pi e A_r \gamma \chi_{\text{amb}}(1-\cos(\Psi))B + i_{\text{amp}}^2B,
\end{equation}

where $\overline{P}_k^r = \eta\mathbf{h}_k^T\mathbf{I}^{\text{DC}}$ is the average of received power at the $k$-th user, $e$ is the elementary charge, $B$ is the modulation bandwidth, $\chi_{\text{amb}}$ is the ambient light photo-current, and $i_{\text{amp}}$ is the pre-amplifier noise current density.

§ SECRECY ENERGY EFFICIENCY

We now analyze the total consumed power at LED transmitters, which can be expressed by

\begin{equation}
\label{eqn:power_consumption}
P_{\text{total}} = P_{\text{DC}}+ P_{\text{AC}},
\end{equation}

where $P_{\text{DC}}$ and $P_{\text{AC}}$ are the powers of the DC and AC currents, respectively.
The DC current power includes the power used for illumination by LEDs denoted as  $P_{\text{DC, LEDs}}$ and that used by other circuit components denoted as $P_{\text{DC, circuitry}}$. While  $P_{\text{DC, circuitry}}$ can be considered to be fixed, the power consumption of LEDs can be adjusted depending on the required dimming level. However, under a specific usage when the illumination level is not changed, it is reasonable to assume that $P_{\text{DC, LEDs}}$ is fixed as well.

The DC power consumption is then given as

\begin{equation}
\label{eqn:DC_consumption}
P_{\text{DC}} = P_{\text{DC, LEDs}}+ P_{\text{DC,circuitry}},
\end{equation}

with $P_{\text{DC, LEDs}}$ being calculated as

\begin{equation}
\label{eqn:DC_led_consumption}
P_{\text{DC, LEDs}} = \sum_{n=1}^{N_T} \ U_{\text{LED}} \ I_n^{\text{DC}},
\end{equation}

where $U_{\text{LED}}$ being the forward voltage of the LEDs.

The AC currents comes from the output current (or voltage) from the precoders of LED drivers, therefore can be calculated as

\begin{equation}
\label{eqn:LED_drive_power_consumption}
P_{\text{AC}} = r \sum_{k=1}^K \sigma_d^2 \lVert\mathbf{w}_k\rVert^2,
\end{equation}

where $r$ is the equivalent resistance of the AC circuit.
Without loss of generality, we denote $\xi = r \sigma_d^2$ as equivalent resistance, then rewrite the total power consumption as

\begin{equation}
\label{eqn:power_consumption2}
P_{\text{total}} = P_{\text{DC}} + \xi \sum_{k=1}^K \lVert\mathbf{w}_k\rVert^2.
\end{equation}

A lower bound on the achievable secrecy rate of each user was derived in [10] assuming the same noise power for all users. This assumption can be valid in the case that thermal noise is dominant over the signal-dependent shot noise. Nonetheless, at the high transmit optical power regime where shot noise is non-negligible, it is of necessity to characterize an achievable secrecy rate taking into account the signal-dependent noise when users' noise powers are not necessarily the same.

A lower bound on the achievable secrecy rate of the $k$-th user is given by

\begin{equation}
\label{eqn:secrecy_rate}
\begin{split}
R_{s,k}(\mathbf{W}) =  \frac{1}{2}\log_2\left(\frac{1+\sum_{i=1}^K a_k(\mathbf{h}_k^T \mathbf{w}_i)^2}{1+\sum_{i=1,i\neq k}^K b_k(\mathbf{h}_k^T \mathbf{w}_i)^2}\right)
- \frac{1}{2}\log_2\left(1+\sum_{i=1,i\neq k}^K b_i(\mathbf{h}_i^T \mathbf{w}_k)^2\right)
\end{split}
\end{equation}

where $\mathbf{W}=\begin{bmatrix}\mathbf{w}_1 \ \mathbf{w}_2 \ \cdots \ \mathbf{w}_K\end{bmatrix}$, $a_k=\frac{2^{2h_d}}{2\pi e \overline{\sigma}_k^2}$ and $b_k=\frac{\sigma_d^2}{\overline{\sigma}_k^2}$ with $h_d = -\int_{-1}^1f(d)\log_2f(d){\rm{d}}x$ being the differential entropy of $d_k$ and $\overline{\sigma}_k^2 = \frac{\sigma^2_k}{\left(\gamma\eta\right)^2}$.

The proof is given in the Appendix.

The SEE of the considered system is therefore calculated as

\begin{align}
\label{eqn:energy_efficiency}
\Phi(\mathbf{W})=\frac{\sum_{k=1}^K R_{s,k}(\mathbf{W})}{P_{\text{DC}} + \xi \text{Tr}\left(\mathbf{W}\mathbf{W}^T\right)}.
\end{align}

For precoding designs to maximize the SEE, it is assumed that the CSI and the noise power are estimated at the receiver (e.g., using pilot sequence) then are fed back to the LED transmitter through an uplink using Wi-Fi or infrared.

§ SEE MAXIMIZATION

§.§ Optimal Precoding

In this section, we focus on designing linear precoding schemes that maximize the  SEE defined in (<ref>)  given constraints on LEDs' input signals and targeted users' secrecy rates. Mathematically, our considered problem can be formulated as follows

𝐖maximize          Φ(𝐖)

subject to

R_s,k(𝐖)  ≥λ_k,   k = 1, 2, ⋯, K,

‖[𝐖]_n, :‖_1  ≤min(I_n^DC,I_max-I_n^DC),    n = 1, 2, ⋯, N_T,

where $\lambda_k$ is the minimum threshold for the  secrecy rate of the $k$-th user. Note that we use (<ref>) to represent (<ref>) for the sake of notational consistency.
Due to (<ref>), the above problem is a non-concave fractional programming problem, which may be difficult to solve the optimal solution.

A traditional approach to tackle fractional programming problems is to employ the Dinkelbach algorithm [37], which is proved to converge in super-linear time. The principle of the algorithm is parameterizing the objective function to get rid of its fractional form.

In particular to our problem, we first write $\Phi\left(\mathbf{W}\right) = \frac{N(\mathbf{W})}{D(\mathbf{W})}$,
where $N(\mathbf{W})=\sum_{k=1}^K R_{s,k}(\mathbf{W})$  and $D(\mathbf{W}) = P_{\text{DC}} + \xi \text{Tr}\left(\mathbf{W}\mathbf{W}^T\right)$. Then for some $\mu \geq 0$, which represents the value of $\Phi\left(\mathbf{W}\right)$, we consider the following parameterized problem

𝐖maximize          N(𝐖) - μD(𝐖)

subject to

R_s,k(𝐖)  ≥λ_k,   k = 1, 2, ⋯, K,

‖[𝐖]_n, :‖_1  ≤min(I_n^DC,I_max-I_n^DC),    n = 1, 2, ⋯, N_T.

It was shown in [37] that for some $\mu^* \geq 0$, $\mu^* = \frac{N(\mathbf{W}^*)}{D(\mathbf{W}^*)} = \underset{\mathbf{W}}{\text{maximize}}~\frac{N(\mathbf{W})}{D(\mathbf{W})}$ if and only if $\underset{\mathbf{W}}{\text{maximize}}~N(\mathbf{W}^*) - \mu^*D(\mathbf{W}^*) = 0$. Based on this, we present a Dinkelbach-type algorithm, which is described in Algorithm <ref>, for solving (<ref>).

It should be noted that the problem in (<ref>) is non-convex due to the non-concave objective function and the non-convex constraint (<ref>). In the following, we propose two approaches to solve the problem. The first approach relies on several variable transformations and first-order Taylor approximations to approximately convexify (<ref>). The CCCP is then utilized to successively solve the approximated problem. In the second approach, we make use of the SDR method to reduce the number of slack variables and approximations needed to approximate the original problem. The approximated problem, in this case, is formulated as a DC programming problem, which can also be solved efficiently using the CCCP.

Choose the maximum number of iterations $L_{\text{max}, 1}$ and the error tolerance $\epsilon_{1}$.

Initialize $\mu>0$, $l=0$ .

$\text{convergence}==\bold{False} \ \text{and} \ l \leq {L}_{\text{max}, 1}$
For a given $\mu$, solve (<ref>) to get $\mathbf{W}^{(l)}$.

$N\left(\mathbf{W}^{(l)}\right)-\mu D\left(\mathbf{W}^{(l)}\right) \leq \epsilon_{1}$

Update $\mu = \frac{N\left(\mathbf{W}^{(l)}\right)}{D\left(\mathbf{W}^{(l)}\right)}$;

Return optimal ${\mathbf{W}^*}$ and ${\mu^*}$.
Dinkelbach-type algorithm for solving (<ref>)

§.§.§ CCCP

We first notice that one can transform (<ref>) to be concave by introducing the following slack variables

\begin{align}
r_{1,k} \ \overset{\Delta}{=} \ & \frac{1}{2}\log_2\left(1+\sum_{i=1}^K a_k\left(\mathbf{h}_k^T \mathbf{w}_i\right)^2\right) \label{eqn:r1k} ,\\
p_{1,k} \ \overset{\Delta}{=} \ & \sum_{i=1}^K a_k\left(\mathbf{h}_k^T \mathbf{w}_i\right)^2 \label{eqn:p1k} ,\\
r_{2,k} \ \overset{\Delta}{=} \ & \frac{1}{2}\log_2\left(1+\sum_{i=1,i\neq k}^K b_k\left(\mathbf{h}_k^T \mathbf{w}_i\right)^2\right) \label{eqn:r2k} ,\\
p_{2,k} \ \overset{\Delta}{=} \ & \sum_{i=1,i\neq k}^K b_k\left(\mathbf{h}_k^T \mathbf{w}_i\right)^2 \label{eqn:p2k} ,\\
r_{3,k} \ \overset{\Delta}{=} \ & \frac{1}{2}\log_2\left(1+\sum_{i=1,i\neq k}^K b_i\left(\mathbf{h}_i^T \mathbf{w}_k\right)^2\right) \label{eqn:r3k} ,\\
p_{3,k} \ \overset{\Delta}{=} \ & \sum_{i=1,i\neq k}^K b_i\left(\mathbf{h}_i^T \mathbf{w}_k\right)^2 \label{eqn:p3k} .
\end{align}

Accordingly, (<ref>) can be equivalently reformulated as

r_1, k,r_2, k, r_3, k,
p_1, k, p_2, k, p_3, kmaximize         ∑_k=1^K(r_1,k - r_2,k - r_3,k)-μ(P_DC+ϵTr(𝐖𝐖^T))

subject to

r_1,k ≤1/2log_2(1+p_1,k),

p_1,k ≤∑_i=1^K a_k(𝐡_k^T 𝐰_i)^2,

r_2,k ≥1/2log_2(1+p_2,k),

p_2,k ≥∑_i=1,i≠k^K b_k(𝐡_k^T 𝐰_i)^2,

r_3,k ≥1/2log_2(1+p_3,k),

p_3,k ≥∑_i=1,i≠k^K b_i(𝐡_i^T 𝐰_k)^2,

r_1,k - r_2,k - r_3,k ≥λ_k,

‖[𝐖]_n, :‖_1  ≤min(I_n^DC,I_max-I_n^DC).

In the construction of the above problem, equalities in (<ref>)-(<ref>) are replaced with their corresponding inequalities constraints (<ref>)-(<ref>), respectively. These replacements, however, do not change the optimality of the optimization problem. Indeed, as (<ref>) monotonically increases with $r_{1, k}$, at the optimal solution, $r_{1, k}$ must achieve its maximum value. This implies that both (<ref>) and (<ref>) then hold with equality. Similarly, the validity of (<ref>)-(<ref>) can be verified.

It is seen that constraints (<ref>), (<ref>), (<ref>), and (<ref>) are convex but not (<ref>), (<ref>), and (<ref>). To cope with this situation, the first-order Taylor expansions are employed to approximately linearize these non-convex constraints. The CCCP is then utilized to successively solve the approximated problem. Specifically, the proposed CCCP involves an iterative process, where at the $m$-th iteration, the following Taylor approximations are used

\begin{align}
\sum_{i=1}^K a_k\left(\left(\mathbf{h}_k^T {\mathbf{w}}^{(m-1)}_i\right)^2 + 2\left[{\mathbf{w}}^{(m-1)}_i\right]^T \mathbf{h}_k \mathbf{h}_k^T \left(\mathbf{w}_i -{\mathbf{w}}^{(m-1)}_i\right)\right)
\leq \sum_{i=1}^K a_k\left(\mathbf{h}_k^T \mathbf{w}_i\right)^2,
\label{Taylor1}
\end{align}

\begin{align}
\frac{1}{2}\log_2\left(1+{p}^{(m-1)}_{2,k}\right) + \frac{\left(p_{2,k}-{p}^{(m-1)}_{2,k}\right)}{2\text{ln}(2)\left(1+{p}^{(m-1)}_{2,k}\right)} \geq \frac{1}{2}\log_2\left(1 + p_{2, k}\right),
\label{Taylor2}
\end{align}

\begin{align}
\frac{1}{2}\log_2\left(1+{p}^{(m-1)}_{3,k}\right) + \frac{\left(p_{3,k}-{p}^{(m-1)}_{3,k}\right)}{2\text{ln}(2)\left(1+{p}^{(m-1)}_{3,k}\right)} \geq \frac{1}{2}\log_2\left(1 + p_{3, k}\right).
\label{Taylor3}
\end{align}

Here, $\mathbf{w}^{(m-1)}_i$, $p_{2, k}^{(m-1)}$, and $p_{3, k}^{(m-1)}$ are the values of $\mathbf{w}$, $p_{2, k}$, and $p_{3, k}$ obtained at the $(m-1)$-th iteration, respectively. Now, replacing the non-linear terms in (<ref>), (<ref>), and (<ref>) by their above approximate linear representations results in

r_1, k, r_2, k, r_3, k,
p_1, k, p_2, k, p_3, kmaximize  ∑_k=1^K(r_1,k - r_2,k - r_3,k)-μ(P_DC+ϵTr(𝐖𝐖^T))

subject to

p_1, k ≤∑_i=1^K a_k((𝐡_k^T 𝐰^(m-1)_i)^2 + 2[𝐰^(m-1)_i]^T 𝐡_k 𝐡_k^T (𝐰_i -𝐰^(m-1)_i)),

r_2, k ≥1/2log_2(1+p^(m-1)_2,k) + (p_2,k-p^(m-1)_2,k)/2ln(2)(1+p^(m-1)_2,k),

r_3, k ≥1/2log_2(1+p^(m-1)_3,k) + (p_3,k-p^(m-1)_3,k)/2ln(2)(1+p^(m-1)_3,k),

(<ref>), (<ref>), (<ref>),  (<ref>),

which is a convex optimization and thus can be solved efficiently by using a standard convex solver, such as CVX toolbox [38]. Finally, (<ref>) can be solved by the proposed Algorithm <ref> described  as follows

Choose the maximum number of iteration ${L}_{\text{max}, 2}$ and the error tolerance $\epsilon_{2}>0$.

Choose feasible initial points $\mathbf{W}^{(0)}$, $p^{(0)}_{2,k}$, $p^{(0)}_{3,k}$ to problem (<ref>).

Set $m=0$

$\text{convergence}==\mathbf{False}$ and $m \leq {L}_{\text{max}, 2}$
Use $\mathbf{W}^{(m-1)}$, $p^{(m-1)}_{2,k}$, $p^{(m-1)}_{3,k}$ obtained from the previous iteration to formulate the problem (<ref>).

Solve (<ref>) to get $\mathbf{W}^{(m)}$, $p^{(m)}_{2,k}$, $p^{(m)}_{3,k}$.

$\frac{\norm{\mathbf{W}^{(m)} - \mathbf{W}^{(m-1)}}}{\norm{\mathbf{W}^{(m)}}} \leq \epsilon_{2}$ and $\frac{\left|p^{(m)}_{2, k} - p^{(m-1)}_{2, k}\right|}{p^{(m)}_{2, k}} \leq \epsilon_{2}$ and $\frac{\left|p^{(m)}_{3, k} - p^{(m-1)}_{3, k}\right|}{p^{(m)}_{3, k}} \leq \epsilon_{2}$
$\text{convergence} = \mathbf{True}$

$\mathbf{W}^* = \mathbf{W}^{(m)}$

$p_{2, k}^* = p^{(m)}_{2, k}$

$p_{3, k}^* = p^{(m)}_{3, k}$

$\text{convergence} = \mathbf{False}$

$m = m + 1$.

Return the optimal value $\mathbf{W}^*$, $p_{2, k}^*$, and $p_{3, k}^*$.
CCCP-type algorithm to solve (<ref>).

§.§.§ CCCP combined with SDR

In this section, the SDR method is used to approximate (<ref>) to a convex problem.
In particular, we first define $\mathbf{P}_k=\mathbf{h}_k \mathbf{h}_k^T$ and $\mathbf{Q}_k=\mathbf{w}_k \mathbf{w}_k^T$. This variable transformation requires that $\mathbf{Q}_k \succeq \mathbf{0}$ and $\text{rank}\left(\mathbf{Q}_k\right) = 1$.
Then, based on the fact that $\left(\mathbf{h}_i^T \mathbf{w}_j\right)^2=\text{Tr}\left(\mathbf{h}_i \mathbf{h}_i^T \mathbf{w}_j \mathbf{w}_j^T\right)=\text{Tr}(\mathbf{P}_i \mathbf{Q}_j)$, the secrecy rate of the $k$-th user in (<ref>) can be rewritten as

\begin{equation}
\label{eqn:EE_DC_equattion}
R_{s,k}\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) = f_k\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) - g_k\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) ,
\end{equation}

\begin{equation}
f_k\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) =  \frac{1}{2}\log_2\left(1+\sum_{i=1}^K a_k \text{Tr}\left(\mathbf{P}_k \mathbf{Q}_i\right)\right)  \label{eqn:f_k},
\end{equation}

\begin{equation}
g_k\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) =   \frac{1}{2}\log_2\left(1+\sum_{i=1,i\neq k}^K b_k \text{Tr}\left(\mathbf{P}_k \mathbf{Q}_i\right)\right)
+ \frac{1}{2}\log_2\left(1+\sum_{i=1,i\neq k}^K b_i \text{Tr}\left(\mathbf{P}_i \mathbf{Q}_k\right)\right) . \label{eqn:g_k}
\end{equation}

Through the definition of $\mathbf{Q}_k$'s, we can write (<ref>) as $\Upsilon\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right) = \sum_{k = 1}^K\left( f_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right) \right.$ $- \left. g_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right)\right) - \mu\left(P_{\text{DC}} + \epsilon\sum_{k = 1}^K\text{Tr}\left(\mathbf{Q}_k\right)\right)$. Since $\sum_{k = 1}^K\text{Tr}\left(\mathbf{Q}_k\right)$ is affine, $\Upsilon\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right)$ can be seen as the difference of two concave functions. Moreover, the left-hand side of (<ref>) is rewritten as $f_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right) - g_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right)$, which is also the difference of two concave functions. We, therefore, recast (<ref>) in the form of DC programming, which can be solved efficiently using the CCCP [39]. Specifically, at the $m$-th iteration of the procedure, the concave term $g_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right)$ is approximately linearized as

\begin{align}
{g}_k\left(\mathbf{Q}_1, \cdots,\mathbf{Q}_K\right) \approx  g_k\left({\mathbf{Q}}_1^{(m-1)}, \cdots, {\mathbf{Q}}_K^{(m-1)}\right)  + \sum_{i=1}^K \left\langle \nabla_{g_k}\left({\mathbf{Q}}_i^{(m-1)}\right),\mathbf{Q}_i-{\mathbf{Q}}_i^{(m-1)} \right\rangle_{\text{F}},
\label{eqn:g_k_approximation}
\end{align}

where $\langle \cdot, \cdot \rangle_{\text{F}}$ denotes the Frobenius inner dot product, $\mathbf{Q}^{(m-1)}_i$ is the value of $\mathbf{Q}_i$ at the $(m-1)$-th iteration, $\nabla_{g_k}\left({\mathbf{Q}}_i^{(m-1)}\right)$ is the gradient of $g_k\left(\mathbf{Q}_1, \cdots, \mathbf{Q}_K\right)$ at point ${\mathbf{Q}}_i^{(m-1)}$ and calculated as

\begin{align}
\nabla_{g_k}\left(\mathbf{Q}_i^{(m-1)}\right) =
\begin{cases}
z^1_k b_k \ \mathbf{P}_k^T & i \neq k, \\
z^2_k \sum_{j=1,j \neq k}^K b_j \mathbf{P}_j^T & i=k, \\
\end{cases}
\label{eqn:grad_g_k}
\end{align}

where $z_k^1$ and $z_k^2$ are given by

\begin{align}
z_k^1 = \frac{1}{2\ln(2)\left(1+\sum_{j=1,j\neq k}^K b_k \text{Tr}\left(\mathbf{P}_k {\mathbf{Q}}_j^{(m-1)}\right)\right)} , \label{eqn:z1k}\\
z_k^2 = \frac{1}{2\ln(2)\left(1+\sum_{j=1,j\neq k}^K b_j \text{Tr}\left(\mathbf{P}_j {\mathbf{Q}}_k^{(m-1)}\right)\right)}. \label{eqn:z2k}
\end{align}

Now, we need to represent the constraint $\left\lVert\left[\mathbf{W}\right]_{n, :}\right\rVert_1 = \sum_{k=1}^K \left|w_{n,k}\right|  \leq \min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)$ as a function of $\mathbf{Q}_k$'s. Unfortunately, such an equivalent reformulation is not possible. To overcome this issue, previous studies [19] [40] replaced the constraint with a stricter one using

\begin{align}
\frac{\left(\sum_{k= 1}^K\left|w_{n, k}\right|\right)^2}{K} \leq \sum_{k=1}^K w^2_{n, k} = \sum_{k = 1}^K [\mathbf{Q}_k]_{n, n} \leq \frac{1}{K} \left(\min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)\right)^2.
\end{align}

Using the stricter constraint, however, results in a smaller feasible region. Hence, the optimal objective value of the new problem may be smaller than that of the original one.
On the other hand, we can also replace the original constraint with a looser one as [41][42]

\begin{align}
\sum_{k=1}^K w^2_{n, k} \leq \left(\sum_{k= 1}^K\left|w_{n, k}\right|\right)^2 \leq \left(\min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)\right)^2.
\end{align}

This approach, nonetheless, leads to a larger feasible region. Thus, the optimal point of the new problem may lie outside the feasible region of the original problem.

To address these issues, we introduce a new approach by using a stricter constraint, whose feasible region progressively approaches the optimal point of the original problem. Specifically, we make use of the following inequality

\begin{align}
\label{eqn:quadratic_maximum_constraints}
\sum_{k=1}^K \frac{w_{n,k}^2}{\delta_{n,k}}  \leq \frac{1}{\sum_{k=1}^K \delta_{n,k}} \left(\min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)\right)^2.
\end{align}

As proof, using Cauchy-Schwarz inequality, one has

\begin{align}
\label{eqn:Cauchy-Schwarz-inequality}
\sum_{k=1}^K |w_{n,k}|  \leq \sqrt{\left(\sum_{k=1}^K \frac{w_{n,k}^2}{\delta_{n,k}}\right)\left(\sum_{k=1}^K \delta_{n,k}\right)}
& \leq \sqrt{\frac{\left(\min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)\right)^2}{\sum_{k=1}^K \delta_{n,k}}\left(\sum_{k=1}^K \delta_{n,k}\right)} \nonumber \\
& = \min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right).
\end{align}

Obviously, setting $\delta_{n,k} = 1$ results in the constraint used in [19] [40]. To make the feasible region of (<ref>) approach the optimal point, $\delta_{n,k}$ is updated by the optimal precoder value from the previous iteration, i.e.,  $\delta_{n,k} = |w_{n,k}^{(m-1)}|$. The intuition behind this is that after each iteration, the feasible region of the updated constraint (a $K$-dimensional ellipsoid) extends (or contracts) along the direction of $w_{n,k}$ corresponding to the obtained $w_{n,k}^{(m-1)}$ from the previous iteration. Hence, when the iteration terminates at the optimal solution, the feasible region of the updated constraint must contain the optimal point. Also, the constraint in (<ref>) needs to be approximately represented in terms of $\mathbf{Q}_k$'s as well. It is first seen that $|w_{n,k}|=\sqrt{w_{n,k}^2}= \sqrt{\text{Tr}(\mathbf{E}_{n} \mathbf{Q}_k)}$, where $\mathbf{E}_{n} \in \mathbb{R}^{N_T\times N_T}$ denotes the all-zero matrix except the $(n,n)$-th entry is 1. Then, the constraint can be approximated as follows

\begin{equation}
\sum_{k=1}^K \frac{\text{Tr}\left(\mathbf{E}_{n}\mathbf{Q}_k\right)}{\sqrt{\text{Tr}\left(\mathbf{E}_{n} {\mathbf{Q}}^{(m-1)}_k\right)}}  \leq \frac{\left(\min\left(I_n^{\text{DC}},I_{\text{max}}-I_n^{\text{DC}}\right)\right)^2}{\sum_{k=1}^K \sqrt{\text{Tr}\left(\mathbf{E}_{n} {\mathbf{Q}}^{(m-1)}_k\right)}} . \label{eqn:quadratic_maximum_constraints_3}
\end{equation}

It is observed that the only non-convex terms now is $\text{rank}(\mathbf{Q}_k) = 1$. The SDR technique works by omitting these rank-one constraints, resulting in the following convex problem.

𝐐_1, ⋯, 𝐐_Kmaximize  ∑_k=1^K(f_k(𝐐_1, ⋯, 𝐐_K) - g_k(𝐐_1^(m-1), ⋯, 𝐐_K^(m-1))
- ∑_i = 1^K⟨∇_g_k(𝐐_i^(m-1)),𝐐_i - 𝐐_i^(m-1)⟩_F )
- μ(P_DC+ϵ∑_k = 1^KTr(𝐐_k))

subject to

f_k(𝐐_1, ⋯, 𝐐_K) - g_k(𝐐_1^(m-1), ⋯, 𝐐_K^(m-1)) - ∑_i = 1^K⟨∇_g_k(𝐐_i^(m-1)),𝐐_i - 𝐐_i^(m-1)⟩_F ≥λ_k,

∑_k=1^K Tr(𝐄_n𝐐_k)/√(Tr(𝐄_n 𝐐^(m-1)_k))  ≤(min(I_n^DC,I_max-I_n^DC))^2/∑_k=1^K √(Tr(𝐄_n 𝐐^(m-1)_k)),

𝐐_k ≽0.

Accordingly, a CCCP-type algorithm for solving (<ref>) is described in Algorithm <ref>.

Choose the maximum number of iteration ${L}_{\text{max}, 3}$ and the error tolerance $\epsilon_3>0$.

Choose an initial point $\mathbf{W}^{(0)}$ being feasible to the problem (<ref>) and calculate corresponding matrices $\mathbf{Q}^{(0)}_1, \cdots, \mathbf{Q}^{(0)}_K$.

Set $m=0$.

$\text{Convergence}==\mathbf{False}$ and $m \leq {L}_{\text{max}, 3}$
Use $\mathbf{Q}_1^{(m-1)}$, $\cdots$,$\mathbf{Q}_K^{(m-1)}$ obtained from the previous iteration to formulate the problem (<ref>).

Get $\mathbf{Q}_1^{(m)*}$, $\cdots$,$\mathbf{Q}_K^{(m)*}$ as the optimal solutions to (<ref>).

Obtain ${\mathbf{Q}}_1^{(m)}$, $\cdots$,${\mathbf{Q}}_K^{(m)}$ as rank-one approximations to $\mathbf{Q}_1^{(m)*}$, $\cdots$,$\mathbf{Q}_K^{(m)*}$, respectively.

Retrieve $\mathbf{W}^{(m)}$ from $\mathbf{Q}_1^{(m)}$, $\cdots$,$\mathbf{Q}_K^{(m)}$ .

$\frac{\norm{\mathbf{W}^{(m)} - \mathbf{W}^{(m-1)}}}{\norm{\mathbf{W}^{(m)}}} \leq \epsilon_3$
$\text{Convergence} = \mathbf{True}$.

$\mathbf{W}^* = \mathbf{W}^{(m)}$.

$\text{Convergence} = \mathbf{False}$.

$m = m + 1$.
Return the optimal solution $\mathbf{W}^*$.
CCCP-type algorithm with SDR for solving (<ref>)

It should be noted that the optimal solutions to (<ref>) in each iteration may not satisfy the rank-one constraints. Thus, in the proposed Algorithm <ref>, the optimal solution $\mathbf{Q}_k^{(m)*}$ at the $m$-th iteration is approximated by a rank-one matrix ${\mathbf{Q}}_k^{(m)}$, which is given by

\begin{align}
{\mathbf{Q}}_k^{(m)} = \Lambda^{(m)}_{\text{max}, k} \mathbf{q}^{(m)}_{\text{max}, k}\left[\mathbf{q}^{(m)}_{\text{max}, k}\right]^T,
\label{rank-oneApproximation}
\end{align}

where $\Lambda_{\text{max}, k}$ is the maximum eigenvalue of $\mathbf{Q}^{(m)*}_k$ and $\mathbf{q}^{(m)}_{\text{max}, k}$ is the eigenvector which is associated with $\Lambda_{\text{max}, k}$.

§.§ ZF Precoding

To further simplify the precoding design, we investigate the use of ZF technique in this section. Under the ZF constraint, the precoder of a user is constructed in such as way that it is orthogonal to other users' channels, i.e., $\mathbf{h}_i^T\mathbf{w}_k = 0~~\forall i \neq k$, thus eliminating all inter-user interference. Due to this orthogonality, the degrees of freedom in designing $\mathbf{W}$ can considerably decrease compared to the general case, resulting in a lower computational complexity at the expense of reduced performance.

If we denote $\mathbf{H} = \begin{bmatrix}\mathbf{h}_1 & \mathbf{h}_2 & \cdots & \mathbf{h}_K\end{bmatrix}^T$, the ZF constraint can then be represented as

\begin{align}
\mathbf{H}\mathbf{W} = \text{diag}\left\{\sqrt{\pmb{\rho}}\right\},
\label{ZFConstraint}
\end{align}

where $\pmb{\rho} = \begin{bmatrix}\rho_1 & \rho_2 & \cdots & \rho_K\end{bmatrix}^T$ is the vector of users' channel gains.
As a result, the secrecy rate of the $k$-th user is simplified to

\begin{equation}
R_{s,k}(\mathbf{W}) = \frac{1}{2}\log_2\left(1+ a_k\rho_k\right).
\end{equation}

Then the problem (<ref>) can be rewritten as

𝐖, ρmaximize          ∑_k=1^K 1/2log_2(1+ a_kρ_k) - μ(P_DC + ξTr(𝐖𝐖^T))

subject to

1/2log_2(1 + a_kρ_k)  ≥λ_k,

‖[𝐖]_n, :‖_1  ≤min(I_n^DC,I_max-I_n^DC),

𝐇𝐖 = diag{√(ρ)}.

The above problem is not convex due to the non-convex constraint in $\eqref{eqn:constraint63}$. Hence, we again employ the CCCP to solve the problem by approximately convexing  $\eqref{eqn:constraint63}$ using the first-order Taylor approximation as

\begin{align}
\mathbf{H}\mathbf{W} \approx \text{diag}\left\{\sqrt{\pmb{\rho}^{(m-1)}} + \frac{1}{2\sqrt{\pmb{\rho}^{(m-1)}}}\left(\pmb{\rho} - \pmb{\rho}^{(m-1)}\right)\right\},
\label{eqn:constraint63_approximate}
\end{align}

where $\pmb{\rho}^{(m-1)}$ is the vector of users' channel gains at the $(m-1)$-th iteration.
Finally, (<ref>) can be solved by the following algorithm

Choose the maximum number of iteration ${L}_{\text{max}, 4}$ and the error tolerance $\epsilon_{4}>0$.

Choose an initial point $\pmb{\rho}^{(0)}$ that is feasible to the problem (<ref>).

Set $m=0$

$\text{convergence}==\mathbf{False}$ and $m \leq {L}_{\text{max}, 4}$
Use $\pmb{\rho}^{(m-1)}$ as previous optimal points to formulate (<ref>) where (<ref>) is replaced with (<ref>).

Get $\mathbf{W}^{(m)}$ and $\pmb{\rho}^{(m)}$ as optimal solutions to (<ref>).

$\frac{\norm{\pmb{\rho}^{(m)} - \pmb{\rho}^{(m-1)}}}{\norm{\pmb{\rho}^{(m)}}} \leq \epsilon_4$
$\text{convergence} = \mathbf{True}$

$\pmb{\rho}^* = \pmb{\rho}^{(m)}$

$\mathbf{W}^* = \mathbf{W}^{(m)}$

$\text{convergence} = \mathbf{False}$

$m = m+1$.

Return the optimal values $\mathbf{W}^*$ and $\pmb{\rho}^*$.
CCCP-type algorithm for solving (<ref>).

§.§ Complexity Analysis

All the three proposed precoding schemes to solve (<ref>) utilize the Dinkelbach algorithm (i.e., Algorithm <ref>) as the outer loop and each of $\textbf{Algorithms \ref{alg.2}}$, $\textbf{\ref{alg.3}}$, and $\textbf{\ref{alg.4}}$ as the inner loop. Assume that in each iteration of the inner loop, the Interior Point Algorithm (IPA) is employed to solve its convex sub-problem (i.e., (<ref>), (<ref>), and (<ref>)). As a result, the number of Newton steps counted from the starting to the optimal points is used as the complexity measure. We estimate the total number of Newton steps $N_{s}^{\text{total}}$ being involved in each scheme as follows

\begin{equation}
\label{eqn:newton_steps}
N_s^{\text{total}}=N_{\text{iter}} \times N_s,
\end{equation}

where $N_{\text{iter}}$ is the number of total iterations required for the algorithm's convergence and $N_s$ is the number of Newton steps needed to solve the involved convex sub-problem. It is noted that $N_{\text{iter}}$ is counted from the start to the convergence of the outer loop. In this section, we qualitatively characterize $N_s$ while $N_{\text{iter}}$ is numerically investigated in the next section. According to [43] and [44], the worst-case $N_s$ required for a non-linear convex problem to find one of its local solutions can be estimated by

\begin{equation}
\label{eqn:newton_steps2}
N_s \sim \sqrt{\text{problem size}},
\end{equation}

where “problem size" is the number of scalar variables of the considered optimization problem. Since there are $(N_T + 6)K$, $N^2_TK$, and $(N_T + 1) K$ scalar variables in (<ref>), (<ref>), and (<ref>), respectively,
we have $N_s^{\textbf{Algorithm 2}} \sim \sqrt{(N_T + 6)K}$, $N_s^{\textbf{Algorithm 3}} \sim \sqrt{N_T^2K}$, and $N_s^{\textbf{Algorithm 4}} \sim \sqrt{(N_T + 1)K}$. It is obvious that the use of ZF constraint in (<ref>) reduces  the design complexity, especially when $N_T$ is small. In addition, we have

\begin{align}
\frac{N_s^{\textbf{Algorithm 2}}}{N_s^{\textbf{Algorithm 3}}} & \sim \sqrt{\frac{N_T+6}{N_T^2}},
%\frac{N_s^{\text{CCP-ZF}}}{N_s^{\text{SDR-general}}} & \sim \sqrt{\frac{N_T+2}{N_T^2}} \sim \frac{1}{\sqrt{N_T}}
\end{align}

which indicates that in the worst-case scenario, the complexity of solving each iteration of Algorithm 3 is higher than that of Algorithm 2 theoretically. Also, this increased complexity is proportional to the number of LED luminaries $N_T$.
However, in the next section, we show that under our parameter setting when $N_T$ is set practically small, i.e., $N_T \leq 9$, the software implementation of Algorithm 3 achieves lower average running time than that of Algorithm <ref> does.

§ SIMULATION RESULTS AND DISCUSSIONS

In this section, numerical results regarding feasibilities, convergence behaviors, and performances of the three design approaches are presented. To specify positions of LED luminaries and users, a Cartesian coordinate system whose origin is the center of the floor is used. Additionally, users are placed on the same plane at the height of $0.5$ m as described in Fig. <ref>.

If not otherwise specified, we also assume that $N_T = 4$, the average emitted optical power of each LED luminary $\overline{P}_n^s = \eta I_{n}^{\text{DC}} = 30$ dBm, $P_{\text{DC,circuitry}}=8 ~\text{W}$, $\xi=3~\Omega$, $\lambda_k$'s $=0.5$ bits/s/Hz. Other parameters for simulations are summarized in Table <ref> [15]. For an explicit expression of the lower bound in (<ref>), users' data symbols $d_k$'s are assumed to be uniformly distributed over $[-1, 1]$.

System Parameters

2cRoom and LED configurations

Room dimension  (Length $\times$ Width $\times$ Height)      5 (m) $\times$ 5 (m) $\times$ 3 (m)

LED positions     luminary 1 : $\left(-\sqrt{2}, -\sqrt{2}, 3\right)$, luminary 2 : $\left(\sqrt{2}, -\sqrt{2}, 3\right)$
luminary 3 : $\left(\sqrt{2}, \sqrt{2}, 3\right)$, luminary 4 : $\left(-\sqrt{2}, \sqrt{2}, 3\right)$

LED bandwidth, $B$     20 MHz

LED beam angle, $\phi$     $120^\circ$  (LED Lambertian order is 1)

LED conversion factor, $\eta$     2.00 W/A

2cReceiver photodetectors

Active area, $A_r$     1 $\text{cm}^2$

Responsivity, $\gamma$     0.54 A/W

Field of view (FOV), $\Psi$     $60^\circ$

Optical filter gain, $T_s(\psi)$     1

Refractive index of the concentrator, $\kappa$     1.5

2cOther parameters

Ambient light photocurrent, $\chi_{\text{amp}}$     10.93 $\text{A}/(\text{m}^2 \cdot \text{Sr}$)

Preamplifier noise current density, $i_{\text{amp}}$     5 $\text{pA}/\text{Hz}^{-1/2}$

Feasibility probability of proposed designs versus minimum threshold of secrecy rate.

First, the feasibility of the three design approaches is studied. For this purpose, 10,000 channel realizations are randomly generated. Then, the feasibility probability for each approach is calculated based on the number of realizations that the corresponding optimization problem is feasible.
With respect to the minimum threshold of the secrecy rate $\lambda_k$, we illustrate in Fig. <ref> the feasibility probabilities for different values of the average optical power $\overline{P}_n^s$. Essentially, this is to evaluate the influences of (<ref>) and (<ref>) on the feasibility of the design problem. Our obtained result shows almost no difference among the three designs, especially for larger values of $\overline{P}_n^s$. The slightly lower feasibility probability of the ZF design is due to the additional constraint in (<ref>).

We also observed that to ensure $\lambda_k$'s $=0.5$ bits/s/Hz and a feasibility of above 90%, $\overline{P}_n^s$ should be chosen to be 30 dBm or higher.

Furthermore, the feasibility probability is investigated for different numbers of LED transmitters and users. To do so, we examine three different configurations of LED luminaries, which are $2 \times 2$, $2 \times 3$, and $3 \times 3$ as illustrated in Fig. <ref>. The results, as displayed in Fig. <ref>, clearly show that the feasibility increases following an increase in the number of transmitters and a decrease in the number of users. Here, to guarantee a feasibility of above 90%, the following settings can be employed: $(N_T, K ) = (4, 3)$, $(6, 4)$, and $(9, 6)$.

$2\times2$ layout.

$2\times3$ layout.

$3\times3$ layout.

Different layouts of LED transmitters.

Feasibility probability of proposed designs among different numbers of $N_T$ and $K$.

Second, we evaluate in Fig. <ref> the convergence behaviors of three precoding schemes. Since all these proposed schemes employ Algorithm <ref> as the outer loop, we differentiate them by referring to the algorithms used in the inner loop (i.e., Algorithms <ref>, <ref>, or <ref>). In this simulation, the results are obtained through averaging those from 10000 randomly generated channel realizations and are normalized with respect to the optimal value solved by CCCP.

Furthermore, the maximum numbers of iterations in all scenarios are set to be 100. Because the size of the system (i.e., numbers of transmitters and users) influences the convergence of the proposed algorithms as analyzed in Section <ref>, we examine three scenarios: $(N_T, K) = (4, 3)$, $(6, 4)$, and $(9, 6)$,  which guarantee at least 90% feasibility.

Convergence behaviors of proposed algorithms for different numbers of $N_T$ and $K$.

Average execution time (seconds) of each iteration.

SettingApproach     CCCP (Algorithm <ref>)      CCCP+SDR (Algorithm <ref>)    ZF  (Algorithm <ref>)

$N_T=4$ & $K=3$      0.4975     0.3473     0.3236

$N_T=6$ & $K=4$      0.5996     0.4328     0.3756

$N_T=9$ & $K=6$      0.7884     0.6437     0.4594

For the CCCP (Algorithm <ref>) and the CCCP combined with SDR  approaches (Algorithm <ref>), we first examine the case that the initial precoder $\mathbf{W}^{(0)}$ is chosen to be random. With this initialization, it is seen that both algorithms require considerably large numbers of iterations for the normalized energy efficiency to converge (i.e., approximately 60, 80, and 100 iterations for the $(N_T, K) = $(4, 3), (6, 4), and (9, 6) settings, respectively).
Previous studies usually assumed random initial points, which might result in a large number of iterations, especially in our considered problem whose solution involves a nested iterative procedure. As we notice from the simulation results that the optimal solutions to (<ref>) and (<ref>) are usually near-ZF precoders, one thus can choose the initial point as a ZF precoder to speed up the convergence rate. Using this initialization strategy, we observe significant improvements in the convergence rate. For $(N_T, K) = $(4, 3), (6, 4), and (9, 6), the CCCP approach needs about 8, 15, and 30 iterations, respectively. For the case of CCCP combined with SDR, these required iterations are all approximately 5. In terms of required CCCP iterations, this confirms the low complexity offered by using the SDR technique. However, one should emphasize that this reduced complexity comes from the relaxation of the rank constraints, which renders the approach an approximation in general.
Moreover, it is shown in Table <ref> that in practice for each iteration, the $\textbf{Algorithm \ref{alg.3}}$ achieves relatively lower average execution times compared to those of the $\textbf{Algorithm \ref{alg.2}}$ over all three settings[Simulations were conducted using Matlab 2015a, CVX with Mosek solver version 9.1.9 running on a Windows 10 computer with Intel Xeon E5-1603 v3, 16 GB RAM.].
In the case of ZF precoding, the same convergence behavior to that of the CCCP combined with the SDR approach using the ZF initial point is seen. Also, its average execution time is shown to be the lowest.

SEE versus average optical power $\overline{P}_n^s$ of each LED luminary.

In Fig. <ref>, we investigate the average SEE with respect to each LED luminary's average optical power. Despite the relaxation of the rank constraint, the CCCP combined with the SDR approach achieves virtually the same performance as the CCCP approach. This reveals that the solution to (<ref>) can be either very near or rank-one. We also notice that the superiority of these two designs over the ZF one vanishes as $\overline{P}_n^s$ increases. As such, at the high optical power regime, one should employ the ZF design for the sake of low complexity while still being able to offer comparable performance to the optimal approaches.
In addition to the three described precoding algorithms, we consider a ZF precoder selection scheme, which simply chooses the best precoder from a number of random ZF ones. The advantage of this scheme is that it does not involve any optimization procedure, thus further reducing the design complexity. The use of this selection scheme seems reasonable when $\overline{P}_n^s$ is high (e.g., $35$ dBm or higher) since its performance losses concerning the CCCP (and CCCP combined with SDR) are only 2.6%, 3.9%, and 5.2% when the numbers of random samples are 100,000, 10,000, and 1,000, respectively.

SEE versus average optical power $\overline{P}_n^s$ of each LED luminary for different $P_{\text{DC, circuitry}}$.

Figure <ref> illustrates the impact of the circuitry power consumption on the SEE performance. In this simulation, results of the CCCP combined with the SDR approach are used. As can also be seen from the previous figure, there exists optimal $\overline{P}^s_n$ where the SEE achieves its maximum value. This is because, at the low $\overline{P}^s_n$ region, the increase in the achievable secrecy sum-rate is dominant in improving the SEE. When $\overline{P}^s_n$ further increases, it becomes the dominant factor that reduces the SEE as the achievable secrecy sum-rate only logarithmically increases with $\overline{P}^s_n$.
Interestingly, it is observed that the optimal $\overline{P}^s_n$ shifts to larger value as $P_{\text{DC, circuitry}}$ increases. For example, the optimal $\overline{P}^s_n$ values are 26, 28, 29, and 31 dBm when $P_{\text{DC, circuitry}} = $ 2, 4, 8, and 16 Watts, respectively.

SEE versus minimum threshold of secrecy rate.

Finally, assume that all users require the same minimum threshold of secrecy rate (i.e., $\lambda_k$'s $= \lambda$) we assess in Fig. <ref> the average SEE in accordance with $\lambda$ for different values of the average optical power per LED luminary. Within the examined value range, the average SEE increases with an increase in $\lambda$. The rate of increase, nonetheless, is different with $\overline{P}^s_n$. For instance, when $\lambda$ increases from 1 to 3 bits/s/Hz, the average SEE increases by 40.5%, 26.2%, and 10.5% for $\overline{P}^s_n = $ 25, 30, and 35 dBm, respectively. This result is intuitively obvious since at the high transmit optical power regime, increasing $\lambda$ (hence increasing the achievable secrecy sum-rate) does not lead to much EE improvement.

§ CONCLUSION

In this paper, we have addressed the problem of maximizing the SEE of multi-user VLC systems where users' messages are kept mutually confidential. Two different techniques based on the CCCP and SDR were presented to tackle the design problem. A sub-optimal low complexity design using ZF precoding was also investigated. Simulation results showed that although the CCCP combined with SDR had the highest theoretical worst-case complexity, it, in fact, achieves lower complexity in practice compared with the CCCP approach while providing the same performance. We also observed that the sub-optimal ZF design performance approached those of the CCCP and CCCP combined with SDR at the high transmit optical power regime. This motivated a random ZF scheme, which offered reasonable performance yet did not involve any optimization procedure.

In this appendix, we prove the lower bound of secrecy rate of MISO VLC system with confidential messages presented in (<ref>).

For demodulation, the DC term in (<ref>) is filtered out. The remaining AC term can then be expressed as

\begin{align}
\overline{y}_k = \mathbf{h}_k^T \mathbf{w}_k d_k + \mathbf{h}_k\sum_{i=1,\ i\neq k}^{K} \mathbf{w}_i d_i + \mathbf{h}_k \mathbf{I}^{\text{DC}} + \overline{n}_k,
\end{align}

where $\overline{n}_k = \frac{n_k}{\gamma\eta}$. According to [10], a lower bound on the secrecy capacity of the $k$-th user in an MISO VLC wiretap system can be derived as follows:

\begin{equation}
\label{eqn:Ck1}
\begin{split}
C_k \geq \left[h\left(\overline{y}_k\right)-h\left(\overline{y}_k|d_k\right)-h\left(\hat{\mathbf{y}}_k|\hat{\mathbf{d}}_k\right)+h\left(\hat{\mathbf{y}}_k|d_k;\hat{\mathbf{d}}_k\right)\right]^+  ,
\end{split}
\end{equation}

where $\hat{\mathbf{y}}_k=\begin{bmatrix}\overline{y}_1&\cdots&\overline{y}_{k-1}&\overline{y}_{k+1}& \cdots&\overline{y}_K\end{bmatrix}^T$ and $\hat{\mathbf{d}}_k=\begin{bmatrix}d_1 & \cdots & d_{k-1} & d_{k+1} & \cdots & d_k\end{bmatrix}^T$ are the vector $\overline{\mathbf{y}}$ and $\mathbf{d}$ having the $k$-th element removed.

In order to get a lower bound of $C_k$, we need to derive lower bounds of $h\left(\overline{y}_k\right)$ and $h\left(\hat{\mathbf{y}}_k|d_k;\hat{\mathbf{d}}_k\right)$, and  upper bounds of $h\left(\overline{y}_k|d_k\right)$ and $h\left(\hat{\mathbf{y}}_k|\hat{\mathbf{d}}_k\right)$.

Firstly, using the entropy power inequality [45], we obtain the lower bound of $h(\overline{y}_k)$ as

\begin{align}
\label{eqn:y_k}
h(\overline{y}_k) & = h\left(\sum_{i=1}^K \mathbf{h}_k^T \mathbf{w}_i d_i + \overline{n}_k\right)  \geq \frac{1}{2} \log_2 \left(\sum_{i=1}^{K} 2^{2h\left(\mathbf{h}_k^T \mathbf{w}_i d_i\right)} + 2^{2h(\overline{n}_k)}\right) \nonumber \\
& = \frac{1}{2} \log_2 \left(\sum_{i=1}^K \left(\mathbf{h}_k^T \mathbf{w}_i\right)^2 2^{2h_{d}} + 2\pi e \overline{\sigma}_k^2\right) \nonumber \\
& = \frac{1}{2} \log_2 \left(\sum_{i=1}^K \left(\mathbf{h}_k^T \mathbf{w}_i\right)^2 \frac{2^{2h_{d}}}{2\pi e \overline{\sigma}_k^2} + 1\right) + \frac{1}{2} \log_2 \left(2\pi e \overline{\sigma}_k^2\right).
\end{align}

Next, an upper bound of $h(\overline{y}_k|d_k)$ can be given by

\begin{align}
\label{eqn:y_k|d_k}
h\left(\overline{y}_k|d_k\right) & = h\left(\sum_{i=1,i\neq k}^K \mathbf{h}_k^T \mathbf{w}_i d_i + \overline{n}_k\right)
\stackrel{(a)}{\leq} \frac{1}{2} \log_2 \left( 2\pi e \left( \sum_{i=1,i\neq k}^{K} \left(\mathbf{h}_k^T \mathbf{w}_i\right)^2 \sigma^2_d + \overline{\sigma}_k^2 \right) \right) \nonumber \\
& = \frac{1}{2} \log_2 \left(\sum_{i=1,i\neq k}^K \left(\mathbf{h}_k^T \mathbf{w}_i\right)^2  \frac{2\pi e \sigma_d^2}{2\pi e \overline{\sigma}_k^2} + 1\right) + \frac{1}{2} \log_2 \left(2\pi e \overline{\sigma}_k^2\right).
\end{align}

where $(a)$ is due to the well-known result that $h(X) \leq \frac{1}{2}\log_2 \left(2\pi e \sigma^2_X\right)$ with $\sigma_X^2$ being the variance of the scalar random variable $X$.

Similarly, using the inequality $h(\mathbf{X}) \leq \log_2 \left(\left(2\pi e\right)^n \det\left(\mathbf{K}\right)\right)$ with $\mathbf{X} = \begin{bmatrix}x_1 & x_2 & \cdots & x_n\end{bmatrix}^T$ and $\mathbf{K}=\mathbb{E}\left[\mathbf{X}\mathbf{X}^T\right]$ being the covariance matrix of $\mathbf{X}$, we have

\begin{align}
\label{eqn:y_k|d_k2}
& = h\left(\hat{\mathbf{H}}_k \mathbf{w}_k d_k + \hat{\mathbf{n}}_k\right) \nonumber \\
& \leq \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \det\left(\mathbb{E}\left[\left(\hat{\mathbf{H}}_k \mathbf{w}_k d_k + \hat{\mathbf{n}}_k\right)\left(\hat{\mathbf{H}}_k \mathbf{w}_k d_k + \hat{\mathbf{n}}_k\right)^T\right] \right)\right) \nonumber \\
& = \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \det\left(\mathbb{E}\left[d_k^2\left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)\left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)^T + \hat{\mathbf{n}}_k \hat{\mathbf{n}}_k^T\right]\right)\right) \nonumber \\
& = \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \det\left(\sigma_d^2 \left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)\left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)^T + \text{diag}\left\{\hat{\pmb{\sigma}}^2_k\right\}\right)\right).
\end{align}

with $\hat{\mathbf{H}}_k = \begin{bmatrix}\mathbf{h}_1&\cdots&\mathbf{h}_{k-1}&\mathbf{h}_{k+1}&\cdots&\mathbf{h}_K\end{bmatrix}^T$, $\hat{\mathbf{n}}_k = \begin{bmatrix}\overline{n}_1&\cdots&\overline{n}_{k-1}&\overline{n}_{k+1}&\cdots&\overline{n}_K\end{bmatrix}^T$, and $\hat{\pmb{\sigma}}^2_k = \begin{bmatrix}\overline{\sigma}_1 & \cdots & \overline{\sigma}_{k-1}& \overline{\sigma}_{k+1}&\cdots&\overline{\sigma}_K\end{bmatrix}^T$. Then, using the matrix determinant lemma, i.e., $\det\left(\mathbf{A}+\mathbf{u}\mathbf{v}^T\right) = \det(\mathbf{A})\left(1+\mathbf{v}^T\mathbf{A}^{-1}\mathbf{u}\right)$,  for some invertible square matrix $\mathbf{A}$ and vectors $\mathbf{u}$, $\mathbf{v}$, results in

\begin{align}
\label{eqn:y-k|u-k}
& \leq \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \text{det}\left(\text{diag}\left\{\hat{\pmb{\sigma}}^2_k\right\}\right)\right.
\left.\left(1 + \sigma_d^2 \left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)^T \text{diag}\left\{\hat{\pmb{\sigma}}^{-2}_k\right\} \left(\hat{\mathbf{H}}_k \mathbf{w}_k\right)\right) \right) \nonumber \\
% & = \frac{1}{2} \log_2\left((2\pi e)^{K-1} \prod_{i=1,i \neq k}^K \sigma_i^2 \left(1 + \sum_{i=1,i\neq K}^K \left(\mathbf{h}_i^T \mathbf{w}_k\right)^2 \frac{\sigma_u^2}{\sigma_i^2}\right)\right) \\
& = \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \prod_{i=1,i \neq k}^K \overline{\sigma}_i^2\right) + \frac{1}{2} \log_2 \left(1 + \sum_{i=1,i\neq k}^K \left(\mathbf{h}_i^T \mathbf{w}_k\right)^2 \frac{\sigma_d^2}{\overline{\sigma}_i^2}\right).
\end{align}

Finally, $h\left(\hat{\mathbf{y}}_k|d_k,\hat{\mathbf{d}}_k\right)$ is straightforwardly given by

\begin{equation}
\label{eqn:y-k|u-k,d_k}
\begin{split}
h\left(\hat{\mathbf{y}}_k|d_k,\hat{\mathbf{d}}_k\right) & = h(\hat{\mathbf{n}}_k)  = \frac{1}{2} \log_2 \left((2\pi e)^{K-1} \prod_{i=1,i \neq k}^K \overline{\sigma}_i^2\right).
\end{split}
\end{equation}

By combining (<ref>), (<ref>), (<ref>), and (<ref>) and denoting $a_k = \frac{2^{2h_{d}}}{2\pi e \overline{\sigma}_k^2}$ and $b_k = \frac{\sigma_d^2}{\overline{\sigma}_k^2}$, we obtain the expression in (<ref>).

[1]
M. A. Arfaoui, M. D. Soltani, I. Tavakkolnia, A. Ghrayeb, M. Safari, C. Assi,
and H. Haas, “Physical layer security for visible light communication
systems: A survey,” IEEE Commun. Surveys Tuts., 2020.

[2]
Y. Zou, J. Zhu, X. Wang, and L. Hanzo, “A survey on wireless security:
Technical challenges, recent advances, and future trends,” Proceedings
of the IEEE, vol. 104, no. 9, pp. 1727–1765, 2016.

[3]
A. Mukherjee, S. A. A. Fakoorian, J. Huang, and A. L. Swindlehurst,
“Principles of physical layer security in multiuser wireless networks: A
survey,” IEEE Commun. Surveys Tuts., vol. 16, no. 3, pp. 1550–1573,

[4]
J.-Y. Wang, C. Liu, J.-B. Wang, Y. Wu, M. Lin, and J. Cheng, “Physical-layer
security for indoor visible light communications: Secrecy capacity
analysis,” IEEE Trans. Commun., vol. 66, no. 12, pp. 6423–6436,

[5]
X. Zhao, H. Chen, and J. Sun, “On physical-layer security in multiuser visible
light communication systems with non-orthogonal multiple access,” IEEE
Access, vol. 6, pp. 34 004–34 017, 2018.

[6]
A. Mostafa and L. Lampe, “Physical-layer security for miso visible light
communication channels,” IEEE J. Sel. Areas Commun., vol. 33, no. 9,
pp. 1806–1818, 2015.

[7]
——, “Optimal and robust beamforming for secure transmission in miso
visible-light communication links,” IEEE Trans. Signal Process.,
vol. 64, no. 24, pp. 6501–6516, 2016.

[8]
S. Ma, Z.-L. Dong, H. Li, Z. Lu, and S. Li, “Optimal and robust secure
beamformer for indoor miso visible light communication,” J. Lightw.
Technol., vol. 34, no. 21, pp. 4988–4998, 2016.

[9]
T. V. Pham and A. T. Pham, “Secrecy sum-rate of multi-user miso visible light
communication systems with confidential messages,” Optik, vol. 151,
pp. 65–76, 2017.

[10]
M. A. Arfaoui, A. Ghrayeb, and C. M. Assi, “Secrecy performance of multi-user
miso vlc broadcast channels with confidential messages,” IEEE Trans.
Wireless Commun., vol. 17, no. 11, pp. 7789–7800, 2018.

[11]
S. Cho, G. Chen, and J. P. Coon, “Zero-forcing beamforming for active
and passive eavesdropper mitigation in visible light communication systems,”
IEEE Transactions on Information Forensics and Security, vol. 16, pp.
1495–1505, 2021.

[12]
G. Pan, J. Ye, and Z. Ding, “On secure vlc systems with spatially random
terminals,” IEEE Communications Letters, vol. 21, no. 3, pp.
492–495, 2017.

[13]
S. Cho, G. Chen, and J. P. Coon, “Physical layer security in visible
light communication systems with randomly located colluding eavesdroppers,”
IEEE Wireless Communications Letters, vol. 7, no. 5, pp. 768–771,

[14]
——, “Securing visible light communication systems by beamforming in the
presence of randomly distributed eavesdroppers,” IEEE Transactions on
Wireless Communications, vol. 17, no. 5, pp. 2918–2931, 2018.

[15]
T. V. Pham, T. Hayashi, and A. T. Pham, “Artificial-noise-aided
precoding design for multi-user visible light communication channels,”
IEEE Access, vol. 7, pp. 3767–3777, 2019.

[16]
F. Wang, C. Liu, Q. Wang, J. Zhang, R. Zhang, L. Yang, and
L. Hanzo, “Optical jamming enhances the secrecy performance of the
generalized space-shift-keying-aided visible-light downlink,” IEEE
Transactions on Communications, vol. 66, no. 9, pp. 4087–4102, 2018.

[17]
M. A. Arfaoui, H. Zaid, Z. Rezki, A. Ghrayeb, A. Chaaban, and
M. Alouini, “Artificial noise-based beamforming for the miso vlc wiretap
channel,” IEEE Transactions on Communications, vol. 67, no. 4, pp.
2866–2879, 2019.

[18]
S. Cho, G. Chen, and J. P. Coon, “Enhancement of physical layer security
with simultaneous beamforming and jamming for visible light communication
systems,” IEEE Transactions on Information Forensics and Security,
vol. 14, no. 10, pp. 2633–2648, 2019.

[19]
T. V. Pham and A. T. Pham, “Energy efficient artificial noise-aided
precoding designs for secured visible light communication systems,”
IEEE Transactions on Wireless Communications, vol. 20, no. 1, pp.
653–666, 2021.

[20]
L. Yin and H. Haas, “Physical-layer security in multiuser visible light
communication networks,” IEEE J. Sel. Areas Commun., vol. 36, no. 1,
pp. 162–174, 2017.

[21]
M. Kashef, M. Ismail, M. Abdallah, K. A. Qaraqe, and E. Serpedin, “Energy
efficient resource allocation for mixed rf/vlc heterogeneous wireless
networks,” IEEE J. Sel. Areas Commun., vol. 34, no. 4, pp. 883–893,

[22]
A. Khreishah, S. Shao, A. Gharaibeh, M. Ayyash, H. Elgala, and
N. Ansari, “A hybrid rf-vlc system for energy efficient wireless access,”
IEEE Transactions on Green Communications and Networking, vol. 2,
no. 4, pp. 932–944, 2018.

[23]
H. Zhang, N. Liu, K. Long, J. Cheng, V. C. M. Leung, and L. Hanzo,
“Energy efficient subchannel and power allocation for software-defined
heterogeneous vlc and rf networks,” IEEE Journal on Selected Areas in
Communications, vol. 36, no. 3, pp. 658–670, 2018.

[24]
Y. Hsiao, Y. Wu, and C. Lin, “Energy-efficient beamforming design for
mu-miso mixed rf/vlc heterogeneous wireless networks,” IEEE
Transactions on Signal Processing, vol. 67, no. 14, pp. 3770–3784, 2019.

[25]
S. Aboagye, A. Ibrahim, T. M. N. Ngatched, and O. A. Dobre, “Vlc in
future heterogeneous networks: Energy^^e2^^80^^93 and
spectral^^e2^^80^^93efficiency optimization,” in ICC 2020 - 2020 IEEE
International Conference on Communications (ICC), 2020, pp. 1–7.

[26]
S. Aboagye, A. Ibrahim, T. M. N. Ngatched, A. R. Ndjiongue, and O. A.
Dobre, “Design of energy efficient hybrid vlc/rf/plc communication system
for indoor networks,” IEEE Wireless Communications Letters, vol. 9,
no. 2, pp. 143–147, 2020.

[27]
S. Ma, T. Zhang, S. Lu, H. Li, Z. Wu, and S. Li, “Energy
efficiency of siso and miso in visible light communication systems,”
Journal of Lightwave Technology, vol. 36, no. 12, pp. 2499–2509,

[28]
L. An, H. Shen, J. Wang, Y. Zeng, and R. Ran, “Energy efficiency
optimization for mimo visible light communication systems,” IEEE
Wireless Communications Letters, vol. 9, no. 4, pp. 452–456, 2020.

[29]
X. Liu, Y. Wang, F. Zhou, S. Ma, R. Q. Hu, and D. W. K. Ng, “Beamforming
design for secure miso visible light communication networks with slipt,”
IEEE Trans. Commun, 2020.

[30]
S. T. Duong, T. V. Pham, C. T. Nguyen, and A. T. Pham,
“Energy-efficient precoding for multi-user visible light communication with
confidential messages,” 2021 IEEE 93rd Vehicular Technology
Conference, 2021. [Online]. Available:

[31]
A. L. Yuille and A. Rangarajan, “The concave-convex procedure (cccp),”
Neural Comput., vol. 15, no. 4, pp. 915–936, 2003.

[32]
T. Lipp and S. Boyd, “Variations and extension of the convex–concave
procedure,” Optimization and Engineering, vol. 17, no. 2, pp.
263–287, 2016.

[33]
Y. Nesterov, “Semidefinite relaxation and nonconvex quadratic optimization,”
Optimization methods and software, vol. 9, no. 1-3, pp. 141–160,

[34]
Z. Luo, W. Ma, A. M. So, Y. Ye, and S. Zhang, “Semidefinite
relaxation of quadratic optimization problems,” IEEE Signal Processing
Magazine, vol. 27, no. 3, pp. 20–34, 2010.

[35]
T. Komine and M. Nakagawa, “Fundamental analysis for visible-light
communication system using led lights,” IEEE transactions on Consumer
Electronics, vol. 50, no. 1, pp. 100–107, 2004.

[36]
L. Zeng, D. C. O'Brien, H. Le Minh, G. E. Faulkner, K. Lee, D. Jung, Y. Oh, and
E. T. Won, “High data rate multiple input multiple output (mimo) optical
wireless communications using white led lighting,” IEEE Journal on
Selected Areas in Communications, vol. 27, no. 9, pp. 1654–1662, 2009.

[37]
W. Dinkelbach, “On nonlinear fractional programming,” Management
science, vol. 13, no. 7, pp. 492–498, 1967.

[38]
M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming, version 2.1,” <http://cvxr.com/cvx>, Mar. 2014.

[39]
R. Horst and N. V. Thoai, “Dc programming: Overview,” Journal of
Optimization Theory and Applications, vol. 103, pp. 1–43, 1999.

[40]
R. Feng, M. Dai, H. Wang, B. Chen, and X. Lin, “Linear precoding for multiuser
visible-light communication with field-of-view diversity,” IEEE
Photonics Journal, vol. 8, no. 2, pp. 1–8, 2016.

[41]
S. Ma, R. Yang, H. Li, Z.-L. Dong, H. Gu, and S. Li, “Achievable rate with
closed-form for siso channel and broadcast channel in visible light
communication networks,” Journal of Lightwave Technology, vol. 35,
no. 14, pp. 2778–2787, 2017.

[42]
C. Du, F. Zhang, S. Ma, Y. Tang, H. Li, H. Wang, and S. Li, “Secure
transmission for downlink noma visible light communication networks,”
IEEE Access, vol. 7, pp. 65 332–65 341, 2019.

[43]
M. A. Arfaoui, H. Zaid, Z. Rezki, A. Ghrayeb, A. Chaaban, and M.-S. Alouini,
“Artificial noise-based beamforming for the miso vlc wiretap channel,”
IEEE Transactions on Communications, vol. 67, no. 4, pp. 2866–2879,

[44]
Y. Nesterov and A. Nemirovskii, Interior-point polynomial algorithms in
convex programming.1em plus 0.5em minus 0.4emSIAM, 1994.

[45]
C. E. Shannon, “A mathematical theory of communication,” Bell Labs
Technical Journal, vol. 27, no. 3, pp. 379–423, 1948.

[1]
M. A. Arfaoui, M. D. Soltani, I. Tavakkolnia, A. Ghrayeb, M. Safari, C. Assi,
and H. Haas, “Physical layer security for visible light communication
systems: A survey,” IEEE Commun. Surveys Tuts., 2020.

[2]
Y. Zou, J. Zhu, X. Wang, and L. Hanzo, “A survey on wireless security:
Technical challenges, recent advances, and future trends,” Proceedings
of the IEEE, vol. 104, no. 9, pp. 1727–1765, 2016.

[3]
A. Mukherjee, S. A. A. Fakoorian, J. Huang, and A. L. Swindlehurst,
“Principles of physical layer security in multiuser wireless networks: A
survey,” IEEE Commun. Surveys Tuts., vol. 16, no. 3, pp. 1550–1573,

[4]
J.-Y. Wang, C. Liu, J.-B. Wang, Y. Wu, M. Lin, and J. Cheng, “Physical-layer
security for indoor visible light communications: Secrecy capacity
analysis,” IEEE Trans. Commun., vol. 66, no. 12, pp. 6423–6436,

[5]
X. Zhao, H. Chen, and J. Sun, “On physical-layer security in multiuser visible
light communication systems with non-orthogonal multiple access,” IEEE
Access, vol. 6, pp. 34 004–34 017, 2018.

[6]
A. Mostafa and L. Lampe, “Physical-layer security for miso visible light
communication channels,” IEEE J. Sel. Areas Commun., vol. 33, no. 9,
pp. 1806–1818, 2015.

[7]
——, “Optimal and robust beamforming for secure transmission in miso
visible-light communication links,” IEEE Trans. Signal Process.,
vol. 64, no. 24, pp. 6501–6516, 2016.

[8]
S. Ma, Z.-L. Dong, H. Li, Z. Lu, and S. Li, “Optimal and robust secure
beamformer for indoor miso visible light communication,” J. Lightw.
Technol., vol. 34, no. 21, pp. 4988–4998, 2016.

[9]
T. V. Pham and A. T. Pham, “Secrecy sum-rate of multi-user miso visible light
communication systems with confidential messages,” Optik, vol. 151,
pp. 65–76, 2017.

[10]
M. A. Arfaoui, A. Ghrayeb, and C. M. Assi, “Secrecy performance of multi-user
miso vlc broadcast channels with confidential messages,” IEEE Trans.
Wireless Commun., vol. 17, no. 11, pp. 7789–7800, 2018.

[11]
S. Cho, G. Chen, and J. P. Coon, “Zero-forcing beamforming for active
and passive eavesdropper mitigation in visible light communication systems,”
IEEE Transactions on Information Forensics and Security, vol. 16, pp.
1495–1505, 2021.

[12]
G. Pan, J. Ye, and Z. Ding, “On secure vlc systems with spatially random
terminals,” IEEE Communications Letters, vol. 21, no. 3, pp.
492–495, 2017.

[13]
S. Cho, G. Chen, and J. P. Coon, “Physical layer security in visible
light communication systems with randomly located colluding eavesdroppers,”
IEEE Wireless Communications Letters, vol. 7, no. 5, pp. 768–771,

[14]
——, “Securing visible light communication systems by beamforming in the
presence of randomly distributed eavesdroppers,” IEEE Transactions on
Wireless Communications, vol. 17, no. 5, pp. 2918–2931, 2018.

[15]
T. V. Pham, T. Hayashi, and A. T. Pham, “Artificial-noise-aided
precoding design for multi-user visible light communication channels,”
IEEE Access, vol. 7, pp. 3767–3777, 2019.

[16]
F. Wang, C. Liu, Q. Wang, J. Zhang, R. Zhang, L. Yang, and
L. Hanzo, “Optical jamming enhances the secrecy performance of the
generalized space-shift-keying-aided visible-light downlink,” IEEE
Transactions on Communications, vol. 66, no. 9, pp. 4087–4102, 2018.

[17]
M. A. Arfaoui, H. Zaid, Z. Rezki, A. Ghrayeb, A. Chaaban, and
M. Alouini, “Artificial noise-based beamforming for the miso vlc wiretap
channel,” IEEE Transactions on Communications, vol. 67, no. 4, pp.
2866–2879, 2019.

[18]
S. Cho, G. Chen, and J. P. Coon, “Enhancement of physical layer security
with simultaneous beamforming and jamming for visible light communication
systems,” IEEE Transactions on Information Forensics and Security,
vol. 14, no. 10, pp. 2633–2648, 2019.

[19]
T. V. Pham and A. T. Pham, “Energy efficient artificial noise-aided
precoding designs for secured visible light communication systems,”
IEEE Transactions on Wireless Communications, vol. 20, no. 1, pp.
653–666, 2021.

[20]
L. Yin and H. Haas, “Physical-layer security in multiuser visible light
communication networks,” IEEE J. Sel. Areas Commun., vol. 36, no. 1,
pp. 162–174, 2017.

[21]
M. Kashef, M. Ismail, M. Abdallah, K. A. Qaraqe, and E. Serpedin, “Energy
efficient resource allocation for mixed rf/vlc heterogeneous wireless
networks,” IEEE J. Sel. Areas Commun., vol. 34, no. 4, pp. 883–893,

[22]
A. Khreishah, S. Shao, A. Gharaibeh, M. Ayyash, H. Elgala, and
N. Ansari, “A hybrid rf-vlc system for energy efficient wireless access,”
IEEE Transactions on Green Communications and Networking, vol. 2,
no. 4, pp. 932–944, 2018.

[23]
H. Zhang, N. Liu, K. Long, J. Cheng, V. C. M. Leung, and L. Hanzo,
“Energy efficient subchannel and power allocation for software-defined
heterogeneous vlc and rf networks,” IEEE Journal on Selected Areas in
Communications, vol. 36, no. 3, pp. 658–670, 2018.

[24]
Y. Hsiao, Y. Wu, and C. Lin, “Energy-efficient beamforming design for
mu-miso mixed rf/vlc heterogeneous wireless networks,” IEEE
Transactions on Signal Processing, vol. 67, no. 14, pp. 3770–3784, 2019.

[25]
S. Aboagye, A. Ibrahim, T. M. N. Ngatched, and O. A. Dobre, “Vlc in
future heterogeneous networks: Energy^^e2^^80^^93 and
spectral^^e2^^80^^93efficiency optimization,” in ICC 2020 - 2020 IEEE
International Conference on Communications (ICC), 2020, pp. 1–7.

[26]
S. Aboagye, A. Ibrahim, T. M. N. Ngatched, A. R. Ndjiongue, and O. A.
Dobre, “Design of energy efficient hybrid vlc/rf/plc communication system
for indoor networks,” IEEE Wireless Communications Letters, vol. 9,
no. 2, pp. 143–147, 2020.

[27]
S. Ma, T. Zhang, S. Lu, H. Li, Z. Wu, and S. Li, “Energy
efficiency of siso and miso in visible light communication systems,”
Journal of Lightwave Technology, vol. 36, no. 12, pp. 2499–2509,

[28]
L. An, H. Shen, J. Wang, Y. Zeng, and R. Ran, “Energy efficiency
optimization for mimo visible light communication systems,” IEEE
Wireless Communications Letters, vol. 9, no. 4, pp. 452–456, 2020.

[29]
X. Liu, Y. Wang, F. Zhou, S. Ma, R. Q. Hu, and D. W. K. Ng, “Beamforming
design for secure miso visible light communication networks with slipt,”
IEEE Trans. Commun, 2020.

[30]
S. T. Duong, T. V. Pham, C. T. Nguyen, and A. T. Pham,
“Energy-efficient precoding for multi-user visible light communication with
confidential messages,” 2021 IEEE 93rd Vehicular Technology
Conference, 2021. [Online]. Available:

[31]
A. L. Yuille and A. Rangarajan, “The concave-convex procedure (cccp),”
Neural Comput., vol. 15, no. 4, pp. 915–936, 2003.

[32]
T. Lipp and S. Boyd, “Variations and extension of the convex–concave
procedure,” Optimization and Engineering, vol. 17, no. 2, pp.
263–287, 2016.

[33]
Y. Nesterov, “Semidefinite relaxation and nonconvex quadratic optimization,”
Optimization methods and software, vol. 9, no. 1-3, pp. 141–160,

[34]
Z. Luo, W. Ma, A. M. So, Y. Ye, and S. Zhang, “Semidefinite
relaxation of quadratic optimization problems,” IEEE Signal Processing
Magazine, vol. 27, no. 3, pp. 20–34, 2010.

[35]
T. Komine and M. Nakagawa, “Fundamental analysis for visible-light
communication system using led lights,” IEEE transactions on Consumer
Electronics, vol. 50, no. 1, pp. 100–107, 2004.

[36]
L. Zeng, D. C. O'Brien, H. Le Minh, G. E. Faulkner, K. Lee, D. Jung, Y. Oh, and
E. T. Won, “High data rate multiple input multiple output (mimo) optical
wireless communications using white led lighting,” IEEE Journal on
Selected Areas in Communications, vol. 27, no. 9, pp. 1654–1662, 2009.

[37]
W. Dinkelbach, “On nonlinear fractional programming,” Management
science, vol. 13, no. 7, pp. 492–498, 1967.

[38]
M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming, version 2.1,” <http://cvxr.com/cvx>, Mar. 2014.

[39]
R. Horst and N. V. Thoai, “Dc programming: Overview,” Journal of
Optimization Theory and Applications, vol. 103, pp. 1–43, 1999.

[40]
R. Feng, M. Dai, H. Wang, B. Chen, and X. Lin, “Linear precoding for multiuser
visible-light communication with field-of-view diversity,” IEEE
Photonics Journal, vol. 8, no. 2, pp. 1–8, 2016.

[41]
S. Ma, R. Yang, H. Li, Z.-L. Dong, H. Gu, and S. Li, “Achievable rate with
closed-form for siso channel and broadcast channel in visible light
communication networks,” Journal of Lightwave Technology, vol. 35,
no. 14, pp. 2778–2787, 2017.

[42]
C. Du, F. Zhang, S. Ma, Y. Tang, H. Li, H. Wang, and S. Li, “Secure
transmission for downlink noma visible light communication networks,”
IEEE Access, vol. 7, pp. 65 332–65 341, 2019.

[43]
M. A. Arfaoui, H. Zaid, Z. Rezki, A. Ghrayeb, A. Chaaban, and M.-S. Alouini,
“Artificial noise-based beamforming for the miso vlc wiretap channel,”
IEEE Transactions on Communications, vol. 67, no. 4, pp. 2866–2879,

[44]
Y. Nesterov and A. Nemirovskii, Interior-point polynomial algorithms in
convex programming.1em plus 0.5em minus 0.4emSIAM, 1994.

[45]
C. E. Shannon, “A mathematical theory of communication,” Bell Labs
Technical Journal, vol. 27, no. 3, pp. 379–423, 1948.