# The discrete mKdV equation revisited: a Riemann-Hilbert approach
00footnotetext:

Junyi Zhu, Xianguo Geng and Yonghui Kuang
School of Mathematics and Statistics, Zhengzhou University,
Zhengzhou, Henan 450001, People’s Republic of China Email<EMAIL_ADDRESS>
###### Abstract

We study the plus and minus type discrete mKdV equation. Some different
symmetry conditions associated with two Lax pairs are introduced to derive the
matrix Riemann-Hilbert problem with zero. By virtue of regularization of the
Riemann-Hilbert problem, we obtain the complex and real solution to the plus
type discrete mKdV equation respectively. Under the gauge transformation
between the plus and minus type, the solutions of minus type can be obtained
in terms of the given plus ones.

##  1 Introduction

The discrete mKdV (dmKdV) equation [3]

$u_{t}(n,t)=\left(1\pm u^{2}(n,t)\right)[u(n+1,t)-u(n-1,t)],$ (1.1)

is an integrable equation in mathematical physics, and it is an important
member of the discrete Ablowitz-Ladik equations [1, 2, 3, 4]. For specific
purpose, We call equation (1.1) the plus and minus type dmKdV. In this paper,
we study the plus type dmKdV equation with the help of the Riemann-Hilbert
(RH) method following [5], then the solutions of the minus type can be
obtained by virtue of a gauge transformation.

The plus type dmKdV equation (1.1) admits the following Lax pair formulation
[6]:

$\psi(n+1,t)=\gamma_{n}(I+Q_{n})Z\psi(n,t),\quad\psi_{t}(n,t)=(k\sigma_{3}+\tilde{Q}_{n})\psi(n,t),$
(1.2)

where $I$ is the identity matrix, $\gamma_{n}=\sqrt{\det(I+Q_{n})}^{-1}$, and
the matrices $Q_{n},Z,\tilde{Q}_{n}$ take the form

$\displaystyle Q_{n}=\left(\begin{matrix}0&u(n,t)\\\
-u(n,t)&0\end{matrix}\right),\quad Z=\left(\begin{matrix}z&0\\\
0&1\end{matrix}\right),$ (1.3)
$\displaystyle\tilde{Q}_{n}=Q_{n}+Z^{-1}Q_{n-1}Z,\quad
k=\frac{1}{2}(z-z^{-1}),$

with $z$ is a spectral parameter. We note that the Lax pair formulation (1.2)
can be rewritten as

$\psi(n+1,t)=(I+Q_{n})Z\psi(n,t),\quad\psi_{t}(n,t)=(k\sigma_{3}+\tilde{Q}_{n}-Q_{n}Q_{n-1})\psi(n,t),$
(1.4)

which are the ones in [6].

The plus type dmKdV equation (1.1) admits another Lax pair formulation [3]

$\varphi(n+1,t)=\gamma_{n}(E+Q_{n})\varphi(n,t),\quad\varphi_{t}(n,t)=(\omega\sigma_{3}+\hat{Q}_{n})\varphi(n,t),$
(1.5)

where $Q_{n}$ is defined as (1.3) and

$\displaystyle E=\left(\begin{matrix}z&0\\\
0&z^{-1}\end{matrix}\right),\quad\gamma_{n}=\frac{1}{\sqrt{\det(E+Q_{n})}},$
(1.6)
$\displaystyle\hat{Q}_{n}=EQ_{n}+Q_{n-1}E,\quad\omega=\frac{1}{2}(z^{2}-z^{-2}).$

We note that the Lax pair (1.5) can be rewritten in a similar form as (1.4)
which are the ones as in [3].

It is known that self-dual network can also be reduced to the discrete
analogue of the mKdV equation [7]. we note that there are many other
differential-difference equations which can be transformed into the dmKdV
equation [8, 9, 10, 11, 12, 13, 14, 15]. The dmKdV equation has widely
applications in the fields as plasma physics, electromagnetic waves in
ferromagnetic, antiferromagnetic or dielectric systems, and can be solved by
the method of inverse scattering transform, Hirota bilinear, Algebro-geometric
approach and others [3, 16, 17, 18, 19, 20, 6, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32].

In this paper, we firstly consider the Lax pair (1.2) and introduce a special
symmetry condition, which imply that $u(n,t)$ in (1.1) can be extended to
complex value. For simplicity, we suppose that $u(n,t)\neq\pm i$ and
$\gamma_{n}$ is chosen as one of branches. It is noted that the complex
solution in this paper is different from the complexiton solution introduced
by W.X. Ma [33, 34, 35, 36], in which the complexiton solutions are obtained
in the sense of complex eigenvalues, and are still real. Next we consider the
linear system (1.5) under the usual symmetry condition which confine the
potential $u(n,t)$ to be real. We note that to obtain the real solution of the
plus type dmKdV equation, one needs to introduce some constraint condition.
For one soliton solution as an example, we assume that $(2n+1)\eta+2\cosh
2\xi\sin 2\eta t=m\pi$, where $m=0,\pm 1,\pm 2,\cdots$ and the discrete
spectrum for $N=1$ is defined as $z_{1}=e^{\xi+i\eta}$.

The organization of this paper is as follows. In section 2, we derive the
complex solution of the plus type dmKdV equation by virtue of RH problem
associated with linear system (1.2). In section 3, we derive the real solution
of the plus type dmKdV equation by virtue of RH problem associated with linear
system (1.5). In section 4, we study the gauge transformation between the plus
and minus type dmKdV equation, from which the solution of minus type can be
obtained in terms of the given plus one.

##  2 Complex solution of the dmKdV equation

###  2.1 The spectral analysis

For the sake of convenience, we write the spectral equation (1.2) in terms of
the matrix

$J(n)=\psi(n)Z^{-n}e^{-k\sigma_{3}t}.$

Hence, the dmKdV equation allows the Lax representation:

$J(n+1)=\gamma_{n}(I+Q_{n})ZJ(n)Z^{-1},$ (2.7)

and

$J_{t}(n)=k[\sigma_{3},J(n)]+\tilde{Q}_{n}J(n).$ (2.8)

Here and after we suppress the variables dependence for simplicity of
notation.

Now we introduce matrix Jost functions $J_{\pm}(n,z)$ of the spectral equation
(2.7) obeying the asymptotic conditions

$J_{\pm}(n,z)\rightarrow I,\quad n\rightarrow\pm\infty.$ (2.9)

Then there exists the scattering matrix $S(z)$ admitting

$J_{-}(n,z)=J_{+}(n,z)Z^{n}S(z)Z^{-n},\quad
S(z)=\left(\begin{matrix}a_{+}(z)&-b_{-}(z)\\\
b_{+}(z)&a_{-}(z)\end{matrix}\right).$ (2.10)

Here we assume that the Jost functions and the scattering matrix satisfy the
symmetry condition

$J_{\pm}^{T}(n,z^{-})=J_{\pm}^{-1}(n,z),\quad\det J_{\pm}(n,z)=1,$ (2.11)

and

$S^{T}(n,z^{-})=S^{-1}(n,z),\quad\det S(n,z)=1,\quad z^{-}=z^{-1}.$ (2.12)

In the following, we consider the asymptotic behavior of the solution
$J(n,z)$. To this end, we first let

$J(n,z)=J^{(0)}(n)+z^{-1}J^{(1)}(n)+O(z^{-2}),\quad z\rightarrow\infty,$
(2.13)

and substitute it into the spectral equation (2.7). This yields

$J_{11}^{(0)}(n+1)=\gamma_{n}J_{11}^{(0)}(n),\quad
J_{22}^{(0)}(n+1)=\gamma_{n}^{-1}J_{22}^{(0)}(n),$ (2.14)

and $J_{12}^{(0)}(n)=0$,

$J_{12}^{(1)}(n)=-u(n)J_{22}^{(0)}(n),\quad
J_{21}^{(0)}(n+1)=-u(n)J_{11}^{(0)}(n+1).$ (2.15)

Next, according to the symmetry (2.11), we let

$J^{-1}(n,z)=\tilde{J}_{(0)}+z\tilde{J}_{(1)}+O(z^{2}),\quad z\rightarrow 0,$

and find similarly $\tilde{J}_{(0)21}(n)=0$,

$\displaystyle\tilde{J}_{(0)11}(n+1)=\gamma_{n}\tilde{J}_{(0)11}(n),\quad\tilde{J}_{(0)22}(n+1)=\gamma_{n}^{-1}\tilde{J}_{(0)22}(n),$
(2.16)
$\displaystyle\tilde{J}_{(0)12}(n+1)=-u(n)\tilde{J}_{(0)11}(n+1),\quad\tilde{J}_{(1)21}(n)=-u(n)\tilde{J}_{(0)22}(n).$

We will now discuss the analytic of the Jost solutions. The spectral equation
(2.7), as a iterative relation, can be written as

$J(n,z)=\nu_{+}(n)\lim\limits_{N\rightarrow\infty}\prod\limits_{l=n}^{N}(Z^{-1}(I-Q_{l}))J(N+1,z)Z^{N-n+1},$
(2.17)

here and after we introduce two new functions $\nu_{\pm}(n)$ as following

$\nu_{+}(n)=\prod\limits_{l=n}^{\infty}\gamma_{l},\quad\nu_{-}(n)=\prod\limits_{l=-\infty}^{n-1}\gamma_{l}.$
(2.18)

We note that the first column of the matrix equation (2.17) involves two
positive power series in $z$, while the second column involves two negative
power series in $z$. Thus the first column $J_{+}^{[1]}(n,z)$ of the Jost
function $J_{+}$ is analytical for $|z|<1$, denoted by ${\mathbb{C}}_{I}$, and
the second column $J_{+}^{[2]}(n,z)$ is analytical for $|z|>1$ or
(${\mathbb{C}}_{O}$). By the same way one can show that the column
$J_{-}^{[1]}(n,z)$ is analytical for $|z|>1$ or (${\mathbb{C}}_{O}$). We
introduce a matrix function

$\Phi_{+}(n,z)=\left(J_{-}^{[1]},J_{+}^{[2]}\right)$

which is analytical in ${\mathbb{C}}_{O}$ and solves the spectral equation
(2.7).

It follows from the symmetry condition (2.11) that the rows
$(J_{-})_{[1]}^{-1}$ and $(J_{+})_{[2]}^{-1}$ are analytical in
${\mathbb{C}}_{I}$. Thus the matrix function

$\Phi_{-}^{-1}(n,z)=\left(\begin{array}[]{c}(J_{-})_{[1]}^{-1}\\\
(J_{+})_{[2]}^{-1}\end{array}\right)$

is analytical in ${\mathbb{C}}_{I}$ and solves the adjoint spectral problem of
(2.7).

By virtue of the definition (2.10) of the scattering matrix , we find

$\Phi_{+}(n,z)=J_{\pm}Z^{n}S_{\pm}Z^{-n},$ (2.19)

where

$S_{+}=\left(\begin{matrix}a_{+}&0\\\ b_{+}&1\end{matrix}\right),\quad
S_{-}=\left(\begin{matrix}1&b_{-}\\\ 0&a_{+}\end{matrix}\right).$

Hence, on use of (2.11), we obtain

$\det\Phi_{+}=\det J_{\pm}\det S_{\pm}=a_{+}(z).$ (2.20)

Following the same procedure as the one used for $\Phi_{+}$, one obtains

$\displaystyle\Phi_{-}^{-1}(n,z)=Z^{n}T_{\pm}Z^{-n}J_{\pm}^{-1},\quad\det\Phi_{-}^{-1}(n,z)=a_{-}(z),$
(2.21) $\displaystyle T_{+}=\left(\begin{matrix}a_{-}&b_{-}\\\
0&1\end{matrix}\right),\qquad T_{-}=\left(\begin{matrix}1&0\\\
b_{+}&a_{-}\end{matrix}\right).$

Asymptotic formulae for these sectionally analytic functions can be derived
from equations (2.13) to (2.16),

$\Phi_{+}(n,z)\rightarrow\Phi_{+}^{(0)}(n)=\left(\begin{matrix}\nu_{-}(n)&0\\\
-u(n-1)\nu_{-}(n)&\nu_{+}(n)\end{matrix}\right),\quad z\rightarrow\infty,$
(2.22)

and

$\Phi_{-}^{-1}(n,z)\rightarrow\Phi_{-(0)}^{-1}(n)=\left(\begin{matrix}\nu_{-}(n)&-u(n-1)\nu_{-}(n)\\\
0&\nu_{+}(n)\end{matrix}\right),\quad z\rightarrow 0.$ (2.23)

where $\nu_{\pm}(n)$ defined in (2.18). Indeed, to obtain equation (2.22), we
know, from the definition of $\Phi_{+}$ and the asymptotic expansion (2.13),
that

$\Phi_{+}^{(0)}(n)=\left(\begin{matrix}J_{-11}^{(0)}(n)&0\\\
J_{-21}^{(0)}(n)&J_{+22}^{(0)}(n)\end{matrix}\right),$

Iterating the relations (2.14) and (2.15), we find

$\displaystyle
J_{-11}^{(0)}(n)=\gamma_{n-1}J_{-11}^{(0)}(n-1)=\cdots=\nu_{-}(n),$
$\displaystyle J_{-21}^{(0)}(n)=-u(n-1)J_{-11}^{(0)}(n)=-u(n-1)\nu_{-}(n),$
$\displaystyle
J_{+22}^{(0)}(n)=\gamma_{n}J_{+22}^{(0)}(n+1)=\cdots=\nu_{+}(n),$

which give (2.22) in terms of the boundary condition (2.9). Equation (2.23)
can be obtained from (2.16) in a same way.

We note that the symmetry condition about these sectionally analytic functions
can be obtained from that of the Jost solutions as

$\Phi_{+}^{T}(n,z^{-})=\Phi_{-}^{-1}(n,z).$ (2.24)

In addition, equations (2.19) and (2.20) imply that $a_{+}(z)$ and $a_{-}(z)$
are analytical in the domain of ${\mathbb{C}}_{O}$ and ${\mathbb{C}}_{I}$
respectively. Furthermore they admit the following asymptotic behavior

$a_{+}(z)\rightarrow\nu,~{}z\rightarrow\infty;\quad
a_{-}(z)\rightarrow\nu,~{}z\rightarrow 0,$ (2.25)

where $\nu=\nu_{+}(n)\nu_{-}(n)=\prod_{l=-\infty}^{\infty}\gamma_{l}$.

It is noted that the potential $u(n)$ can be reconstructed by the analytic
functions. Indeed, from the first equation of (2.15), we find

$u=-\frac{J_{12}^{(1)}(n)}{J_{22}^{(0)}(n)}=-\lim\limits_{z\rightarrow\infty}\frac{(z\Phi_{+})_{12}}{(\Phi_{+})_{22}}=-\frac{\Phi_{+12}^{(1)}}{\Phi_{+22}^{(0)}},$
(2.26)

while the second equation of (2.15) is an identity.

###  2.2 RH problem and its regularization

Now we can introduce the RH problem

$\displaystyle\Phi_{-}^{-1}(n,z)\Phi_{+}(n,z)=Z^{n}G(z)Z^{-n},\quad|z|=1,$
(2.27) $\displaystyle
G(z)=T_{+}S_{+}=T_{-}S_{-}=\left(\begin{matrix}1&b_{-}(z)\\\
b_{+}(z)&1\end{matrix}\right).$

The normalization of the RH problem is given by (2.22) which is noncanonical.
Hence the dmKdV potential can be retrieved by virtue of the solution of RH
problem.

In order to obtain the soliton solutions of the dmKdV equation, we take
$G(z)=I$ and suppose $a_{+}(z)$ has simple zeros at
$z_{j}\in{\mathbb{C}_{O}},~{}j=1,\cdots,N$. From the symmetry (2.24), we know
that
$\det\Phi_{+}(z_{j})=0,~{}\det\Phi_{-}^{-1}(z^{-}_{l})=0,~{}j,l=1,\cdots,N$.
In this case, problem (2.27) is called the RH one with zeros which can be
solved by virtue of its regularization.

To obtain the relevant regular problem, we introduce a rational matrix
function

$\chi_{j}^{-1}=I+\frac{z_{j}-z^{-}_{j}}{z-z_{j}}P_{j},\quad
P_{j}=\frac{|y_{j}\rangle\langle\tilde{y}_{j}|}{\langle\tilde{y}_{j}|y_{j}\rangle},$

where the eigenvector $|y_{j}\rangle$ solves
$\Phi_{+}(n,z_{j})|y_{j}\rangle=0$. Since $\Phi_{+}(n,z_{j})$ admits the
linear system (2.7) and (2.8), then we have

$\displaystyle\Phi_{+}(n,t,z_{j})Z^{-1}(z_{j})|y_{j}\rangle(n+1,t)=0,$
$\displaystyle\Phi_{+}(n,t,z_{j})(|y_{j}\rangle_{t}-k_{j}\sigma_{3}|y_{j}\rangle)(n,t)=0,$

which imply that

$|y_{j}\rangle(n,t)=Z^{n}(z_{j})e^{k_{j}\sigma_{3}t}|y_{j}\rangle_{0},\quad
k_{j}=(z_{j}-z_{j}^{-})/2,$ (2.28)

where $|y_{j}\rangle_{0}$ is an arbitrary constant vector. In addition, one
finds that $\langle\tilde{y}_{j}|=|y_{j}\rangle^{T}$ satisfies
$\langle\tilde{y}_{j}|\Phi_{-}^{-1}(n,z_{l}^{-})=0$.

Therefore the product $\Phi_{+}(z)\chi_{j}^{-1}(z)$ is regular at the point
$z_{j}$ and $\chi_{l}(z)\Phi_{-}^{-1}(z)$ is regular at $z^{-}_{l}$, where

$\chi_{l}=I-\frac{z_{l}-z^{-}_{l}}{z-z^{-}_{l}}P_{l}.$ (2.29)

The regularization of all the other zeros is performed similarly and
eventually we obtain the following representation for the analytic solution

$\Phi_{\pm}=\phi_{\pm}\Gamma,\quad\Gamma=\chi_{N}\chi_{N-1}\cdots\chi_{1},$
(2.30)

where the holomorphic matrix functions $\phi_{\pm}$ solve the regular RH
problem

$\phi_{-}^{-1}(n,z)\phi_{+}(n,z)=I.$ (2.31)

We note that the soliton matrix $\Gamma$ can be decomposed into simple
fractions

$\Gamma=I-\sum\limits_{j,l=1}^{N}\frac{1}{z-z^{-}_{l}}|y_{j}\rangle(D^{-1})_{jl}\langle\tilde{y}_{l}|,\quad
D_{lj}=\frac{\langle\tilde{y}_{l}|y_{j}\rangle}{z_{j}-z_{l}^{-}}.$ (2.32)

In the following, we will establish the relationship between the solution of
dmKdV equation and the soliton matrix. Taking into account the asymptotic
formula (2.22) and the expression of $\Gamma$ (2.32), we choose
$\Phi_{+}^{(0)}=\phi_{+}$. Then, in view of (2.30), we find

$\Gamma(n,z)=I+z^{-1}\Gamma^{(1)}(n)+O(z^{-2}),\quad z\rightarrow\infty.$
(2.33)

and $\Phi_{+12}^{(1)}(n)=\nu_{-}(n)\Gamma_{12}^{(1)}(n)$. In addition, the
assumption $G(z)=I$ implies that $b_{\pm}(z)=0$ and then $a_{+}(z)a_{-}(z)=1$
in view of (2.12). From (2.25), we know that $\nu=\nu_{+}(n)\nu_{-}(n)=1$.
Hence the potential $u(n)$ can be rewritten as

$u(n)=-\frac{\nu_{-}(n)\Gamma_{12}^{(1)}(n)}{\nu_{+}(n)}=-\nu_{-}^{2}(n)\Gamma_{12}^{(1)}(n).$
(2.34)

Next we will establish the relationship between $\nu_{-}$ and $\Gamma$. Since
$G(z)=I$, the RH problem (2.27) reduces to $\Phi_{+}=\Phi_{-}$, from which we
can consider the asymptotic behavior of $\Phi_{+}$ near $z=0$. Indeed, the
asymptotic formulae (2.21) and (2.22) imply that

$\Phi_{+}\rightarrow\Phi_{-(0)}(n)=\left(\begin{matrix}\nu_{+}(n)&u(n-1)\nu_{-}(n)\\\
0&\nu_{-}(n)\end{matrix}\right),\quad z\rightarrow 0.$

Thus from (2.30) we obtain

$\displaystyle\Gamma(n,z)|_{z=0}$
$\displaystyle=(\Phi_{+}^{(0)})^{-1}(n)\Phi_{+}|_{z=0}$ (2.35)
$\displaystyle=\left(\begin{matrix}\nu_{+}^{2}(n)&u(n-1)\\\
u(n-1)&\nu_{-}^{2}(n-1)\end{matrix}\right),$

which implies that $\nu_{-}^{2}(n)=\Gamma_{22}(n+1,z=0)$. As a result, the
potential $u(n)$ takes the form

$u(n)=-\Gamma_{12}^{(1)}(n)\Gamma_{22}(n+1,z=0).$ (2.36)

###  2.3 Complex soliton solutions

In this section, we will derive the soliton solutions of the dmKdV equation
(1.1). To this end, we let

$z_{j}=e^{\xi_{j}+i\eta_{j}},~{}\xi_{j}>0,\quad|y_{j}\rangle_{0}=\left(\begin{array}[]{c}e^{a_{j}+i\alpha_{j}}\\\
1\end{array}\right).$

Hence the vector $|y_{j}\rangle,~{}(j=1,2,\cdots,N)$ take the form

$|y_{j}\rangle=e^{\frac{1}{2}(\theta_{j}(n)+i\phi_{j}(n))}\left(\begin{array}[]{c}e^{\frac{1}{2}(X_{j}(n,t)+i\varphi_{j}(n,t))}\\\
e^{-\frac{1}{2}(X_{j}(n,t)+i\varphi_{j}(n,t))}\end{array}\right),$ (2.37)

where

$\displaystyle X_{j}(n,t)=n\xi_{j}+2\sinh\xi_{j}\cos\eta_{j}t+a_{j},$ (2.38)
$\displaystyle\varphi_{j}(n,t)=n\eta_{j}+2\cosh\xi_{j}\sin\eta_{j}t+\alpha_{j},$

with $\theta_{j}(n)=X_{j}(n,0),\phi_{j}(n)=\varphi_{j}(n,0)$.

In particularly, for $N=1$, equation (2.32) reduces to

$\Gamma(n,z)=I-\frac{z_{1}-z_{1}^{-}}{z-z_{1}^{-}}\frac{|y_{1}\rangle\langle\tilde{y}_{1}|}{\langle\tilde{y}_{1}|y_{1}\rangle},$
(2.39)

from which we have the complex form of one-soliton solution to dmKdV equation

$u(n,t)=\frac{z_{1}^{2}-1}{2}{\rm
sech}\\{X_{1}(n+1,t)+i\varphi_{1}(n+1,t)\\}.$ (2.40)

For $N=2$, we find

$\displaystyle\Gamma^{(1)}(n,z)=$ $\displaystyle-\frac{1}{\det
D}\left\\{D_{22}|y_{1}\rangle\langle\tilde{y}_{1}|-D_{21}|y_{2}\rangle\langle\tilde{y}_{1}|\right.$
(2.41)
$\displaystyle\qquad\left.+D_{11}|y_{2}\rangle\langle\tilde{y}_{2}|-D_{12}|y_{1}\rangle\langle\tilde{y}_{2}|\right\\},$

and

$\displaystyle\Gamma(n,0)=I+$ $\displaystyle\frac{1}{\det
D}\left\\{z_{1}\left[D_{22}|1\rangle\langle\tilde{1}|B-D_{21}|2\rangle\langle\tilde{1}|B\right]\right.$
(2.42)
$\displaystyle\qquad\left.+z_{2}\left[D_{11}|2\rangle\langle\tilde{2}|B-D_{12}|1\rangle\langle\tilde{2}|B\right]\right\\},$

where $\det D$ is obtained according to the definition of (2.32) and (2.37) as

$\det
D=\Xi\Omega_{2}(n),\quad\Xi=\left(\prod\limits_{j,l=1}^{2}(z_{j}-z_{l}^{-})\right)^{-1}\frac{2e^{\theta_{1}+\theta_{2}+i(\phi_{1}+\phi_{2})}}{z_{1}z_{2}},$
(2.43)

and

$\displaystyle\Omega_{2}(n)=$
$\displaystyle(z_{1}-z_{2})^{2}\cosh\\{\vartheta_{1}(n,t)+\vartheta_{2}(n,t)\\}$
(2.44)
$\displaystyle+(z_{1}z_{2}-1)^{2}\cosh\\{\vartheta_{1}(n,t)-\vartheta_{2}(n,t)\\}-(z_{1}^{2}-1)(z_{2}^{2}-1),$

with

$\vartheta_{j}(n,t)=X_{j}(n,t)+i\varphi_{j}(n,t).$

In addition, from (2.37), (2.41) and (2.42), we know that

$\Gamma_{12}^{(1)}(n)=-\frac{V_{2}(n+1)}{\Omega_{2}(n)},\quad\Gamma_{22}(n,0)=\frac{\Omega_{2}(n)}{\Omega_{2}(n+1)},$
(2.45)

where

$\displaystyle V_{2}(n)=$
$\displaystyle(z_{2}-z_{1})(z_{1}z_{2}-1)\left[z_{1}(z_{2}^{2}-1)\cosh\\{\vartheta_{1}(n,t)\\}\right.$
(2.46)
$\displaystyle\quad\left.-z_{2}(z_{1}^{2}-1)\cosh\\{\vartheta_{2}(n,t)\\}\right].$

Hence the solution of the plus type dmKdV equation for $N=2$ can be given by

$u(n)=\frac{V_{2}(n+1)}{\Omega_{2}(n+1)}.$ (2.47)

It is noted that the solution $u(n)$ can also be derived through
$\Gamma_{12}(n+1,z=0)$ by (2.35). It is verified that the representations of
solution by $\Gamma_{12}(n+1,z=0)$ are same as the ones in (2.40) and (2.47).

##  3 Real solutions of the dmKdV equation

###  3.1 The spectral analysis

In the section, we consider the inear system (1.5) and assume that the
solution $u(n,t)$ is a real function. After the transformation

$J(n)=\varphi(n)E^{-n}e^{-\omega\sigma_{3}t},$

the dmKdV equation allows the Lax representation:

$J(n+1)=\gamma_{n}(E+Q_{n})J(n)E^{-1},$ (3.1)

and

$J_{t}(n)=\omega[\sigma_{3},J(n)]+\hat{Q}_{n}J(n).$ (3.2)

We assume that the function $J(n,z)$ admits the following symmetry conditions

$J^{\dagger}(n,\bar{z})=J^{-1}(n,z),\quad\sigma_{3}J(n,-z)\sigma_{3}=J(n,z),$
(3.3)

where $\bar{z}=(z^{*})^{-1}$ with $z^{*}$ denotes the complex conjugate of
$z$.

The Jost functions $J_{\pm}(n,z)$ and the scattering matrix $S(z)$ can be
introduced in the same way as in (2.9) and (2.10). It is readily verified that
the matrices $J_{\pm}(n,z)$ and $S(z)$ are unimodular, and satisfy the
symmetry conditions (3.3).

We note that similar considerations apply to the asymptotic behavior of the
Jost functions $J_{\pm}(n,z)$, one find

$\displaystyle J(n,z)=J^{(0)}(n)+z^{-1}J^{(1)}(n)+O(z^{-2}),\quad
z\rightarrow\infty,$ (3.4) $\displaystyle
J(n,z)=J_{(0)}(n)+zJ_{(1)}(n)+O(z^{2}),\quad z\rightarrow 0,$

where the diagonal matrices $J^{(0)}(n)$ and $J_{(0)}(n)$ admit the following
iterative relations

$J^{(0)}(n+1)=\left(\begin{matrix}\gamma_{n}&0\\\
0&\gamma_{n}^{-1}\end{matrix}\right)J^{(0)}(n),\quad
J_{(0)}(n+1)=\left(\begin{matrix}\gamma_{n}^{-1}&0\\\
0&\gamma_{n}\end{matrix}\right)J_{(0)}(n).$ (3.5)

In addition, the solution can be constructed by

$u(n)=-\frac{J_{12}^{(1)}(n)}{J_{22}^{(0)}(n)}.$ (3.6)

The analytical properties of the Jost functions $J_{\pm}(n,z)$ are the same as
the ones in complex section above, and can be used to define the same
sectionally holomorphic $\Phi_{+}(n,z)$ and $\Phi_{-}^{-1}(n,z)$ as (2.19) and
(2.21) with $Z$ replaced by $E$. Furthermore, we have the symmetry condition
about $\Phi_{\pm}(n,z)$

$\Phi_{+}^{\dagger}(n,\bar{z})=\Phi^{-1}(n,z),$ (3.7)

and the asymptotic behavior

$\Phi_{+}(n,z)\rightarrow\Phi_{+}^{(0)}(n)=\left(\begin{matrix}\nu_{-}(n)&0\\\
0&\nu_{+}(n)\end{matrix}\right),\quad z\rightarrow\infty,\\\ $ (3.8)

and

$\Phi_{-}^{-1}(n,z)\rightarrow\tilde{\Phi}_{-}^{(0)}(n)=\left(\begin{matrix}\nu_{-}(n)&0\\\
0&\nu_{+}(n)\end{matrix}\right),\quad z\rightarrow 0,$ (3.9)

where the real functions $\nu_{\pm}(n)$ are defined as in (2.18).

###  3.2 The regularization of the RH problem and the soliton solutions

The RH problem associated with $\Phi_{\pm}(n,z)$ can be constructed as in
(2.27), while the normalization of the RH problem is given by (3.8).

In order to obtain the real soliton solutions of the dmKdV equation, we take
$G(z)=I$ and suppose $a_{+}(z)=\det\Phi_{+}(n,z)$ has simple zeros at $\pm
z_{j}\in{\mathbb{C}_{O}},~{}j=1,\cdots,N$. From the symmetries (3.7), we know
that $\det\Phi_{+}(\pm
z_{j})=0,~{}\det\Phi_{-}^{-1}(\pm\bar{z}_{l})=0,~{}j,l=1,\cdots,N$.

For convenience, we introduce the notations

$k_{2j}=z_{j},\quad k_{2j-1}=-z_{j}.$

Then the soliton matrix can be written in the form

$\Gamma(n,z)=\chi_{2N}\chi_{2N-1}\cdots\chi_{2}\chi_{1},$ (3.10)

where

$\chi_{l}=I-\frac{k_{l}-\bar{k}_{l}}{z-\bar{k}_{l}}P_{l},\quad\chi_{j}^{-l}=I+\frac{k_{j}-\bar{k}_{j}}{z+k_{j}}P_{j},\quad
P_{j}=\frac{|j\rangle\langle j|}{\langle j|j\rangle},$ (3.11)

with the eigenvector $\langle j|=|j\rangle^{\dagger}$ and $|j\rangle$ solves
$\Phi_{+}(n,k_{j})|j\rangle=0$. In this case, one may find that
$|2j\rangle=\sigma_{3}|2j-1\rangle$ and $P_{2j}=\sigma_{3}P_{2j-1}\sigma_{3}$.
Hence the product $\Phi_{+}(z)\chi_{j}^{-1}(z)$ is regular at the point
$k_{j}$ and $\chi_{l}(z)\Phi_{-}^{-1}(z)$ is regular at $\bar{k}_{l}$.

Since $\Phi_{+}(n,z)$ solves the linear system (3.1) and (3.2), we know that
the eigenvector $|j\rangle$ takes the form

$|j\rangle(n,t)=E^{n}(k_{j})e^{\omega_{j}\sigma_{3}t}|j_{0}\rangle,\quad\omega_{j}=\omega(k_{j}).$
(3.12)

The regular RH problem can be derived similarly as (2.31) and (2.30), where
the soliton matrix $\Gamma$ has the following decomposition

$\Gamma=I-\sum\limits_{j,l=1}^{2N}\frac{1}{z-\bar{k}_{l}}|j\rangle(D^{-1})_{jl}\langle
l|,\quad D_{lj}=\frac{\langle l|j\rangle}{k_{j}-\bar{k}_{l}}.$ (3.13)

For $N=1$, we take $z_{1}=e^{\xi+i\eta}$, that is $k_{2}=z_{1},k_{1}=-z_{1}$,

$|2\rangle=\left(\begin{array}[]{c}e^{\theta+i\phi}\\\
e^{-(\theta+i\phi)}\end{array}\right),\quad|1\rangle=\sigma_{3}|2\rangle,$
(3.14)

where

$\theta=n\xi+\sinh 2\xi\cos 2\eta t+\alpha,\quad\phi=n\eta+\cosh 2\xi\sin
2\eta t+\beta.$ (3.15)

In this case, the soliton matrix takes the form

$\Gamma(n,t,z)=I-\frac{D_{-}}{z-\bar{z}_{1}}-\frac{D_{+}}{z+\bar{z}_{1}},$
(3.16)

where

$D_{-}=\frac{z_{1}-\bar{z}_{1}}{2}\left(\begin{matrix}\frac{e^{2\theta}}{z_{1}e^{-2\theta}+\bar{z}_{1}e^{2\theta}}&\frac{e^{2i\phi}}{z_{1}e^{-2\theta}+\bar{z}_{1}e^{2\theta}}\\\
\frac{e^{-2i\phi}}{z_{1}e^{2\theta}+\bar{z}_{1}e^{-2\theta}}&\frac{e^{-2\theta}}{z_{1}e^{2\theta}+\bar{z}_{1}e^{-2\theta}}\end{matrix}\right),\quad
D_{+}=-\sigma_{3}D_{-}\sigma_{3}.$ (3.17)

From (3.16) and (3.17), we have

$\Gamma(n,z=0)=\frac{z_{1}}{\bar{z}_{1}}\left(\begin{matrix}\frac{z_{1}e^{2\theta}+\bar{z}_{1}e^{-2\theta}}{z_{1}e^{-2\theta}+\bar{z}_{1}e^{2\theta}}&0\\\
0&\frac{z_{1}e^{-2\theta}+\bar{z}_{1}e^{2\theta}}{z_{1}e^{2\theta}+\bar{z}_{1}e^{-2\theta}}\end{matrix}\right),$
(3.18)

and

$\Gamma^{(1)}(n)=-(z_{1}^{2}-\bar{z}_{1}^{2})\left(\begin{matrix}0&\frac{e^{2i\phi}}{z_{1}e^{-2\theta}+\bar{z}_{1}e^{2\theta}}\\\
\frac{e^{-2i\phi}}{z_{1}e^{2\theta}+\bar{z}_{1}e^{-2\theta}}&0\end{matrix}\right),$
(3.19)

where $\Gamma^{(1)}(n)$ is defined by the asymptotic behavior

$\Gamma(n,z)=I+z^{-1}\Gamma^{(1)}(n)+O(z^{-2}),\quad z\rightarrow\infty.$
(3.20)

Next we will give the solution of dmKdV equation (1.1) for $N=1$. To this end,
we take the asymptotic behavior of the sectionally holomorphic $\Phi_{+}(n,z)$
as

$\Phi_{+}(n,z)=\Phi_{+}^{(0)}(n)+z^{-1}\Phi_{+}^{(1)}(n)+O(z^{-2}),\quad
z\rightarrow\infty,$

which together with $\Phi_{+}(n,z)=\Phi_{+}^{(0)}(n)\Gamma(n,z)$ imply that
$\Phi_{+}^{(1)}(n)=\Phi_{+}^{(0)}(n)\Gamma^{(1)}(n)$. Note that the solution
$u(n,t)$ can be rewritten as

$u(n,t)=-\frac{\Phi_{+12}^{(1)}(n)}{\Phi_{+22}^{(0)}(n)}=-\frac{\nu_{-}(n)}{\nu_{+}(n)}\Gamma_{12}^{(1)}(n),$
(3.21)

in view of (3.6), (3.8) and the definition of $\Phi_{+}(n,z)$. On the other
hand, since $S^{\dagger}(\bar{z})=S^{-1}(z),\det S(z)=1$, and
$a_{+}(z)=\det\Phi_{+}(n,z)\rightarrow\nu_{+}(n)\nu_{-}(n),~{}z\rightarrow\infty$,
as well as $G(z)=I$ in the RH problem as in (2.27), then

$\nu=\nu_{+}(n)\nu_{-}(n)=1.$ (3.22)

Taking notice of the RH problem reduces to $\Phi_{+}(n,z)=\Phi_{-}(n,z)$,
which allows us to discuss the asymptotic behavior of $\Phi_{+}(n,z)$ near
$z=0$,

$\Phi_{+}(n,z)\rightarrow\tilde{\Phi}_{-}^{(0)}(n),\quad z\rightarrow 0,$
(3.23)

where $\tilde{\Phi}_{-}^{(0)}(n)$ defined in (3.9). Now using again
$\Phi_{+}(n,z)=\Phi_{+}^{(0)}(n)\Gamma(n,z)$, we know that

$\Gamma(n,z=0)=\left(\begin{matrix}\nu_{-}^{-2}(n)&0\\\
0&\nu_{+}^{-2}(n)\end{matrix}\right).$ (3.24)

Hence the solution $u(n,t)$ can be reconstructed from (3.24), (3.22) and
(3.26)

$u(n,t)=-\Gamma_{22}(n,z=0)\Gamma_{12}^{(1)}(n),$ (3.25)

For $N=1$, we have the one soliton solution of dmKdV equation (1.1) from
(3.18) and (3.19) as

$u(n,t)=e^{2\xi}\sinh 2\xi{\rm sech}(2\theta+\xi),$ (3.26)

in terms of the assumption $\eta+2\phi=0$, where $\theta$ and $\phi$ are
defined in (3.15).

##  4 Gauge transformation

In this section, we discuss the Gauge transformation between the plus type
dmKdV equation and the minus one. Here we confine ourselves to the system
(1.4). In [6], the minus type dmKdV equation is the compatibility condition of
the Lax pair

$\chi(n+1)=\tilde{L}_{n}\chi(n),\quad\chi_{t}(n)=\tilde{M}_{n}\chi(n),$ (4.1)

where $\tilde{L}_{n}=(I+U_{n})Z,\tilde{M}_{n}=k\sigma_{3}+\tilde{U}_{n}$ and

$\displaystyle U_{n}=\left(\begin{matrix}0&\tilde{u}(n)\\\
\tilde{u}(n)&0\end{matrix}\right),$ (4.2)
$\displaystyle\tilde{U}_{n}=U_{n}+Z^{-1}U_{n-1}Z-U_{n}U_{n-1}.$

We let $\chi(n)=G_{n}\psi(n)$. If $\psi(n)$ and $Q_{n}$ solve the linear
equations (1.4), then

$\displaystyle\tilde{L}_{n}=G_{n+1}L_{n}G_{n}^{-1},$ (4.3)
$\displaystyle\tilde{M}_{n}=G_{n,t}G_{n}^{-1}+G_{n}M_{n}G_{n}^{-1},$

and

$\tilde{u}(n)=-iu(n-2),$ (4.4)

where $L_{n}=(I+Q_{n})Z,M_{n}=k\sigma_{3}+\tilde{Q}_{n}-Q_{n}Q_{n-1}$ and

$G_{n}=\rho_{n}^{\pm}\left(\begin{matrix}1-iu(n-1)\tilde{u}(n)\lambda&-u(n-1)-i\tilde{u}(n)\lambda\\\
-\tilde{u}(n)\lambda+iu(n-1)\lambda^{2}&u(n-1)\tilde{u}(n)\lambda+i\lambda^{2}\end{matrix}\right),$
(4.5)

with

$\rho_{n}^{+}=\prod\limits_{k=n}^{+\infty}\frac{1+u^{2}(k)}{1-\tilde{u}^{2}(k)},\quad\rho_{n}^{-}=\prod\limits_{k=-\infty}^{n-1}\frac{1-\tilde{u}^{2}(k)}{1+u^{2}(k)}.$

Indeed, $\tilde{L}_{n}=G_{n+1}L_{n}G_{n}^{-1}$ implies that the matrix $G_{n}$
can be represented in the following form

$G_{n}=\left(\begin{matrix}a_{0}(n)+a_{1}(n)\lambda&b_{0}(n)+b_{1}(n)\lambda\\\
c_{1}(n)\lambda+c_{2}(n)\lambda^{2}&d_{1}(n)\lambda+d_{2}(n)\lambda^{2}\end{matrix}\right).$

Then we have a set of equations

$a_{0}(n+1)(1+u^{2}(n))=a_{0}(n)(1-\tilde{u}^{2}(n)),$ (4.6a)
$b_{0}(n+1)=-u(n)a_{0}(n+1),\quad c_{1}(n)=-\tilde{u}(n)a_{0}(n),$ (4.6b)
$\tilde{u}(n)b_{0}(n)=u(n)c_{1}(n+1)+d_{1}(n+1)-d_{1}(n),$ (4.6c)
$d_{2}(n+1)(1+u^{2}(n))=d_{2}(n)(1-\tilde{u}^{2}(n)),$ (4.6d)
$b_{1}(n)=-\tilde{u}(n)d_{2}(n),\quad c_{2}(n+1)=u(n)d_{2}(n+1),$ (4.6e)
$a_{1}(n+1)-a_{1}(n)-u(n)b_{1}(n+1)=\tilde{u}(n)c_{2}(n),$ (4.6f)
$b_{0}(n)=u(n)a_{1}(n+1)+b_{1}(n+1)-\tilde{u}(n)d_{1}(n),$ (4.6g)
$c_{2}(n)=-\tilde{u}(n)a_{1}(n)+c_{1}(n+1)-u(n)d_{1}(n+1).$ (4.6h)

Hence, equation (4.6a) implies $a_{0}(n)=\rho_{n}^{\pm}$, then $b_{0},c_{1}$
and $d_{1}$ can be obtained from (4.6b) and (4.6c). We take
$d_{2}(n)=\alpha\rho_{n}^{\pm}$ by (4.6d), then by (4.6e) and (4.6f),
$b_{1},c_{2}$ and $a_{1}$ is at hand, where $\alpha$ is some constant. Thus
$G_{n}$ in (4.5) is obtained. In addition, the last two equation (4.6g) and
(4.6h) product

$\alpha\tilde{u}(n+1)=u(n-1),\quad\tilde{u}(n+1)=-\alpha u(n-1),$ (4.7)

in view of the identity
$\rho_{n+1}^{\pm}(1+u^{2}(n))=\rho_{n}^{\pm}(1-\tilde{u}^{2}(n))$. Thus (4.4)
is proven. It is remarked that the second equation of (4.3) is valid for
$G_{n}$ given by (4.5), since the gauge transformation between (1.4) and (4.1)
implies

$\tilde{L}_{n,t}-\tilde{M}_{n+1}\tilde{L}_{n}+\tilde{L}_{n}\tilde{M}_{n}=G_{n+1}(L_{n,t}-M_{n+1}L_{n}+L_{n}M_{n})G_{n}^{-1}.$

In this equation, $L_{n,t}-M_{n+1}L_{n}+L_{n}M_{n}=0$ implies the plus type
equation of (1.1), while
$\tilde{L}_{n,t}-\tilde{M}_{n+1}\tilde{L}_{n}+\tilde{L}_{n}\tilde{M}_{n}=0$
gives the minus one.

We note that the gauge transformation about the similar problem of (1.5) gives
rise to

$\tilde{u}(n)=-iu(n-1).$ (4.8)

Hence the solutions of minus type dmKdV equation can be obtained by (4.4) or
(4.8) from the given plus ones. It is interesting to remark that the soliton
solutions to minus type dmKdV equation can be also obtained without needing to
consider the nonvanishing boundary conditions [31], and the solutions are
complex value by equation (4.4) or (4.8) for above two cases.

## Acknowledgments

Projects 11001250 and 11171312 are supported by the National Natural Science
Foundation of China. The work of JY Zhu is partially supported by the
Foundation for Young Teachers in Colleges and Universities of Henan Province.

## References

  * [1] M. J. Ablowitz and J. F. Ladik, Nonlinear differential-difference equations, J. Math. Phys. 16, 598-603 (1975).
  * [2] M. J. Ablowitz and J. F. Ladik, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math. 55, 213-229 (1976).
  * [3] M.J. Ablowitz and J.F. Ladik, Nonlinear differential-difference equations and Fourier analysis, J. Math. Phys. 17, 1011-1018 (1976).
  * [4] M. J. Ablowitz and J. F. Ladik, On the solution of a class of nonlinear partial difference equations, Stud. Appl. Math. 57, 1-12 (1977).
  * [5] E. V. Doktorov, N. P. Matsuka and V. M. Rothos, Perturbation-induced radiation by the Ablowitz-Ladik soliton, Phys. Rev. E 68, 066610 (2003).
  * [6] X.G. Geng and D. Gong, Quasi-periodic solutions of the discrete mKdV hierarchy, Int. J. Geom. Methods M. 10, 1250094 (2013).
  * [7] A. Noguchi, H. Watanabe and K. Sakai, J. Phys. Soc. Jpn. 43, 1441-1446 (1977).
  * [8] K. Narita, Decoupling of the N-soliton solution for a new discrete MKdV equation, Chaos, Solitons Fractals 4, 2237-2244 (1994).
  * [9] S. Palit and A. R. Chowdhury, On the Miura map for discrete integrable systems J. Phys. A: math. Gen. 29, 2861-2867 (1996).
  * [10] C. Chandre, A comparison of two discrete mKdV equations, Physica Scripta 55, 129-130 (1997).
  * [11] K. Narita, Miura transformations between Sokolov-Shabat’s equation and the discrete MKdV equation, J. Phys. Soc. Jpn. 66, 4047-4048 (1997).
  * [12] K. Narita, Multiple shock wave solution for new highly nonlinear difference-differential equation related to the discrete MKdV equation, J. Phys. Soc. Jpn. 68, 1505-1507 (1999).
  * [13] A. Mukaihira, Y. Nakamura, Integrable discretization of the modified KdV equation and applications, Inverse Problems 16, 413-424 (2000).
  * [14] P. Liu, M. Jia and S. Y. Lou, Lax pair and exact solutions of a discrete coupled system related to coupled KdV and coupled mKdV equations, Physica Scripta 76, 674-679 (2007).
  * [15] K. Narita, N-soliton solution of a lattice equation related to the discrete MKdV equation, J. Math. Anal. Appl. 381, 963-965 (2011).
  * [16] M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform (SIAM, Philadelphia, 1981).
  * [17] N. N. Bogolyubov, A. K. Prikarpatskii and V. G. Samoilenko, Discrete periodic problem for the modified nonlinear Korteweg-de Vries equation , Soviet Physics - Doklady 26, 490-492 (1981).
  * [18] Q. Li, Q. Y. Duan and J. B. Zhang, Soliton solutions of the mixed discrete modified Korteweg-de Vries hierarchy via the inverse scattering transform, Physica Scripta 86, 065009 (2012).
  * [19] Q. Li, D. J. Chen, J. B. Zhang and S. T. Chen, Solving the non-isospectral Ablowitz-Ladik hierarchy via the inverse scattering transform and reductions, Chaos Solitons Fractals 45, 1479-1485 (2012).
  * [20] F. W. Nijhoff, G. R. W. Quispel and H. W. Capel, Direct linearization of nonlinear difference-difference equations, Phys. Lett. A 97, 125-128 (1983).
  * [21] J. f. Zhang and C. Q. Dai, Jacobian elliptic function method for nonlinear differential-difference equations, Chaos, Solitons Fractals 27, 1042-1047 (2006).
  * [22] Q. Yang and H. J. Zhang, Exact two-soliton solutions for discrete mKdV equation, Commun. Theor. Phys. 49, 1553-1556 (2008).
  * [23] X. Y. Wen and Y. T. Gao, Darboux transformation and explicit solutions for discretized modified Korteweg-de Vries lattice equation, Commun. Theor. Phys. 53, 825-830 (2010).
  * [24] Y. F. Zhang, J. Q. Mei and Y. C. Hon, Exact soliton solutions of the discrete modified Korteweg-de Vries (mKdV) equation, Phys. Essays 23, 276-284 (2010).
  * [25] Z. Wang and W. X. Ma, Discrete Jacobi sub-equation method for nonlinear differential-difference equations, Math. Method Appl. Sci. 33, 1463-1472 (2010).
  * [26] A. Ankiewicz, N. Akhmediev and J. M. Soto-Crespo, Discrete rogue waves of the Ablowitz-Ladik and Hirota equations, Phys. Review E 82, 026602 (2010).
  * [27] Y. Feng and H. Q. Zhang, Theta function solutions for two discrete equations, Commun. Nonlinear Sci. 15, 2267-2271 (2010).
  * [28] Z. N. Zhu, H. Q. Zhao and X. N. Wu, On the continuous limits and integrability of a new coupled semidiscrete mKdV system, J. Math. Phys. 52, 043508 (2011).
  * [29] K. Narita, Decoupling of the N-soliton solution for a new discrete mKdV equation, Chaos, Solitons Fractals 4, 2237-2244 (1994).
  * [30] S. F. Shen, J. Zhang and Y. X. Wang, New Jocobi-elliptic function solutions of the semi-discrete coupled mKdV system, Phys. Lett. A 343, 148-152 (2005).
  * [31] E. C. Shek and K. W. Chow, The discrete modified Korteweg-de Vries equation with non-vanishing boundary conditions: Interactions of solitons, Chaos Solitons Fractals 36, 296-302 (2008).
  * [32] C. M. Ormerod, Reductions of lattice mKdV to q-P-VI, Phys. Lett. A 376, 2855-2859 (2012).
  * [33] W. X. Ma, Complexiton solutions to the Korteweg-de Vries equation, Phys. Lett. A 301, 35-44 (2002).
  * [34] W. X. Ma and K. Maruno, Complexiton solutions of the Toda lattice equation, Physica A 343, 219-237 (2004).
  * [35] W. X. Ma, Complexiton solutions of the Korteweg-de Vries equation with self-consistent sources, Chaos Solitons Fractals 26, 1453-1458 (2005).
  * [36] W. X. Ma, Complexiton solutions to integrable equations, nonlinear Anal. 63, e2461-e2471 (2005).


# Age-of-Updates Optimization for UAV-assisted Networks

Mouhamed Naby Ndiaye  School of Computer Science, Mohammed VI Polytechnic
University (UM6P), Benguerir, Morocco,
emails<EMAIL_ADDRESS>El Houcine
Bergou  School of Computer Science, Mohammed VI Polytechnic University (UM6P),
Benguerir, Morocco,
emails<EMAIL_ADDRESS>Mounir Ghogho
TICLab, International University of Rabat (UIR), Rabat, Morocco,
email<EMAIL_ADDRESS>Hajar El Hammouti  School of Computer Science,
Mohammed VI Polytechnic University (UM6P), Benguerir, Morocco,
emails<EMAIL_ADDRESS>
###### Abstract

Unmanned aerial vehicles (UAVs) have been proposed as a promising technology
to collect data from IoT devices and relay it to the network. In this work, we
are interested in scenarios where the data is updated periodically, and the
collected updates are time-sensitive. In particular, the data updates may lose
their value if they are not collected and analyzed timely. To maximize the
data freshness, we optimize a new performance metric, namely the Age-of-
Updates (AoU). Our objective is to carefully schedule the UAVs hovering
positions and the users’ association so that the AoU is minimized. Unlike
existing works where the association parameters are considered as binary
variables, we assume that devices send their updates according to a
probability distribution. As a consequence, instead of optimizing a
deterministic objective function, the objective function is replaced by an
expectation over the probability distribution. The expected AoU is therefore
optimized under quality of service and energy constraints. The original
problem being non-convex, we propose an equivalent convex optimization that we
solve using an interior-point method. Our simulation results show the
performance of the proposed approach against a binary association.

###### Index Terms:

3D placements of UAVs, Age-of-Updates, Convex optimization, Users’
association.

##  I Introduction

In the last few years, unmanned aerial vehicles (UAVs) have been considered
for a variety of applications in the telecommunications industry [1]. In
particular, UAVs are used as base stations to extend the network coverage in
dense and out-of-reach areas [2]. They are also deployed as data collectors to
relay information from Internet of Things (IoT) devices to sink nodes [3]. In
this work, we address the problem of UAVs acting as data collectors. We are
interested in the context where the collected data is updated periodically and
the data updates are time-sensitive [4]. In particular, the data updates may
lose their meaning and value if they are not collected and analyzed timely. As
a consequence, it is important to carefully design the UAV’s flight and
stopping points so that the collected updates are kept as fresh as possible.
For this purpose, a new performance metric, namely Age-of-Updates (AoU), has
been introduced. This metric captures the time since the last update was
collected. In this paper, we answer the question: How to schedule the hovering
locations of the drones and determine the selected devices, from which the
data is collected, so that the total AoU is minimized?

###  I-A Related Work

The majority of UAV-enabled networks research works focus on optimizing
metrics such as the sum-rate and energy efficiency [5, 6, 7, 8, 9, 10]. For
example, in [8], the authors design the trajectory of a UAV to maximize the
rate between a drone and a ground base station. The authors consider the
limited energy of the drone, and the dynamic nature of the communication
channel. They decompose the problem into a sequence of control optimizations
that are solved using control theory. Another example can be found in [9]
where the authors propose an energy-efficient drone control policy to ensure
efficient and fair communication coverage for IoT devices using deep
reinforcement learning.

However, only a handful of papers address the problem of data timeliness in
UAV-enabled networks. In [11], the authors optimize a similar concept, called
Age-of-Information (AoI). AoI is defined as the time since the last relevant
information arrived at its destination. The objective of the work in [11] is
to optimize the data collection mode and the UAV trajectory to reduce the
average AoI of IoT devices. In the same context, in [12], the authors propose
a UAV trajectory planning to minimize the maximum AoI of a UAV-enabled
wireless sensor network. Another work in [13] consists of optimizing the
service time allocation and the flight of a UAV used as a mobile relay. The
previously cited works are limited to a single drone application and optimize
only the 2D location of the drone, ignoring the impact of its altitude.

The AoU has been introduced recently in the context of federated learning (FL)
[14]. An age measure of this type has been shown to improve the notion of data
freshness in a variety of applications and has already been used, primarily in
the context of networking (see, for example, the study in [15]). In [14], the
authors propose a device scheduling policy to maximize the parameters’ updates
freshness, and therefore, improve the accuracy of the FL model. The studied
problem is combinatorial. It involves mixed-integer and continuous variables
and it is solved using a decomposition method that addresses the channel
allocation and devices association subproblems separately.

In fact, UAV networks related optimizations are often formulated as mixed-
integer non-linear programming, where the scheduling/association parameters
are considered as binary variables, whereas the positions of UAVs are modelled
as continuous variables. As a result, the optimization problems are non-
convex, most of the time NP-hard, and therefore, difficult to solve. Moreover,
the proposed algorithms are heuristics with no performance guarantee and, in
general, require a long time to converge.

In our paper, we deal with the problem of AoU optimization in UAV-assisted
networks differently. We claim that, the binary association constraint can be
relaxed. The association variables can be seen as probabilities. In
particular, a user decides whether to connect to a UAV or not based on a
probability distribution. As a consequence, instead of optimizing a
deterministic objective function, the objective function is replaced by an
expectation over the probability distribution. This approach has two main
advantages. First, it allows to formulate the problem as a continuous
optimization which can be solved numerically with some guarantees on the final
solution. Such an approach can even result in continuous and convex problems
which can be solved efficiently using polynomial-time algorithms such as
interior-point methods [16]. Second, since the association is probabilistic
(i.e., a device decides to send its updates or not by randomly drawing samples
from a probability distribution), it allows the participation of a large
number of devices, which leads to some fairness among devices.

###  I-B Contribution

In this paper, we consider a UAV-assisted network, where UAVs act as relays
between IoT devices and a macro base station (BS). Due to the limited
communication and energy resources, we assume that only a subset of IoT
devices send their updates to the UAVs. Similarly, only a subset of UAVs
interact with the BS to send the collected updates. The main objective of this
work is to schedule the UAVs’ hovering locations and optimize the association
probabilities (i.e., the probabilities between IoT devices and UAVs, and
between UAVs and the BS) so that the collected updates are kept as fresh as
possible. The contributions of our paper are summarized as follows.

  * •

We formulate the joint association probabilities and 3D locations problem as a
non-linear programming (NLP) optimization where the objective is to minimize
the expected AoU under energy and quality of service constraints.

  * •

To solve the underlying optimization efficiently, we propose a reformulation
of the problem. We show that the new formulation of the problem is convex.
Hence, we propose an interior-point based algorithm which provides a near-
optimal solution to the studied optimization.

  * •

Finally, we validate our proposed approach with simulation experiments. In
particular, we show that the proposed joint probability association and 3D
position scheduling outperforms a benchmark that adopts binary association.

###  I-C Organization

The remainder of the paper is organized as follows. First, the studied system
model is described in Section II. The mathematical formulation of the problem
is given in Section III. In Section IV, we propose a convex reformulation of
the problem. In Section V, we describe the joint association probability and
3D locations to solve the underlying optimization efficiently. Next, in
Section VI, we show the performance of our algorithm for various scenarios.
Finally, concluding remarks are provided in Section VII.

##  II System Model

###  II-A Communication Model

We consider a set $\mathcal{I}$ of $I$ IoT devices within an area
$\mathcal{A}$. The IoT devices collect data from their surrounding environment
and send data updates periodically to a server located at the BS. Since IoT
devices have a limited communication range, a set $\mathcal{U}$ of UAVs are
deployed to collect data on a regular basis from IoT devices and re-transmit
it to the BS. An example of the studied system model is illustrated in Fig. 1.

Figure 1: System Model

We suppose that the IoT devices generate data updates that are transmitted
periodically. At the end of each time interval $k$ of duration $\tau$ (i.e.,
$\left[k\tau,(k+1)\tau\right]$, with $k\in\mathcal{K}=\\{0,1,\dots,K\\}$, and
$KT$ the time horizon), an IoT device $i$ sends its data updates to UAV $u$
with some probability, which in turn transmits the updates to the macro BS in
the following time interval $k+1$ with another probability. Each UAV
communicates with both IoT devices and the BS via the air-to-ground channel.
To model the air-to-ground channel, we assume a Rician fading distribution
$\widehat{\Delta}_{i,u}[k]$

$\widehat{\Delta}_{i,u}[k]=\left(\sqrt{\frac{\Omega}{\Omega+1}}\bar{\Delta}_{i,u}[k]+\sqrt{\frac{1}{\Omega+1}}\widetilde{\Delta}_{i,u}[k]\right),$
(1)

where $\Omega$ is the Rician factor, $\bar{\Delta}_{i,u}[k]$ is the line-of-
sight (LoS) component with $\left|\bar{\Delta}_{i,u}[k]\right|=1$, and
$\widetilde{\Delta}_{i,u}[k]$ the random non-line-of-sight (NLoS) component
where $\widetilde{\Delta}_{i,u}[k]\sim\mathcal{CN}(0,1)$ where
$\mathcal{CN}(0,1)$ is the complex normal distribution.

Each IoT device $i$ transmits its updates with a power $P_{i}^{D}[k]$ during
time interval $k$. Thus, the received power at UAV $u$ is
$\left|\widehat{\Delta}_{i,u}[k]\right|^{2}\beta_{0}\left(d_{i,u}[k]\right)^{-2}P_{i}^{D}[k]$,
where $\beta_{0}$ is the average channel power gain at a reference distance
$d_{0}=1\mathrm{~{}m}$, and $d_{i,u}[k]$ is the distance between device $i$
and UAV $u$ during time slot $k$,

$d_{i,u}[k]=\left\|\textbf{q}_{i}-\textbf{w}_{u}[k]\right\|_{2},\quad
k=1,2,\ldots,K,$ (2)

where $\left\|.\right\|_{2}$ is the two-norm,
$\textbf{q}_{i}=\left[x_{i},y_{i},0\right]$ is the position vector of IoT
device $i$, and $\textbf{w}_{u}[k]=\left[x_{u}[k],y_{u}[k],h_{u}[k]\right]$ is
the 3D location vector of UAV $u$ during time interval $k$.

The signal-to-noise ratio (SNR) of IoT device $i$ with respect to UAV $u$ is
given by

$\Upsilon_{i,u}^{D}[k]=P_{i}^{D}[k]\left|\widehat{\Delta}_{i,u}[k]\right|^{2}\beta_{0}\left(d_{i,u}[k]\right)^{-2}/\sigma^{2},$
(3)

where $\sigma^{2}$ is the thermal noise power. Therefore, the rate of IoT
device $i$ when it is associated with UAV $u$ during time slot $k$ can be
written

$R_{i,u}^{D}[k]=\Lambda_{i,u}^{D}[k]\log_{2}\left(1+\Upsilon_{i,u}^{D}[k]\right),$
(4)

where $\Lambda_{i,u}^{D}[k]$ is the allocated bandwidth between device $i$ and
UAV $u$ during time slot $k$.

Similarly, the rate of UAV $u$ when it transmits the data updates to the BS
during time slot $k$ is given by

$R_{u}^{U}[k]=\Lambda_{u}^{U}[k]\log_{2}\left(1+\Upsilon_{u}^{U}[k]\right),$
(5)

where $\Lambda_{u}^{U}[k]$ is the bandwidth of UAV $u$ during time slot $k$,
and $\Upsilon_{u}^{U}[k]$ is the SNR of UAV $u$ when it transmits to the BS
during time slot $k$ which is given by

$\Upsilon_{u}^{U}[k]=P_{u}^{U}[k]\left|\widehat{\Delta}_{u}[k]\right|^{2}\beta_{0}\left(d_{u}[k]\right)^{-2}/\sigma^{2},$
(6)

with $P_{u}^{U}[k]$ the transmit power of UAV $u$ and $d_{u}[k]$ the distance
between UAV $u$ and the BS during time interval $k$.

###  II-B Energy Consumption Model

Due to their limited battery budget, it is important to account for the energy
that a UAV consumes during its mission. In our work, we account for three
types of energy. First, the flying energy, which is the energy to travel from
one location to another. Second, the hovering energy, which is the energy that
the UAV remains aloft and supports its movement. And finally, the
communication energy which is the transmission energy. Accordingly, the total
consumed energy by a UAV $u$ during time interval $k$ is given by

$\displaystyle E_{u}[k]$
$\displaystyle=E_{u}^{f}[k]+E_{u}^{h}[k]+E_{u}^{c}[k]$
$\displaystyle=P_{u}^{f}\tau_{u}^{f}[k]+(P_{u}^{h}+P_{u}^{U}[k])(\tau-\tau_{u}^{f}[k]),$
(7)

where $P_{u}^{f}$ and $P_{u}^{h}$ are the propulsion power of a rotary-wing
UAV and the hovering power respectively, that we assume constant over time.
$\tau_{u}^{f}[k]$ is the time the UAV $u$ makes to travel from one location to
another during time slot $k$, it is given by
$\tau_{u}^{f}[k]=v^{-1}\left\|\textbf{w}_{u}[k]-\textbf{w}_{u}[k+1]\right\|_{2}$,
where $v$ is the speed of the UAV that we assume constant.

###  II-C Age-of-Updates Metric

The objective of this work is to maximize the total data collected, while
maximizing the freshness of the data and ensure that its UAVs complete their
misson. In particular, our aim is to minimize the AoU.

Let $A_{i,u}$ be the probabilistic event that IoT device $i$ is associated
with UAV $u$. We denote by $a_{i,u}$ the probability that $A_{i,u}=1$. Also,
let $B_{u}$ be the probabilistic event that the UAV $u$ is selected by the BS
and $b_{u}$ is the probability that $B_{u}=1$. More formally, we have:

$A_{i,u}=\left\\{\begin{array}[]{lc}1,&\text{ with probability }a_{i,u}\\\
0,&\text{ with probability }1-a_{i,u}\end{array}\right.$
$B_{u}=\left\\{\begin{array}[]{lc}1,&\text{ with probability }b_{u}\\\
0,&\text{ with probability }1-b_{u}\end{array}\right.$

To model the updates’ freshness, we use the same definition as in [17].
Particularly, if device $i$ is associated with UAV $u$ during time interval
$k$, its AoU evolves as follows,

$T_{i,u}[k+1]=\left(T_{i,u}[k]+1\right)\left(1-A_{i,u}[k]\right),$ (8)

where $T_{i,u}[0]=0$. $T_{i,u}[k]$ is the age of the updates received by drone
$u$ and collected from device $i$ during time interval $k$. Specifically, when
the updates of device $i$ are not collected during time interval $k$ (i.e.,
$A_{i,u}[k]=0$), the AoU is increased by one unit of time. Inversely, when the
updates are transmitted, the AoU is reinitialized to zero.

Therefore, the AoU of IoT devices whose data transit through UAV $u$ and is
received by the BS is given by

$T_{u}[k+1]=\left(\sum_{i\in\mathcal{I}}T_{i,u}[k-1]A_{i,u}[k-1]+1\right)\left(1-B_{u}[k]\right).$
(9)

Finally, the global AoU at the BS, during time interval $k+1$, with $k>1$, can
be written

${AoU}[k+1]=\\!\\!\sum_{u\in\mathcal{U}}\\!\\!\left(\\!\\!\sum_{i\in\mathcal{I}}\left(T_{i,u}[k-1]A_{i,u}[k-1]+1\right)\left(1-B_{u}[k]\right)\\!\\!\right)$
(10)

Due to the random behavior of events $A_{i,u}$ and $B_{u}$, we target to
minimize the expected value of the global AoU as perceived by the BS. In Lemma
12, we give the expression of the studied metric.

$\begin{array}[]{rcl}&\mathbb{E}(AoU[k+1])=\mathbb{E}\Bigg{[}\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\bigg{(}\left(T_{i,u}[k-1]A_{i,u}[k-1]+1\right)\left(1-B_{u}[k]\right)\bigg{)}\Bigg{]}\\\
&=\mathbb{E}\Bigg{[}\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\bigg{(}T_{i,u}[k-1]A_{i,u}[k-1]+(1-B_{u}[k])-T_{i,u}[k-1]A_{i,u}[k-1]B_{u}[k]\bigg{)}\Bigg{]}\\\
&=\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\Bigg{(}\mathbb{E}\bigg{[}T_{i,u}[k-1]A_{i,u}[k-1]\bigg{]}+\mathbb{E}\bigg{[}(1-B_{u}[k])\bigg{]}-\mathbb{E}\bigg{[}T_{i,u}[k-1]A_{i,u}[k-1]B_{u}[k]\bigg{]}\Bigg{)}\\\
&=\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\Bigg{(}T_{i,u}[k-1]\mathbb{E}\Big{[}A_{i,u}[k-1]\Big{]}+\bigg{(}1-\mathbb{E}\Big{[}B_{u}[k]\Big{]}\bigg{)}-T_{i,u}[k-1]\mathbb{E}\Big{[}A_{i,u}[k-1]\Big{]}\mathbb{E}\Big{[}B_{u}[k]\Big{]}\Bigg{)}\\\
&=\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\left(T_{i,u}[k-1]a_{i,u}[k-1]+(1-b_{u}[k])-T_{i,u}[k-1]a_{i,u}[k-1]b_{u}[k]\right)\end{array}\\\
$ (11)

######  Lemma 1.

Under the assumption that $A_{i,u}[k]$ and $B_{u}[k]$ are independent for all
$k$, $u$ and $i$, the expected value of the global AoU, during time interval
$k+1$ is given by

$\displaystyle\mathbb{E}(AoU[k+1])=\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\left(T_{i,u}[k-1]a_{i,u}[k-1]+(1-b_{u}[k])\right.$
(12) $\displaystyle\left.-T_{i,u}[k-1]a_{i,u}[k-1]b_{u}[k]\right).$

###### Proof.

This expression is obtained by developing the AoU and applying the properties
of the expectation, as shown in equation (11).

∎

Next, we formulate the expected AoU optimization as a NLP. Our target is to
find the optimal probability association and 3D positioning of UAVs under
energy and quality of service constraints.

##  III Problem Formulation

In this section, we first describe the constraints of the studied system.
Then, we propose a mathematical formulation of the problem.

First, to ensure that the updates are transmitted timely, it is required that
the expected uplink rate is above a predefined threshold. Let $\bar{R_{I}}$
and $\bar{R_{U}}$ be the rate thresholds for IoT devices and UAVs,
respectively. As a consequence, the quality of service constraints are
expressed as follows

$a_{i,u}[k]R_{i,u}^{D}[k]\geq\bar{R_{I}},\quad\forall i\in\mathcal{I},\forall
k\in\mathcal{K}$ (13)

$b_{u}[k]R_{u}^{U}[k]\geq\bar{R_{U}},\quad\forall u\in\mathcal{U},\forall
k\in\mathcal{K}.$ (14)

Additionally, UAVs must be able to accomplish their mission without exceeding
their energy budget $E_{u}^{\rm max}$. Therefore, the expected consumed energy
should satisfy at each time slot the following constraint

$b_{u}[k]E_{u}[k]\leq E_{u}^{\text{max}},\;\forall u\in\mathcal{U},\forall
k\in\mathcal{K}.$ (15)

Moreover, to avoid any collision between UAVs, we set a minimum distance
$d_{min}$ between any two UAVs, thus

$\left\|\textbf{w}_{u}[k]-\textbf{w}_{v}[k]\right\|_{2}\geq d_{min}\ \forall
u\in\mathcal{U},v\in\mathcal{U},u\neq v,\forall k\in\mathcal{K}.$ (16)

Finally, our target is to optimize the expected AoU over time. Therefore, the
optimization problem is formulated as follows.

$\displaystyle\begin{split}\begin{aligned}
&\mathcal{P}:\min_{\mathbf{a},\mathbf{b},\mathbf{w}}\sum_{k\in\mathcal{K}}\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\left(T_{i,u}[k-1]a_{i,u}[k-1]+(1-b_{u}[k])-\right.\\\
&\left.T_{i,u}[k-1]a_{i,u}[k-1]b_{u}[k]\right)\end{aligned}\end{split}$ (17a)
$\displaystyle\begin{split}\text{s .
t.}\quad(\ref{third:b})-(\ref{dist})\end{split}$ (17b)
$\displaystyle\begin{split}\sum_{u\in\mathcal{U}}a_{i,u}[k]\leq 1,\;\forall
i\in\mathcal{I},\forall k\in\mathcal{K}\end{split}$ (17c)
$\displaystyle\begin{split}0\leq a_{i,u}[k]\leq 1,\forall
i\in\mathcal{I},u\in\mathcal{U},\forall k\in\mathcal{K}\end{split}$ (17d)
$\displaystyle\begin{split}0\leq b_{u}[k]\leq 1,\forall
u\in\mathcal{U},\forall k\in\mathcal{K}\end{split}$ (17e)
$\displaystyle\begin{split}x_{u\text{min}}\leq x_{u}[k]\leq
x_{u\text{max}},\forall u\in\mathcal{U},\forall k\in\mathcal{K}\end{split}$
(17f) $\displaystyle\begin{split}y_{u\text{min}}\leq y_{u}[k]\leq
y_{u\text{max}},\forall u\in\mathcal{U},\forall k\in\mathcal{K}\end{split}$
(17g) $\displaystyle\begin{split}z_{u\text{min}}\leq z_{u}[k]\leq
z_{u\text{max}},\forall u\in\mathcal{U},\forall k\in\mathcal{K}\end{split}$
(17h)

where $\mathbf{a}$ and $\mathbf{b}$ are the association probability vectors
and $\mathbf{w}$ is the 3D position vector over time. In this formulation, the
constraint (17c) ensures that the association probabilities over UAVs do not
exceed $1$. The constraints (17d) and (17e) ensure that the probability
vectors $\mathbf{a}$ and $\mathbf{b}$ are well-defined. Finally, the
constraints (17f), (17g) and (17h) restrict the movement of the UAVs to a
limited $3D$ space.

The Problem $\mathcal{P}$ includes only continuous variables, which makes it a
continuous non-linear and non-convex programming problem. It is not convex,
because constraints (13), (14) and the objective function are not convex. In
the next section, we propose a convex reformulation of Problem $\mathcal{P}$
that can be solved efficiently using methods of convex programming [16].

##  IV A Convex Reformulation of AoU Optimization

In this section, we propose a new convex reformulation of Problem
$\mathcal{P}$. For all $k$, to deal with the non-convexity factor
$a_{i,u}[k-1]b_{u}[k]$ in the objective function, we introduce a slack
variable $t_{k}=a_{i,u}[k-1]b_{u}[k]$. Since $t_{k}$ is the multiplication of
two probabilities, it should be between $0$ and $1$. The objective function
using the new variables then becomes

$\displaystyle\sum_{k\in\mathcal{K}}\sum_{u\in\mathcal{U}}\sum_{i\in\mathcal{I}}\left(T_{i,u}[k-1]a_{i,u}[k-1]+(1-b_{u}[k])-\right.$
$\displaystyle\left.T_{i,u}[k-1]t_{k}\right),$

which is a linear function of the decision variables. The following lemma
deals with the non convexity of the constraint (13). The constraint (14) is
similar to (13), thus the same reformulation applies.

######  Lemma 2.

For each $k\in\mathcal{K}$, the IoT device’s rate constraint (13), can be
rewritten as

$\frac{e^{\left(\frac{\bar{R_{I}}}{a_{i,u}[k]\Lambda_{i,u}^{D}}\right)}-1}{P_{i}^{D}/\sigma^{-2}\left|\widehat{\Delta}_{i,u}[k]\right|^{2}\beta_{0}}-t_{1}\leq
0,\text{ where }$ (18)
$\frac{1}{(x_{u}[k]-x_{i})^{2}+(y_{u}[k]-y_{i})^{2}+(z_{u}[k])^{2}}-t_{1}=0.$
(19)

The functions defining (18) and (19) are both convex.

###### Proof.

After applying the exponential function, the constraint (13) becomes

$\Upsilon_{i,u}^{D}[k]\geq
e^{\left(\frac{\bar{R_{I}}}{b_{u}[k]\Lambda_{i,u}^{D}}\right)}-1,$

where
$\Upsilon_{i,u}^{D}[k]=P_{i}^{D}[k]\left|\widehat{\Delta}_{i,u}[k]\right|^{2}\beta_{0}\sigma^{-2}\left(d_{i,u}[k]\right)^{-2}$.
Now, by letting
$t_{1}=\frac{1}{(x_{u}[k]-x_{i})^{2}+(y_{u}[k]-y_{i})^{2}+(z_{u}[k])^{2}},$ we
get

$\frac{e^{\left(\frac{\bar{R_{I}}}{a_{i,u}[k]\Lambda_{i,u}^{D}}\right)}-1}{P_{i}^{D}[k]/\sigma^{-2}\left|\widehat{\Delta}_{i,u}[k]\right|^{2}\beta_{0}}-t_{1}\leq
0.$

Finally, note that the functions $h_{1}(x)=e^{\frac{1}{x}}$, and
$h_{2}(x)=\frac{1}{x^{2}}$ are both convex.

∎

Note all the remaining constraints in the definition of Problem $\mathcal{P}$
are convex.

Now that we showed the convex reformulation of Problem $\mathcal{P}$. In the
next section, we outline the main steps of our proposed approach to perform
the Joint probability Association and 3D position Scheduling (JAS) of IoT
devices, UAVs and the BS.

##  V Joint probability Association and position Scheduling for AoU
Minimization

Figure 2: Comparison between
JAS and Deterministic approach

Figure 3: Expected AoU over iterations

Figure 4: Energy consumption by each UAV

Figure 5: UAV trajectory and association for $k=5$

Now that we have succeeded in getting a convex reformulation of $\mathcal{P}$,
we use convex optimization methods to solve it efficiently, namely the
interior-point method.

After obtaining the optimal probabilities, in practice at each time interval
$k$, we use the procedure U2B to decide the association between the UAVs and
the BS, and I2U to decide the association of IoT devices with the UAVs.

U2B procedure: Let $\textbf{b}^{*}$ be the optimal UAV-BS association vector
obtained by solving the convex optimization. To link a given UAV $u$ to the BS
during time interval $k$, we draw a random number between $0$ and $1$. If this
number is lower than the obtained optimal association probability
$b_{u}[k]^{*}$, UAV $u$ is associated with the BS during time interval $k$ and
therefore is able to transmit the collected updates to the BS.

I2U procedure: Let $\textbf{a}^{*}$ be the optimal IoT device-UAV association
matrix obtained by solving the convex optimization problem. We assume that
during a time interval $k$, a device $i$ is allowed to send its updates to at
most one UAV. The question that arises here is how to select at most one UAV
according to the association probability vector $\textbf{a}^{*}$? For ease of
notations, we remove the dependency on $k$. The probability selection vector
of device $i$ with respect to the UAVs is given by
$\textbf{a}_{i}^{*}=(a_{i,1}^{*},\dots,a_{i,U}^{*})$. To pick at most one UAV
according to the probability vector $\textbf{a}_{i}^{*}$, we construct a
virtual set of UAVs, where each UAV is represented according to its
probability of selection. For example, let $U=3$ and
$\textbf{a}_{i}^{*}=(0.3,0.2,0.1)$. To select its associated UAV, device $i$
has to pick uniformly a random one UAV from the set
$\\{\text{UAV}_{1},\text{UAV}_{1},\text{UAV}_{1},\text{UAV}_{2},\text{UAV}_{2},\text{UAV}_{3},\text{UAV}_{0},\text{UAV}_{0},\\\
\text{UAV}_{0},\text{UAV}_{0}\\}$ with $\text{UAV}_{0}$ is a hypothetical UAV
which simply indicates that device $i$ will not be associated during that time
interval.

Algorithm 1 summarizes the different steps of our approach.

Algorithm 1 Joint probability Association and position Scheduling (JAS)

1:Input:  The maximum energy budget: $E_{u}^{max}\,~{}\forall u\in\mathcal{U}$
, the minimum IoT rate: $\bar{R_{I}}$, the minimum UAV rate: $\bar{R_{U}}$,
the minimum distance between UAV: $d_{min}$. $K$ number of communication
rounds.

2:Optimization: Solve the convex reformulation of Problem $\mathcal{P}$ using
a convex programming method, such as the interior-point method, to obtain the
optimal association probabilities $\mathbf{a}^{*}$ and $\mathbf{b}^{*}$; and
the optimal UAVs’ positions over time.

3:Association:

4:while $k\leq K$ do

5: Associate IoT devices to UAVs using I2U procedure

6: Connect UAVs to BS using U2B procedure

7: $k=k+1$

##  VI Simulation Results

To evaluate the performance of our approach, we consider the following
scenario. We assume $K=10$, the number of UAVs $U=4$, the number of IoT
devices $I=25$, the speed to $v=15m/s$, the minimum distance between two UAVs
$d_{min}=15m$, and a 3D area of $100\times 100\times 100m^{3}$. Initially, the
UAVs are randomly scattered in the 3D space. We compare the proposed approach
with a deterministic approach similar to the one proposed in [17].The results
of this comparison are shown in Fig. 5. As one may expect, our method gives a
smaller value of AoU when compared to the deterministic approach, which in
turn is better than a random feasible solution of the optimization problem.

In Fig. 5, we plot the objective value over iterations. As the plot shows, the
value of the function decreases over the iterations until it reaches the
optimum.

Fig. 5 shows the energy consumed by each UAV against their maximum energy. We
can see that none of the UAVs exceeded their energy limits.

In Fig. 5, we plot the UAVs trajectory and the associations for the last time
interval $k=5$. In this figure, we show the 3D positioning of the UAVs as well
as the trajectory followed by each UAV from $k=1$ to $k=5$. We also represent
the associations between each UAV and the IoT devices as well as the UAVs
selected by the BS for transmission during the last time interval. We can
notice that, the UAVs try to position themselves in order to communicate well
with the IoT devices and the BS. We can also see that most of the IoT devices
are served by the UAVs.

##  VII Conclusion

In this paper, we consider the optimization of the expected AoU in a UAV-
assisted network. Our objective is to find the optimal 3D location scheduling
and the probabilities of association between the IoT devices, the UAVs, and
the BS. We first formulate the problem as a non-linear programming. Then, we
propose a convex reformulation of the studied optimization. We solve the
convex optimization using an interior-point based algorithm. Finally, we show
that the proposed approach outperforms the benchmark algorithm.

## Acknowledgment

Research reported in this work was supported by the European Commission in
partnership with the African Academy of Sciences (AAS) through ARISE-PP Grant
(Grant N°:ARISE-PP-FA-109).

## References

  * [1] H. Price. (2019) Federal Aviation Administration (FAA) Forecast Fiscal Years 2017-2038. [Online]. Available: https://www.faa.gov/news/fact_sheets/news_story.cfm?newsId=22594/
  * [2] H. E. Hammouti and M. Ghogho, “Air-to-ground channel modeling for uav communications using 3d building footprints,” in _Ubiquitous Networking_ , N. Boudriga, M.-S. Alouini, S. Rekhis, E. Sabir, and S. Pollin, Eds. Cham: Springer International Publishing, 2018, pp. 372–383.
  * [3] M. Alzenad, A. El-Keyi, F. Lagum, and H. Yanikomeroglu, “3-D placement of an unmanned aerial vehicle base station (UAV-BS) for energy-efficient maximal coverage,” _IEEE Wireless Communications Letters_ , vol. 6, no. 4, pp. 434–437, 2017.
  * [4] N. Wang, Y. Xin, J. Zheng, J. Wang, X. Liu, X. Hou, and Y. Liu, “Priority-oriented trajectory planning for UAV-aided time-sensitive IoT networks,” in _2020 IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, 2020, pp. 1–7.
  * [5] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 6, pp. 3747–3760, 2017.
  * [6] Y. Zeng, J. Xu, and R. Zhang, “Energy minimization for wireless communication with rotary-wing UAV,” _IEEE Transactions on Wireless Communications_ , vol. 18, pp. 2329–2345, 4 2019.
  * [7] C. Zhan and Y. Zeng, “Completion time minimization for multi-UAV-enabled data collection,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 10, pp. 4859–4872, 2019.
  * [8] D. B. Licea, M. Bonilla, M. Ghogho, S. Lasaulce, and V. S. Varma, “Communication-aware energy efficient trajectory planning with limited channel knowledge,” _IEEE Transactions on Robotics_ , vol. 36, no. 2, pp. 431–442, 2019.
  * [9] C. H. Liu, Z. Chen, J. Tang, J. Xu, and C. Piao, “Energy-efficient UAV control for effective and fair communication coverage: A deep reinforcement learning approach,” _IEEE Journal on Selected Areas in Communications_ , vol. 36, no. 9, pp. 2059–2070, 2018.
  * [10] A. Ranjha and G. Kaddoum, “Urllc-enabled by laser powered uav relay: A quasi-optimal design of resource allocation, trajectory planning and energy harvesting,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 1, pp. 753–765, 2021.
  * [11] Z. Jia, X. Qin, Z. Wang, and B. Liu, “Age-based path planning and data acquisition in UAV-assisted IoT networks,” in _2019 IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, 2019, pp. 1–6.
  * [12] P. Tong, J. Liu, X. Wang, B. Bai, and H. Dai, “UAV-enabled age-optimal data collection in wireless sensor networks,” in _2019 IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, 2019, pp. 1–6.
  * [13] M. A. Abd-Elmagid and H. S. Dhillon, “Average peak age-of-information minimization in UAV-assisted IoT networks,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 2, pp. 2003–2008, 2018.
  * [14] H. H. Yang, A. Arafa, T. Q. Quek, and H. V. Poor, “Age-based scheduling policy for federated learning in mobile edge networks,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 8743–8747.
  * [15] R. D. Yates, Y. Sun, D. R. Brown, S. K. Kaul, E. Modiano, and S. Ulukus, “Age of information: An introduction and survey,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 5, pp. 1183–1210, 2021.
  * [16] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [17] H. H. Yang, A. Arafa, T. Q. S. Quek, and H. V. Poor, “Age-based scheduling policy for federated learning in mobile edge networks,” 10 2019. [Online]. Available: http://arxiv.org/abs/1910.14648