text:
\begin{document}

\vspace{-0.5cm}

\title{Study of BEM-Type Channel Estimation Techniques for 5G Multicarrier Systems}

\author{\IEEEauthorblockN{Joao T. Dias$^{1,2}$,
Rodrigo C. de Lamare$^{2,3}$ and
 Yuriy V. Zakharov$^{3}$} \\%\IEEEauthorrefmark{3}}
\vspace{5mm}
\IEEEauthorblockA{$^{1}$Federal Center for Technological Education of Rio de Janeiro (CEFET-RJ)},\\
$^{2}$Centre for Telecommunications Studies (CETUC) of Pontifical Catholic University of Rio de Janeiro, Brazil and\\
$^{3}$Department of Electronic Engineering, University of York, UK\\
e-mail: joao.dias@cefet-rj.br, $\left\{ \text{joao.dias, delamare}\right\}$@cetuc.puc-rio.br and $\left\{ \text{rodrigo.delamare, yury.zakharov}\right\}$@york.ac.uk}
\maketitle

\thispagestyle{empty}

\begin{abstract} ~ In this paper, we investigate channel estimation techniques for 5G multicarrier systems. Due to the characteristics of the 5G application scenarios, channel estimation techniques have been tested in Orthogonal Frequency Division Multiplexing (OFDM) and Generalized Frequency Division Multiplexing (GFDM) systems. The orthogonality between subcarriers in OFDM systems permits inserting and extracting pilots without interference. However, due to pulse shaping, subcarriers in GFDM are no longer orthogonal and interfere with each other. Due to such interference, the channel estimation for GFDM is not trivial.  A robust and low-complexity channel estimator can be obtained by combining a minimum mean-square error (MMSE) regularization and the basis expansion model (BEM) approach. In this work, we develop a BEM-type channel estimator along with a strategy to obtain the covariance matrix of the BEM coefficients. Simulations show that the BEM-type channel estimation shows performance close to that of the linear MMSE (LMMSE), even though there is no need to know the channel power delay profile, and its complexity is low.
\end{abstract}

\begin{IEEEkeywords}~$5$G waveforms, BEM, Channel estimation, Multicarrier Systems.
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction}

The fifth generation (5G) of mobile systems must handle diverse
scenarios and meet the increasing  demand for high data rates. To
achieve an efficient use of limited spectrum, Cognitive Radio (CR)
has gained a great deal of interest.  The idea is to allow the
exploitation of locally available frequency bands, on a temporary
basis but under a non-interfering constraint. In order to avoid the
interference from an opportunistic user to the primary user, a
waveform used by the secondary user must have a low out of band
radiation (OOB)~\cite{Nicola}. Most modern digital communication
standards use Orthogonal Frequency Division Multiplexing (OFDM)
\cite{Ahmad} as the air interface, because of its flexibility and
robustness in frequency-selective channels. Nevertheless, OFDM
presents some drawbacks such as the high OOB \cite{Jaap} that
affects its application especially in CR systems.

In contrast, GFDM can achieve a low OOB using a spectrally-contained pulse on each subcarrier, and it has higher spectrum efficiency, because it does not need to use virtual subcarriers to avoid adjacent channel interference, and because it reduces the ratio between the guard time interval and the total frame duration~\cite{Nicola, Dias}.

However, the subcarrier filtering results in non-orthogonal subcarriers, leading to inter-symbol interference (ISI) and inter-carrier interference (ICI).

In this work, we develop a basis expansion model (BEM) approach for
estimation of the channel parameters~\cite{Edfor, Zakharov, Yuri} in
OFDM and GFDM systems. We first consider classical least-squares
(LS) and linear minimum mean-square error (LMMSE), and then examine
BEM-type channel estimation algorithms. We then use a strategy to
estimate the covariance matrix of the channel vector containing the
coefficients of the channel and incorporate it into a BEM-type
channel estimation technique that approximates the LMMSE channel
estimator. Simulations show that the BEM estimator presents results
close to those of the LMMSE estimator, but there is no need to know
the power delay profile (PDP) of the channel and the estimator
complexity is feasible.

The paper is organized as follows: in Section II, a description of
the signal model of the OFDM and GFDM systems is given; Section III
describes estimation algorithms known in the literature; the BEM
algorithms are described in Section IV; Section V describes the
covariance matrix approach and aproximated-LMMSE-BEM algorithm; the
results of the simulations are presented in Section VI and in
Section VII, some conclusions are drawn.

\section{Signal Model}
\label{ProbState}

The OFDM transmit signal is given by \cite{Ahmad}
\begin{equation} \label{eq:ofdmtransmit}

x_{\text{OFDM}} [n]=\displaystyle \sum_{k=0}^{K-1} d_{k}e^{j2\pi \frac{k}{K}n},
\end{equation}
where $d_{k}$ is the data symbol at the $k^{th}$ subcarrier and $K$
is the number of subcarriers in the OFDM symbol.

In GFDM, each data symbol is pulse-shaped using a circularly time
and frequency shifted version of a prototype filter $g[n]$ whose
energy over one period $N=MK$ is normalized to one ($M$ is the
number of subsymbols per GFDM block and $K$ is the number of
subcarriers per subsymbol). The GFDM transmit signal is given by
\cite{Nicola}
\begin{equation} \label{eq:gfdmtransmit}

x_{\text{GFDM}} [n]=\displaystyle \sum_{k=0}^{K-1}\displaystyle \sum_{m=0}^{M-1}d_{k}[m] g[(n-mK)~ \text{mod}~ N]e^{j2\pi \frac{k}{K}n},
\end{equation}
where $d_{k}[m]$ is the data symbol at time index $m$ and subcarrier index $k$.

\begin{figure*}[!th]
\centering
\includegraphics[width=.8\textwidth]{fig1.eps}
\caption{Block diagram of the proposed solution for the OFDM and the GFDM systems}
\label{block_diagram_ofdm_gfdm}
\end{figure*}

In order to insert and extract interference-free pilots, we rewrite the GFDM transmit signal by \cite{Shahab}
\begin{equation} \label{eq:gfdmtransmit}
\bold{x}_{\text{GFDM}}=\bold{W}_N^H\displaystyle \sum_{k=0}^{K-1}\bold{P}^{(k)}\bold{G}^{(\delta)}\bold{S}^{(\delta)} \left(\bold{\Gamma} \tilde{\bold{d}}_k+\bold{W}_M \check{\bold{d}}_k \right)  \in \mathbb{C}^{N  \times 1},
\end{equation}

where $\check{\bold{d}}_k$ contains the symbols at data subcarriers with zeros at pilot subcarriers and $\tilde{\bold{d}}_k$ carries information at only the pilot subcarriers with zeros at data subcarriers.
$\bold{W}_M$ is the unitary DFT matrix of size $M\times M$ and $\bold{\Gamma}=\bold{P'}\text{blkdiag}(\lambda\bold{I}_{n_p},\bold{W}_{M-n_p})$. Here, $\bold{P'}$ can be any permutation matrix of compatible size which allocates the pilots to any frequency bin within the pilot subcarriers, and $\text{blkdiag}(\bold{X},...,\bold{Y})$ is a block diagonal matrix according to its matrix entries with $\bold{X}$ being the top-left and $\bold{Y}$ being the bottom-right blocks. The parameter $\lambda$ is a scaling factor that normalizes the pilot energy to one, and $\bold{I}_{n_p}$ is an identity matrix of size $n_p$, where $n_p$ is the number of pilots in each $\tilde{\bold{d}}_k$. $\bold{S}^{(\delta)}=\bold{1}_{\delta,1}\otimes\bold{I}_M$ is $\delta$-fold repetition matrix which concatenates $\delta$ identity matrices of size $M$; $\bold{1}_{i,j}$ is a $i\times j$ matrix of ones, and $\otimes$ is the Kronecker product. The value of $\delta$ is based on the number of non-zero values in the filter frequency response. The subcarrier filter $\bold{G}^{(\delta)}=\text{diag}(\bold{W}_{M\delta}\bold{g}^{\delta})$ is diagonal in the frequency domain. The circulant filter $\bold{g}^{\delta}$ is the version of $\bold{g}=(g[n])_{n=0,...N-1}$ down-sampled by factor $K/\delta$. An up-conversion of the $k^{th}$ subcarrier to its respective subcarrier frequency is performed by the shift matrix $\bold{P}^{(k)}=\bold{\Psi}(\bold{p}^{(k)})\otimes\bold{I}_M$, where $\bold{\Psi}(.)$ returns the circulant matrix based on the input vector, and $\bold{p}^{(k)}$ is a column vector where the $k^{th}$ element is $1$ and all others are zeros. The K subcarriers are summed and transformed back to the time-domain with $\bold{W}_N^H$, where $(.)^H$ is the Hermitian transpose of (.).

To ensure the interference-free pilot insertion, in expression
(\ref{eq:gfdmtransmit}) the plus sign cannot superimpose the
information, i.e., if $k$ is a pilot subcarrier $\bold{W}_M
\check{\bold{d}}_k$ becomes $\bold{0}_M$ and if it belongs to data
subcarriers $\bold{\Gamma} \tilde{\bold{d}}_k$, it would be
$\bold{0}_M$. The pilots need to be located at the frequency bins
where no inter-carrier interference is present.

At the receiver, the received signal is described by
\begin{equation} \label{eq:signalreceived}
\bold{y}'=\bold{H}\bold{x}' +\bold{n}  \in \mathbb{C}^{(KM+L)\times 1}
\end{equation}
where $\bold{H} \in \mathbb{C}^{(KM+L)\times (KM+L)}$ describes the channel, $\bold{x}'$ is the vector $\bold{x}$ with addition of the cyclic prefix (CP), and $\bold{n}$ is a noise vector containing additive white Gaussian noise (AWGN) samples.

The signal is converted from the analog to the digital domain, followed by removal of the CP. Then, we apply
the Fourier transform to extract the pilots in the frequency domain,
to estimate the channel $\bold{\hat{H}}$ and to obtain
$\bold{y}=\bold{\hat{H}}^{-1}\bold{y}'$.

On the demodulator side, the recovered data symbols for the $k^{th}$ subcarrier in GFDM system are given by
\begin{equation} \label{eq:gfdmreceive}
\bold{\hat{d}}_{k}=\bold{W}_M^H(\bold{S}^{(\delta)})^T(\bold{G}^{(\delta)})^{\ast}(\bold{P}^{(k)})^T \bold{W}_N \bold{y} ~~ \in \mathbb{C}^{M  \times 1}.
\end{equation}

In order to mitigate the interference between adjacent data
subcarriers, we use the interference cancellation (IC) algorithm
proposed  in \cite{Datta}. Note that other interference cancellation
techniques
\cite{delamare_mber,rontogiannis,delamare_itic,stspadf,choi,stbcccm,FL11,delamarespl07,jidf,jio_mimo,tds,peng_twc,spa,spa2,jio_mimo,P.Li,jingjing,memd,did,bfidd,mbdf,bfidd,mserrr,mmimo,wence,shaowcl08}.
can also be examined. In this case, the notation needs to be
extended by an iteration index $j=1,...,J$. The vector of received
frequency domain samples $\bold{y}_k$ then corresponds to
$\bold{y}_k^{(j)}$, and the received data symbols
$\bold{\hat{d}}_{k}$ to $\bold{\hat{d}}_{k}^{(j)}$, respectively.
With this notation, $\bold{y}_k^{(0)}$ and
$\bold{\hat{d}}_{k}^{(0)}$ denote vectors on which no interference
cancellation has been performed. The IC algorithm is detailed in
Algorithm 1.

\begin{algorithm}[]
\caption{IC}
\begin{algorithmic}[1]
\STATE{Consider $\bold{y}_k=(\bold{S}^{(\delta)})^T(\bold{G}^{(\delta)})^{\ast}(\bold{P}^{(k)})^T \bold{W}_N \bold{y}$}
\STATE{receive all subcarriers as $\bold{W}_{M}^H\bold{y}_k^{(0)}$}
\STATE{map each symbol to the closest Quadrature Phase Shift Keying (QPSK) point to obtain $\bold{\hat{d}}_{k}^{(0)}$}

    \FOR{$j=1:J$}
      \FOR{$k=0:K-1$}

           \STATE{Remove interference by computing: $\bold{y}_k^{(j)}=\bold{y}_k^{(0)}-\bold{P}^{(k)}\bold{G}^{(\delta)}\bold{S}^{(\delta)}\bold{W}_M\left(\bold{\hat{d}}_{(k-1)~ \text{mod}~ K}^{(j-1)}+\bold{\hat{d}}_{(k+1)~ \text{mod}~ K}^{(j-1)}\right)$}
           \STATE{update the received symbol to the closest QPSK point to obtain $\bold{\hat{d}}_{k}^{(j)}$}
     \ENDFOR
    \ENDFOR
\end{algorithmic}
\end{algorithm}

Essentially, the process consists of detecting the data symbols $\bold{\hat{d}}_{k}^{(j)}$ and using them in the $(j+1)$th iteration to compute the interfering signal $\bold{P}^{(k)}\bold{G}^{(\delta)}\bold{S}^{(\delta)}\bold{W}_M\left(\bold{\hat{d}}_{(k-1)~ \text{mod}~ K}^{(j-1)}+\bold{\hat{d}}_{(k+1)~ \text{mod}~ K}^{(j-1)}\right)$, which is then subtracted from the original signal $\bold{y}_k^{(0)}$ in order to obtain an interference-reduced version of the data estimate from $\bold{y}_k^{(j)}$.

The analysis of the possible error propagation issue in the interference cancellation performed by algorithm 1 is beyond the scope of this work. However, it was possible to observe in the simulations carried out that such error propagation can be considered negligible in the scenario evaluated.

The computational complexity of the systems (OFDM and GFDM) in terms of complex valued multiplications can be expressed as: $\mathcal{O}_{OFDM} (2MK\text{log}_2K)$, where $2$ is due to the fact that this operation occurs in the transmitter  and receiver, $M$ is the number the symbol per block and $K\text{log}_2K$ originates from the $K\times K$ points DFT (see Eq. (\ref{eq:ofdmtransmit})); and $\mathcal{O}_{GFDM} (2(N\text{log}_2N+KLM+KM\text{log}_2M))$, where $KLM$ denotes the matched filtering of the sub-carriers, $N\text{log}_2N$ and $M\text{log}_2M$ originates from the $N\times N$ and  $M\times M$ points DFT, respectively (see Eqs. (\ref{eq:gfdmtransmit}) and (\ref{eq:gfdmreceive})).

Applying the IC algorithm $J$ times to all sub-carriers introduces additional $JKM\text{log}_2M+JKM+JKM\text{log}_2M$ operations for GFDM system, where $JKM\text{log}_2M$ for transforming the estimated data symbols to the frequency domain, $JKM$ for applying the interference filter and another $JKM\text{log}_2M$ for transforming back to the time domain.

The block diagram of the proposed solution for the OFDM and the GFDM systems is shown in Fig.~\ref{block_diagram_ofdm_gfdm}.

 \section{Classical Channel Estimation Algorithms}
 \label{UwProject}
In this section, we describe the channel estimation methods that are utilized for comparison with the BEM based technique.

 At the receiver side, after performing the time-frequency conversion, we can rewrite~(\ref{eq:signalreceived}) in the pilot positions
 , as:
\begin{equation} \label{eq:signaluwreceived}
\bold{\tilde{y}}_{p}= \bold{\tilde{H}}_{p}\bold{\tilde{x}}_{p} + \bold{\tilde{n}}_{p},
\end{equation}
where $\bold{\tilde{y}}_{p} \in \mathbb{C}^{(N_p)\times 1}$, $\bold{\tilde{H}}_{p} \in \mathbb{C}^{(N_p)\times (N_p)}$, $\bold{\tilde{x}}_p \in \mathbb{C}^{(N_p )\times 1}$ and $N_p$ is the number of pilots per symbol.

\subsection{Least Squares}
Least Squares (LS) \cite{Ozdemir}
provides the estimate of the channel frequency response at the pilot positions by
\begin{equation} \label{eq:ls}
\bold{\hat{\tilde{h}}_p}_{\text{LS}}= \text{diag}\left\{\bold{\tilde{x}}_{p}\right\}^{-1}\bold{\tilde{y}}_{p},
\end{equation}
where the matrix $(\text{diag}\left\{\bold{\tilde{x}}_{p}\right\}^{-1}) \in \mathbb{C}^{(N_p)\times (N_p)}$ and the vector $\bold{\tilde{x}}_{p}$ occupies the main diagonal of the matrix.

The LS algorithm has low computational complexity $\mathcal{O} (N_p)$,  but the estimates have a high variance due to the noise presence.

\subsection{Linear Minimum Mean Square Error}
 The Linear Minimum Mean Square Error (LMMSE) \cite{Edfor}

can be written as

 \begin{equation} \label{eq:lmmse2}
 \bold{\hat{\tilde{h}}_p}_{\text{LMMSE}}= \bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}} \left(\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}} + \frac{\beta}{SNR}\bold{I}_{N_p}\right)^{-1}\bold{\hat{\tilde{h}}_p}_{\text{LS}},
\end{equation}
where $\beta=E\{ |d_k|^2\}E\{\frac{1}{|d_k|^2}\}$, with $d_k$ being the constellation points of the mapping, $SNR$ is the signal to noise ratio, and $\bold{I}_{N_p}$ is an identity matrix.
To compute an element of the channel covariance matrix $\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}}=E\{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}^H\}$, we need to know the PDP of the channel.
The need to estimate the channel covariance together with the need to estimate the SNR periodically, makes the use of the LMMSE estimator computationally prohibitive; its computational complexity is $\mathcal{O} (N_p^3)$.

\section{BEM Channel Estimation Algorithms}

The BEM approach represents the channel frequency response using basis functions
\cite{Hrycak}. It should be noted that other subspace techniques
\cite{scharf,bar-ness,pados99,reed98,hua,goldstein,santos,qian,delamarespl07,xutsa,delamaretsp,kwak,xu&liu,delamareccm,wcccm,delamareelb,jidf,delamarecl,delamaresp,delamaretvt,jioel,delamarespl07,delamare_ccmmswf,jidf_echo,delamaretvt10,delamaretvt2011ST,delamare10,fa10,lei09,ccmavf,lei10,jio_ccm,ccmavf,stap_jio,zhaocheng,zhaocheng2,arh_eusipco,arh_taes,dfjio,rdrab,dcg_conf,dcg,dce,drr_conf,dta_conf1,dta_conf2,dta_ls,song,wljio,barc,jiomber,saalt}
exploit the eigenstructure of the input data matrix can be also
considere. Using BEM, the vector $\bold{\tilde{h}}$ of the channel
frequency response can be represented as
\begin{equation}
\bold{\tilde{h}} = \bold{B}\bold{a},
\end{equation}
where $\bold{B}$ is the matrix of basis functions, every column of $\bold{B}$ is a basis function, and $\bold{a}$ is the vector of BEM expansion coefficients.

When the channel estimation is transformed into estimation of the vector $\bold{a}$, its complexity is reduced because the size of $\bold{a}$ is smaller  than $\bold{\tilde{h}}$.

\subsection{LS-BEM}

To get an estimate $\bold{\hat{a}}$, we can rewrite (\ref{eq:signaluwreceived}) by replacing $\bold{\tilde{h}}$ by  $\bold{B}\bold{a}$ and  we compute $\bold{a}$ that minimizes the cost function $J(\bold{a})= (\bold{\tilde{y}}_{p}- \text{diag}(\bold{\tilde{x}}_{p})\bold{B}\bold{a})^H(\bold{\tilde{y}}_{p}- \text{diag}(\bold{\tilde{x}}_{p})\bold{B}\bold{a})$. This procedure is known as Least-Squares BEM (LS-BEM) \cite{Barhumi}. The minimum of $J(\bold{a})$ is found by differentiating it with respect to $\bold{a}$ and setting the derivative equal to a null vector. The resultant LS estimate of $\bold{a}$ is
\begin{equation}

\bold{\hat{a}} = (\bold{B}^H\bold{\tilde{X}}^H\bold{\tilde{X}}\bold{B})^{-1}\bold{B}^H\bold{\tilde{X}}^H\bold{\tilde{y}}_{p},

\end{equation}
where $\bold{\tilde{X}}=\text{diag}(\bold{\tilde{x}}_{p})$.

Since each pilot subcarrier has modulus equal to $1$ ($|\tilde{\text{x}}_{p}(k)| = 1$), then $\bold{\tilde{X}}^H\bold{\tilde{X}}= \bold{I}$ is an identity matrix. As the matrix $\bold{B}$ with the bases is known a priori, the matrix $\bold{B}^H\bold{B}$ can be precomputed. Its inverse $(\bold{B}^H\bold{B})^{-1}$ can also be precomputed, thus avoiding the real-time matrix inversion.

Taking into account that $\bold{\hat{\tilde{h}}_p}_{\text{LS}}= \text{diag}\left\{\bold{\tilde{x}}_{p}\right\}^{-1}\bold{\tilde{y}}_{p}$, the LS-BEM channel estimate can be represented as
\begin{equation}
\bold{\hat{\tilde{h}}_p}_{\text{LS-BEM}} = (\bold{B}^H\bold{B})^{-1}\bold{B}^H\bold{\hat{\tilde{h}}_p}_{\text{LS}},
\end{equation}
where $(\bold{B}^H\bold{B})^{-1}\bold{B}^H$ can also be precomputed to reduce the real-time computation.

The LS-BEM algorithm has computational complexity $\mathcal{O} (N_p N_a)$, where $N_a$ is the length of the vector $\bold{a}$ which is smaller than $N_p$. However, the LS-BEM channel estimates are sensitive to the additive noise.

\subsection{LMMSE-BEM}

To improve the performance, we take the level of the noise (i.e. the noise variance) into account \cite{Tang}. This is usually done by applying the Bayesian Gauss-Markov theorem~\cite{Kay} to minimize a cost function $J(\bold{\epsilon})= E\{|\bold{\epsilon}|^2\}$, with $\epsilon=\bold{a}-\hat{\bold{a}}$, and by obtaining the LMMSE estimator for $\bold{a}$ as given by
\begin{equation} \label{eq:lmmse3}
 \bold{\hat{a}}=  \left(\bold{B}^H\bold{B} + \sigma_{\tilde{n}}^2\bold{R}_{a}^{-1}\right)^{-1}\bold{B}^H\bold{\tilde{X}}^H\bold{\tilde{y}}_{p},
\end{equation}
where $\bold{R}_{a}= E\{\bold{a}\bold{a}^H\}$ is the covariance matrix of the expansion coefficients.

The covariance matrix of the expansion coefficients can be obtained from the channel covariance, assuming that $\tilde{\bold{h}}_{p} = \bold{B}\bold{a}$, i.e., the modeling error is negligible.
We can find $\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}} = E\{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}^H\}$. Then, we write $\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}} = E\{\bold{B}\bold{a}\bold{a}^H\bold{B}^H\} = \bold{B}E\{\bold{a}\bold{a}^H\}\bold{B}^H = \bold{B}\bold{R}_{a}\bold{B}^H$.

By multiplying from the left by $(\bold{B}^H\bold{B})^{-1}\bold{B}^H$ and from the right by $\bold{B}(\bold{B}^H\bold{B})^{-1}$, and also taking into account that $(\bold{B}^H\bold{B})^{-1}\bold{B}^H\bold{B} = \bold{I}$, we obtain
\begin{equation} \label{eq:varch}
\bold{R}_{a} =(\bold{B}^H\bold{B})^{-1}\bold{B}^H\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}}\bold{B}(\bold{B}^H\bold{B})^{-1}.
\end{equation}

The LMMSE-BEM algorithm has computational complexity $\mathcal{O} (N_p N_a^2)$, but the problem to find the channel PDP is still present.

\subsection{BEM designs}

Many traditional BEM designs have been used to model the channels. In this work, we use the Complex Exponentials (CE-BEM)\cite{Tang} and Legendre Polynomials (LP-BEM)\cite{Hrycak}.

The CE-BEM is described by \cite{Tang}
\begin{equation}  \label{eq:b1}
\bold{B}_{\text{CE}}(q,t)= e^{\frac{-j2\pi(p_s.q) (t-1)}{K}},
\end{equation}
where $\left(\frac{p_s.q}{K}\right)$'s are the normalized frequencies in the pilot positions,  $p_s$ is the pilot separation and $t$'s are the basis function index.

The LP-BEM is described by \cite{Hrycak}
\begin{equation}  \label{eq:b2}
\bold{B}_{\text{LP}}(q,t)= \frac{1}{2^t t!} \frac{d^t}{dq^t}[(q^2-1)^t],
\end{equation}
where $\left(\frac{d^t}{dq^t}\right)$'s are the $t$'s derivatives
in $q$.

\section{
Covariance Matrix Approach and Aproximated-LMMSE-BEM Algorithm}
\label{ApCovMat}
As known in the literature, the LMMSE regularization requires the knowledge of the channel covariance matrix, which in turn is defined by the channel PDP. It is difficult to accurately estimate it in practice. In order to overcome this problem, we develop a BEM-type channel estimator for GFDM and other multicarrier systems that considers
that the channel PDP is constant $\frac{1}{L}$ and that the coefficients are evenly distributed along the cyclic prefix (CP), whose length is $L$.

Therefore, we obtain the covariance matrix for the channel approximation as
\begin{equation} \label{eq:acovmat}
\tilde{\bold{R}}_{\bold{h}'\bold{h}'} (n,p)= \sum \limits_{l=0}^{L-1} \left(\frac{1}{L}\right)e^{\frac{-j2\pi p_s(n-p)(l)}{K}},
\end{equation}
where $n,p = 1, 2,\hdots, N_p$ indicates the position of the pilot subcarrier in the symbol of length $K$. This matrix computation is done once and offline.

The Approximated LMMSE (aLMMSE) is described analogously to (\ref{eq:lmmse3}), through replacing $\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}}$ by $\tilde{\bold{R}}_{\bold{h}'\bold{h}'} $ in (\ref{eq:varch}).

The aLMMSE-BEM algorithm is detailed in Algorithm 2.

\begin{algorithm}[]
\caption{aLMMSE-BEM}
\begin{algorithmic}[2]
\STATE{$q=1:N_p$}
\STATE{$t=1:$ number~of~basis~functions}

\STATE{Generate $\bold{B}$ with (\ref{eq:b1}) or (\ref{eq:b2})}
    \FOR{$n=1:N_p$}
      \FOR{$p=1:N_p$}
        \FOR{$l=1:L$}
           \STATE{Generate $\tilde{\bold{R}}_{\bold{h}'\bold{h}'} (n,p)$ with (\ref{eq:acovmat})}
       \ENDFOR
     \ENDFOR
    \ENDFOR
\STATE{Calculate (\ref{eq:varch}) by replacing $\bold{R}_{\tilde{\bold{h}}_{p}\tilde{\bold{h}}_{p}}$ by $\tilde{\bold{R}}_{\bold{h}'\bold{h}'} $}
\STATE{Calculate (\ref{eq:lmmse3}) with the new (\ref{eq:varch})}
\STATE{Calculate $\bold{\hat{\tilde{h}}} = \bold{B}\bold{\hat{a}}$ with the vector $\bold{\hat{a}}$ obtained from (\ref{eq:lmmse3})}
\end{algorithmic}
\end{algorithm}

config: arxiv


TITLE: $\int_{0}^{1} \frac{ \ln^n x}{1+x^{m}} d x$ =?
QUESTION [4 upvotes]: Inspired by the post, I try to generalize the result to
$$
I(m, n)=\int_{0}^{1} \frac{ \ln ^nx}{1+x^{m}} d x
$$
Instead, I first investigate its partner integral $$
\begin{aligned}
I(a): &=\int_{0}^{1} \frac{x^{a}}{1+x^{m}} d x . \\
&=\sum_{k=0}^{\infty}(-1)^{k} \int_{0}^{1} x^{a} \cdot x^{m k} d x \\
&=\sum_{k=0}^{\infty} \frac{(-1)^{k}}{a+m k+1}
\end{aligned}
$$
Differentiating $I(a)$ w.r.t. $a$ by $n$ times at $a=0$ yields $$
\begin{aligned}
I(m, n)&=\left.\frac{{\partial}^{n}}{{\partial }a^{n}} I(a)\right|_{a=0}\\ &=\left.(-1)^{n} n ! \sum_{k=0}^{\infty} \frac{(-1)^{k}}{(a+m k+1)^{n+1}}\right|_{a=0} \\
&=(-1)^{n} n ! \sum_{k=0}^{\infty} \frac{(-1)^{k}}{(m k+1)^{n+1}} \\
&= (-1)^{n} n !\frac{(-1)^{n} n !}{m^{n+1}} \sum_{k=0}^{\infty} \frac{(-1)^{k}}{\left(k+\frac{1}{m}\right)^{n+1}} \\ &= \sum_{k=0}^{\infty} \frac{1}{\left(2 k+\frac{1}{m}\right)^{n+1}}-\sum_{k=0}^{\infty} \frac{1}{\left(2 k+1+\frac{1}{m}\right)^{n+1}}\\& = (-1)^{n} n !\frac{1}{2^{n+1}}\left[\sum_{k=0}^{\infty} \frac{1}{\left(k+\frac{1}{2 m}\right)^{n+1}}-\sum_{k=0}^{\infty} \frac{1}{\left(k+\frac{1}{2}\left(1+\frac{1}{m}\right)\right)^{n+1}}\right]\\
&=\frac{(-1)^{n} n !}{m^{n+1} 2^{n+1}}\left[ \zeta \left(n+1, \frac{1}{2 m}\right)-\zeta \left( n + 1 , \frac { 1 } { 2 } ( 1 + \frac { 1 } { m } ) \right)\right]
\end{aligned}
$$
My Question
I don’t know how to simplify the last answer. Would you help?

REPLY [3 votes]: Take the $n$th derivative of both sides of
$$\int_0^1\frac{y^{s}}{1+y}\mathrm{d}y=\frac12\left[\psi\left(\frac{s+2}{2}\right)-\psi\left(\frac{s+1}{2}\right)\right],\quad \mathfrak{R}(s)>-1$$
w.r.t $s$, we get
$$\int_0^1\frac{y^{s}\ln^n(y)}{1+y}\mathrm{d}y=2^{-n-1}\left[\psi^{(n)}\left(\frac{s+2}{2}\right)-\psi^{(n)}\left(\frac{s+1}{2}\right)\right]$$
let $y=x^m$
$$\int_0^1\frac{x^{sm+m-1}\ln^n(x)}{1+x^m}\mathrm{d}x=(2m)^{-n-1}\left[\psi^{(n)}\left(\frac{s+2}{2}\right)-\psi^{(n)}\left(\frac{s+1}{2}\right)\right]$$
Finally set $s=\frac{1-m}{m}$
$$\int_0^1\frac{\ln^n(x)}{1+x^m}\mathrm{d}x=(2m)^{-n-1}\left[\psi^{(n)}\left(\frac{1+m}{2m}\right)-\psi^{(n)}\left(\frac{1}{2m}\right)\right],\quad n\in \mathbb{Z}^{+},\quad m>0.$$


text:
\begin{document}

\setlength{\abovedisplayskip}{0pt}
\setlength{\belowdisplayskip}{0pt}
\setlength{\abovedisplayshortskip}{0pt}
\setlength{\belowdisplayshortskip}{0pt}

\begin{frontmatter}

\title{Nonsmooth stabilization and its computational aspects}

\author[TUC]{Pavel Osinenko, Patrick Schmidt, Stefan Streif}

\address[TUC]{Technische Universit\"at Chemnitz, Automatic Control and System Dynamics Laboratory, Germany (e-mail: p.osinenko@gmail.com, \{patrick.schmidt, stefan.streif\}@etit.tu-chemnitz.de)}

\begin{abstract}
This work has the goal of briefly surveying some key stabilization techniques for general nonlinear systems, for which, as it is well known, a smooth control Lyapunov function may fail to exist.
A general overview of the situation with smooth and nonsmooth stabilization is provided, followed by a concise summary of basic tools and techniques, including general stabilization, sliding-mode control and nonsmooth backstepping.
Their presentation is accompanied with examples.
The survey is concluded with some remarks on computational aspects related to determination of sampling times and control actions.
\end{abstract}

\begin{keyword}
Stabilization, nonsmooth analysis, control Lyapunov function
\end{keyword}

\end{frontmatter}

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

Nonsmooth tools play central role in nonlinear system stabilization theory as the smooth ones are prone to limitations posed by the celebrated work of \citet{Brockett1983-stabilization}.
He demonstrated that even such simple systems, as a three-wheel robot, do not admit a smooth control Lyapunov function (CLF).
A particular consequence of this fact is that there can be no continuous control law, which depends only on the system's state and which parks the robot into the desired position.
The alternatives are either to consider a \emph{time-varying} (or dynamical) control law or to give up the continuity condition.
The former approach received great attention in the 80s and 90s \citep{Aeyels1985-stabilization,Kawski1989-stabilization-plane,Coron1991-smooth-stabilization,Samson1991-stabilization,Pomet1992-dyn-stabilization,Coron1992-stabilization,Coron1993-stab-time-var,Coron1994-stabilization,Coron1995-stabilization,Khaneja1999-dynamic-feedback,Morin1999-design}.
Although the design of time-varying control laws may happen to be somewhat involved, the great advantage of this approach is that it requires the usual analysis tools and the closed-loop system trajectory is a classical Carathéodory solution which enjoys uniqueness properties.
Contrary to this approach, if the continuity condition of the control law is omitted, care must be taken when defining what a system trajectory actually is.
That is where alternative solutions to the respective initial value problems with discontinuous right-hand side come into play.
One of the most well-known is the Filippov solution \citep{Filippov2013-discont-dyn-sys}.
In brief, it is an absolutely continuous function that satisfies the said initial value problem with the right-hand side interpreted in the sense of a \emph{differential inclusion} of the kind $\dot x \in F(x,t)$, where $F$ is a \emph{set-valued map}.
Sliding-mode control (SMC) makes wide use of Filippov solutions \citep{Slotine1991-nonlin-ctrl,Young1996-slid-mode,Fridman1996-slid-mode,Fridland1999-slid-mode,Perruquetti2002-slid-mode}.
A particular drawback of this kind of a solution in the context of stabilization is that the very same Brockett's conditions apply just as if one were to limit to the classical solutions \citep{Coron1994-stabilization,Ryan1994-stabilization}.
Further solutions include the ones in the sense of Hermes, Krasovskii, Sentis etc. whose good overview was done by \citet{Cortes2008-discont-dyn-sys}.

Of particular interest for this survey is the setting of \emph{sample-and-hold} (S\&H) solutions which are very simple to interpret.
For the original nonlinear system
\begin{equation} \label{eqn:sys}
	\tag{Sys}
	\dot x  = f(x, u), x \in \R^n, u \in \R^m
\end{equation}
if one designs a discontinuous static control law $\kappa: \R^n \ra \R^m$, instead of applying $\kappa$ literally, one can consider its \emph{digital implementation} with a sampling time $\delta > 0$ in the form
\begin{equation} \label{eqn:sys-SH}
	\tag{Sys-SH}
	\begin{aligned}
		& \dot x = f(x, u_k), \\
		& t \in [ k \delta, (k + 1) \delta], u_k \equiv \kappa(x (k \delta)), k \in \N.
	\end{aligned}
\end{equation}
In this case, if $f$ is locally Lipschitz \wrt $x$, classical trajectories of \eqref{eqn:sys-SH} exist (at least locally) and are unique, since $f$ is also measurable in $t$, precisely due to the sample nature of the control.
The only property one has to give up in general when stabilizing \eqref{eqn:sys} in the S\&H mode, is asymptotic stabilizability of the closed loop.
Instead, one has \emph{practical stabilizability} in the sense of
\begin{dfn}[Practical stabilizability]
	\label{dfn:pract-stab-approx-opt}
	A control law $\kappa$ is said to practically stabilize \eqref{eqn:sys} in the sense of S\&H \ie \eqref{eqn:sys-SH} if, for any data $R > r > 0$, there is a sufficiently small $\delta > 0$ such that any closed-loop trajectory $x(t), t \ge 0, x(0) \in \ball_R$ is bounded and enters and stays in the ball $\ball_r$ within a time $T$ depending uniformly on $R, r$.
\end{dfn}
So, by properly selecting the sampling time, one can achieve any desirable precision of stabilization.
Thus, S\&H scenario of stabilization can be justified from both the implementation and usefulness sides.

\begin{table}[h]
	\caption{Overview of stabilization methods}
	\begin{tabular}{|L{0.2\columnwidth}|L{0.2\columnwidth}|L{0.2\columnwidth}|L{0.2\columnwidth}|}
		\hline
		& \textbf{Cont., static} & \textbf{Cont., time-varying} & \textbf{Discont., static}
		\tabularnewline
		\hline
		pros & Classical trajectories, rel. simple design & Classical trajectories & Rel. simple design
		\tabularnewline
		\hline
		cons & Application limited & Sophisticated design & Nonstandard trajectories (however, S\&H)
		\tabularnewline
		\hline
	\end{tabular}
	\label{tab:stabilization-methods}
\end{table}

A brief summary of the discussed stabilization methods can be found in Table \ref{tab:stabilization-methods}.
The focus of this survey is set to discontinuous static control laws in their S\&H realization due to the aforementioned usefulness and practicability.
In general, as said earlier, a nonlinear system does not admit a smooth CLF, so one has to widen the perspective to include nonsmooth tools, such as nonsmooth CLFs, generalized derivatives, subgradients etc., whose brief overview is given in the next section.
Particular stabilization techniques are surveyed in Section \ref{sec:stabilization}.
In Section \ref{sec:backstepping}, we overview a very powerful technique of nonsmooth backstepping with application to dynamically actuated three-wheel robot and dynamical Artstein's circles.
A brief survey of SMC in the S\&H setting is given in Section \ref{sec:SMC}.
The work is concluded with a discussion on computational aspects of stabilization techniques and a related case study in Section \ref{sec:comp-aspects}.

In the following, $\nrm x$ describes the Euclidean norm of $x$ and $\cco (\X)$ is defined as the closure of the convex hull of a set $\X$.
Furthermore, $\ball_R(x)$ denotes a ball with radius $R$ at $x$ \ie $\ball_R(x) := \{ x \in \R^n: \nrm x \leq R \}$ and $\ball_R$ means the same with $x = 0$.
Finally, $\rightrightarrows$ denotes a set-valued mapping.

\section{Basic nonsmooth tools} \label{sec:nonsmooth-basics}

Perhaps, the most central difference between the machinery of continuous and discontinuous stabilization lies at the level of a CLF.
As mentioned earlier, a general nonlinear system does not admit a smooth CLF, so one is forced to consider nonsmooth alternatives.
For this sake, let us first consider
\begin{dfn}[Lower directional generalized derivative]
	\label{dfn:dini-der}
	\textcolor{white}{For} For a locally Lipschitz function $V: \R^n \ra \R$, the function $\D_{\bullet} V(\bullet): \R^n \times \R^n \ra \R$ defined as
	\begin{align} \label{eqn:Dini-der}
		\begin{split}
			\D_{\vartheta} V(x) & \triangleq \liminf_{\mu \ra 0^+} \frac{ V(x + \mu \vartheta) - V(x)}{\mu}.
			\end{split}
		\tag{Der}
	\end{align}
	is called \textit{lower directional generalized derivative} (further, just ``LDGD'').
\end{dfn}

\begin{exm}
	\label{exm:Dini-der}
	Consider a function
	\[
		V(x) = \begin{cases}
			x^2 \sin \frac 1 x, & x < 0, \\
			0, & x = 0, \\
			2x^2 \sin \frac 1 x, & x > 0,
		\end{cases}
	\]
	LDGD at zero along one is $-2$.
	\begin{figure}[H]
		\centering \includegraphics{gfx/Dini-der-exm}
		\caption{Graph of the derivative of $V$}
	\end{figure}
\end{exm}

Using the introduced LDGD, we can now consider one type of a nonsmooth CLF:
\begin{dfn}[CLF in LDGD sense]
	\label{dfn:dini-CLF}
	For the system \eqref{eqn:sys}, a locally Lipschitz, proper and positive-definite function $V: \R^n \ra \R$ is called \textit{CLF in LDGD sense} if there exists a continuous function $w: \R^n \ra \R, x \ne 0 \implies w(x) > 0$ satisfying a \emph{decay condition}: for any compact set $\X \subseteq \R^n$, there exists a compact set $\U_{\X} \subseteq \U$ such that
	\begin{equation}
		\label{eqn:Dini-decay-cond}
			\forall x \in \X \inf_{\vartheta \in \cco(f(x, \U_{\X}))} \; \D_{\vartheta} V(x) \le -w(x).
		\tag{Dec}
	\end{equation}
	The function $w$ is also called \textit{decay function}.
\end{dfn}

The condition \eqref{eqn:Dini-decay-cond} effectively means that $V$ is an \textit{upper minimax solution} of the Hamilton-Jacobi equation
\begin{equation}
	\label{eqn:HJ}
	\tag{HJ}
	\inf_{u \in \U_{\X}}\scal{\nabla V, f(x,u)} + w(x) = 0
\end{equation}
on the respective domain \citep{Camilli2008-CLF-by-Zubov-mtd,Subbotin2013-gen-sol-PDE}.
There exist several techniques of practical stabilization of \eqref{eqn:sys} using an LDGD CLF in the sense of Definition \ref{dfn:dini-CLF}, some of which are reviewed in Section \ref{sec:stabilization}.
Here, it has to be clarified what ``$\nabla V$'' means in the nonsmooth setting.
There are several substitutes for gradients with different contexts, some of which are overviewed here.
First of all, $\nabla V$ for a smooth function $V$ is a unique vector at each point.
If $V$ is nonsmooth, there is no unique vector which describes a descent direction of $V$, so one speaks of a set of those, summarizing them in a \textit{subdifferential}.
Here is the first such subdifferential, which is useful in practical stabilization \citep{Clarke1997-stabilization}:

\begin{dfn}[Proximal subdifferential]
	\label{dfn:prox-subdiff}
	A vector $\zeta \in \R^n$ is called \textit{proximal subgradient} of $V$ at $x \in \R^n$ if there exists a ball $\ball_r(x)$ and $\sigma > 0$ \sut
	\begin{equation}\label{eqn:prox-subgrad}
		\tag{Prox}
		\forall y \in \ball_r(x) \spc V(y) \geq V(x) + \scal{\zeta, y - x} - \sigma \nrm{y - x}^2.
	\end{equation}
	The set of all such vectors is called \textit{proximal subdifferential} and is denoted by $\partial_P V(x)$.

\end{dfn}

It follows straight from the definition of an LDGD, that for any vector $\vartheta$ and any proximal subgradient $\zeta$, it holds that
\begin{equation} \label{eqn:subgradient-zeta-theta}
	\scal{\zeta, \vartheta} \le \D_\vartheta V(x).
\end{equation}
Therefore, a decay condition in the spirit of \eqref{eqn:Dini-decay-cond} can be formulated as
\begin{equation}\label{eqn:prox-CLF}
	\tag{PDec}
	\forall \zeta \in \partial_P V(x) \spc \inf_{u \in \U_\X} \scal{\zeta, f(x,u)} \le - w(x).
\end{equation}
The condition \eqref{eqn:prox-CLF} means that $V$ is a \textit{proximal supersolution} of \eqref{eqn:HJ} \citep{Clarke1995-sys-traj}, or, equivalently, a \textit{viscosity supersolution} thereof \citep{Crandall1983-viscosity-HJB}.

Other subdifferentials exist, for instance:
\begin{dfn}[Limiting subdifferential]
	\label{dfn:limit-subdiff}
	Let $V: \R^n \ra \R$ be a function.
	The set
	\begin{equation} \label{eqn:limit-subdiff}
		\begin{aligned}
			\partial_L V(x) := \{ &\zeta \in \R^n: \zeta = \text{w-lim } \zeta_i, \\
			& \zeta_i \in \partial_P V(x_i), x_i \ra x, V(x_i) \ra V(x)\}
		\end{aligned}
	\end{equation}
	is called \textit{limiting subdifferential} of $V$ at $x$, where w-lim is the weak limit.
\end{dfn}

Here, a sequence $\{ \zeta_i \}_{i = 1, 2, \ldots} \subset \R^n$ is said to converge weakly to $\zeta \in \R^n$, if $\scal{ \zeta_i, \theta } \ra \scal{ \zeta, \theta }$ for all $\theta \in \R^n$.
Such limiting constructions are used \eg in the nonsmooth practical stability analysis of SMC \citep{Clarke2009-slid-mode-stab}.
An overview is presented in Section \ref{sec:SMC}.

\begin{dfn}[Clarke subdifferential]
	\label{dfn:clarke-subdiff}
	Let $\partial_L V(x)$ be the limiting subdifferential of $V$ at $x$.
	Then, the \textit{Clarke subdifferential} $\partial_C V(x)$ is defined as the closed convex hull of \eqref{eqn:limit-subdiff}, i.e.
	\begin{equation}
	\partial_C V(x) = \cco \left( \partial_L V(x) \right).
	\end{equation}
\end{dfn}

\begin{exm}
	Consider $V(x) = \begin{cases}
	x^2, & \text{if } x < 1 \\
	(x-1)^2 + 1, & \text{else}
	\end{cases}$.
	The subdifferentials defined in Definitions \ref{dfn:prox-subdiff}, \ref{dfn:limit-subdiff} and \ref{dfn:clarke-subdiff} are given as
	\begin{align*}
	&\partial_P V(1) = \emptyset, \spc
	\partial_L V(1) = \left\{ \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \end{pmatrix} \right\}, \\
	&\partial_C V(1) = \cco \left( \left\{ \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \end{pmatrix} \right\} \right).
	\end{align*}
		\begin{figure}[H]
		\centering \includegraphics{gfx/Exm-subgradients}
		\caption{Graph of $V$ and the three subdifferentials}
	\end{figure}
\end{exm}

An important property of many practical stabilizing techniques is \textit{semiconcavity}.
In fact, semiconcavity is a ubiquitous property of CLFs \citep{Clarke2011-discont-stabilization, Cannarsa2004-semiconcave}.
\begin{dfn}[Locally semiconcave function]
	\label{dfn:semiconc-fcn}

	A function $V: \X \ra \R$ with $\X \subset \R^n$ being convex, is called \textit{locally semiconcave with linear modulus}, if it is continuous in $\X$ and there exists $C \geq 0$ such that the following inequality holds for all $x, y \in \X$
	\begin{equation}\label{eqn:semiconc}
		\tag{SemiConc}
		V(x) + V(y) - 2 V \left( \frac{x + y}{2} \right) \leq C \nrm{x - y}^2.
	\end{equation}
\end{dfn}

The following theorem states that any locally semiconcave function can be represented as the infimum of a family of $\mathcal C^2$ functions.
The proof can be found in \citep{Cannarsa2004-semiconcave}.

\begin{thm}
	\label{thm:semiconc-minimum-c2}
	Let $V: \X \ra \R, \X \subseteq \R^n$ be a locally semiconcave function with linear modulus according to Definition \ref{dfn:semiconc-fcn}.
	Then, for each compact subset $\set K \subset \X$, there exists a compact set $\Theta \subset \R^{2n}$ and a continuous function $F: \set K \times \Theta \ra \R$, \sut $F(\bullet; \theta)$ is $\mathcal C^2$ for any $\theta \in \Theta$, the gradients $\nabla_x F(\bullet; \theta)$ are equicontinuous, and $V$ can be expressed as
	\begin{equation} \label{eqn:marginal-fcn}
		V(x) = \min_{\theta \in \Theta} F(x; \theta),
	\end{equation}
	for all $x \in \set K$.
	If there exists such a representation of $V$, then $V$ is called \textit{marginal function} \citep{Cannarsa2004-semiconcave}.

\end{thm}

Theorem \ref{thm:semiconc-minimum-c2} will be useful in presenting a particular technique of nonsmooth backstepping in Section \ref{sec:backstepping}.
Marginal functions are also used to define
\begin{dfn}[$F$-disassembled subdifferential]
	\label{dfn:disass-diff}
	\textcolor{white}{Let} Let $V: \X \ra \R, \spc \X \subseteq \R^n$ be a locally semiconcave function and let $F: \set K \times \Theta \ra \R$, where $\set K \subset \R^n$ and $\Theta \subseteq \R^{2n}$ are compact sets.
	The set-valued map $\partial_D^F V: \X \rightrightarrows \R^n$ defined as
	\begin{equation} \label{eqn:disass-diff}
		\partial_D^F V(x) \triangleq \left\{ \pdiff{F(x; \theta)}{x}: \theta \in \argmin_{\theta \in \Theta} F(x;\theta) \right\}
	\end{equation}
	is called \textit{$F$-disassembled subdifferential}.
	A single element of $\partial_D^F V(x)$ is called \textit{$F$-disassembled subgradient}.
\end{dfn}

Such a subdifferential was used in \citep{Cannarsa2004-semiconcave} and was named in \citep{Nakamura2013-asymptotic}.

Using disassembled subdifferentials, the following type of a nonsmooth CLF can be introduced \citep{Nakamura2013-asymptotic}
\begin{dfn}[$F$-disassembled CLF]
	\label{dfn:disassem-clf}
	A proper, positive-definite, locally semiconcave function $V: \R^n \ra \R$ is called \textit{$F$-disassembled CLF} for \eqref{eqn:sys}, if there exists a continuous, positive definite function $w: \R^n \ra \R$ satisfying the following decay condition: for any compact subset $\X \subseteq \R^n$, there exists a compact set $\U_{\X} \subseteq \U$, such that
	\begin{equation}\label{eqn:dissassem-CLF}
		\tag{DisDec}
		\begin{aligned}
			\forall x \in \X \spc \exists \zeta \in \partial_D^F V(x): \spc
			&\min_{u \in \U_{\X}} \scal{\zeta, f(x, u)} \leq -w(x).
		\end{aligned}
	\end{equation}
\end{dfn}

An immediate relation between $F$-disassembled and proximal subdifferentials can be stated in the following lemma.
\begin{lem}
	\label{lem:prox-disass-subgr}
	Let $\X \subseteq \R^n$ be open and let $\Theta \subset \R^{2 n}$ be compact.
	Let $V$ be a locally semiconcave function given as \eqref{eqn:marginal-fcn} for $F$ continuous in $\X \times \Theta$.
	If $\xi$ is a proximal subgradient of $V$ at $x$, then it is also an $F$-disassembled subgradient of $V$ at $x$ \ie
	\begin{equation} \label{eqn:relation-disass-prox-subgr}
		\partial_P V(x) \subseteq \partial_L V(x) \subseteq \partial_D^F V(x).
	\end{equation}
\end{lem}

\begin{pf}
	The first inclusion follows from Definition \ref{dfn:limit-subdiff}.
	The second one is a part of the proof of Theorem 3.4.4 in \citep{Cannarsa2004-semiconcave}.
	\raggedleft \hfill $\blacksquare$
\end{pf}

The following lemma shows the relation between the two different kinds of CLFs.

\begin{lem}
	\label{lem:Prox-CLF-Disass-CLF}
	Consider \eqref{eqn:sys} and a related $F$-disassembled CLF $V: \R^n \ra \R$ with an arbitrary $F$.
	Then, $V$ is also a CLF in LDGD sense.
\end{lem}

\begin{pf}
	Choose a function $F: \X \times \Theta \ra \R$, $F \in \mathcal C^2(\X \times \Theta)$ satisfying \eqref{eqn:marginal-fcn} and let $\partial_D^F V(x)$ be the corresponding $F$-disassembled differential.
	Define for all $x \in \X$
	\begin{align*}
		\Phi (x, \vartheta, \eps; \theta) := \frac{F(x + \eps \vartheta; \theta) - F(x; \theta)}{\eps}
	\end{align*}
	Obviously,
	\begin{equation} \label{eqn:Dini-lim-liminf}
		\begin{aligned}
			&\liminf_{\eps \ra 0^+} \Phi (x, \vartheta, \eps; \theta) \\
			&\le \lim_{\eps \ra 0} \Phi (x, \vartheta, \eps; \theta) = \scal{\pdiff{F(x; \theta)}{x}, \vartheta},
		\end{aligned}
	\end{equation}
	holds for all $x \in \X$, which follows from the definition of $\liminf$.
	Since $\Phi$ is differentiable w.r.t. $x$ and $\theta$, the equality in \eqref{eqn:Dini-lim-liminf} holds for all $\theta \in \Theta$.
	In particular, for $\theta^\star \in \argmin_{\theta \in \Theta} F(x; \theta)$ it holds that
	\begin{equation}
		\begin{aligned}
			\lim_{\eps \ra 0^+} \Phi(x, \vartheta, \eps; \theta^\star)

			= \scal{\left.\pdiff{F(x; \theta)}{x} \right \vert_{\theta = \theta^\star}, \vartheta} = \scal{\zeta, \vartheta}
		\end {aligned}
	\end{equation}
	for some $\zeta \in \partial_D^{F} V(x)$ corresponding to $\theta^\star$ at $x$.
	In turn, one obtains using \eqref{eqn:subgradient-zeta-theta}
	\begin{equation} \label{eqn:Disass-CLF-Dini-CLF}
		\begin{aligned}
			\lim_{\eps \ra 0^+} \Phi (x, \vartheta, \eps; \theta^\star) = \scal{\zeta, \vartheta} \le \D_\vartheta V(x).
		\end {aligned}
	\end{equation}
	Since $\lim_{\eps \ra 0^+} \Phi (x, \vartheta, \eps; \theta^\star)$ is independent of the choice of $F$, \eqref{eqn:dissassem-CLF} $\implies$ \eqref{eqn:Dini-decay-cond} holds with \eqref{eqn:Disass-CLF-Dini-CLF}.
	And since $F$ was chosen arbitrary, \eqref{eqn:Disass-CLF-Dini-CLF} holds for all $F$ that satisfy \eqref{eqn:marginal-fcn}.
	\raggedleft \hfill $\blacksquare$
\end{pf}

Before proceeding to concrete stabilizing techniques, the \textit{inf-convolution} (InfC) should be recalled \citep{Clarke2008-nonsmooth-analys}:
\begin{dfn}[Inf-convolution]
	\label{dfn:InfC}
	Let $V: \X \ra \R, \X \subseteq \R^n$.
	For $\alpha \in (0,1)$, the \textit{inf-convolution} of $V$ at $x$ is defined by
	\begin{equation} \label{eqn:InfC}
		\tag{InfC}
		V_\alpha(x) \triangleq \inf_{y \in \R^n} \left\{ V(y) + \frac{1}{2 \alpha^2} \nrm{y - x}^2 \right\}.
	\end{equation}
\end{dfn}

In classical convex analysis, $V_\alpha$ is known as Yoreau-Mosida regularization of a (convex) function $V$.
Furthermore, if $V$ is a lower semicontinuous function and bounded from below, then $V_\alpha$ is locally Lipschitz and an approximation of $V$ in the sense of $\lim_{\alpha \ra 0} V_\alpha(x) = V(x)$ \citep{Clarke1997-stabilization}.

\section{Stabilization techniques} \label{sec:stabilization}

In the following section, some stabilization techniques are presented.
Some of them are also described in \citep{Braun2017-SH-stabilization-Dini-aim}.

\subsection{Steepest descent}

A steepest descent control law $\kappa: \R^n \ra \U \subseteq \R^m$ at $x \ne 0$ is computed via

\begin{equation} \label{eqn:steepest-descent}
	\kappa(x) \in \argmin_{u \in \U} \D_{f(x, u)} V(x).
\end{equation}

It is shown in \citep{Braun2017-SH-stabilization-Dini-aim} that $\kappa(x)$, computed by \eqref{eqn:steepest-descent} with a semiconcave LDGD CLF $V: \R^n \ra \R_{\geq 0}$, practically asymptotically stabilizes the origin of \eqref{eqn:sys}.
The semiconcavity of $V$ is crucial to guarantee practical stabilizability by steepest descent \citep{Braun2017-SH-stabilization-Dini-aim}.

\subsection{Dini Aiming}

A control law at the state $x \in \R^n$ by Dini Aiming \citep{Kellett2000-Dini-aim} is computed in two steps based on a nondecreasing continuous function $\sigma: \R_{\geq 0} \ra \R_{\geq 0}$.

\begin{itemize}
	\item[1] Identify a direction $\vartheta^\star$ by means of minimizing the LDGD CLF $V$ over a neighborhood of $x$ for a given $r > 0$ \ie
		\begin{equation} \label{eqn:dini-aiming-direction}
			\vartheta^\star \in \argmin_{s \in \bar \ball_r (x)} V(s).
		\end{equation}

	\item[2] Compute an admissible control $\kappa(x) \in \U \cap \bar \ball_{\sigma(\nrm x + r)}$, $\U \cap \bar \ball_{\sigma(\nrm x + r)} = \{ u \in \U: \nrm u \leq \sigma(\nrm x + r) \}$ via
		\begin{equation} \label{eqn:dini-aiming-feeback}
			\kappa(x) = \argmin_{u \in \U \cap \bar \ball_{\sigma(\nrm x + r)} } \frac{ \scal{ x - \vartheta^\star, f(x, u) } }{ \nrm{ x - \vartheta^\star } }
		\end{equation}
\end{itemize}

If the sampling time $\delta$ is chosen small enough (in accordance with $\sigma$ and $r$), the control computed in \eqref{eqn:dini-aiming-feeback} practically asymptotically stabilizes the origin of \eqref{eqn:sys} \citep{Kellett2004-Dini-aim, Braun2017-SH-stabilization-Dini-aim}.
Such a function $\sigma$ always exists, since $V$ is an LDGD CLF.
Furthermore, $\nrm x \ra 0 \implies \nrm u \ra 0$ has to hold.
In contrast to the steepest descent, $V$ does not necessarily need to be semiconcave.

\subsection{Optimization-based control}

In optimization-based control, the given LDGD is directly minimized over the set of admissible constant inputs \citep{Braun2017-SH-stabilization-Dini-aim} \ie
\begin{equation} \label{eqn:optim-based-feedb-optim-1}
	\min_{u \in \U} \int_0^\delta \D_{f(x, u)} V(\varphi(t, x, u)) \spc \mathrm d t
\end{equation}
or, equivalently,
\begin{equation} \label{eqn:optim-based-feedb-optim-2}
	\min_{u \in \U} \abs{ V(\varphi(\delta, x, u)) - V(x)}.
\end{equation}
At every step, $\varphi(\delta, x, u)$, as a solution of \eqref{eqn:sys}, is computed over the sampling time period $[0, \delta]$.
The one step optimization-based control can be defined as
\begin{equation} \label{eqn:optim-based-feedb}
	\kappa_\delta(x) \in \argmin_{u \in \U} V(\varphi(\delta, x, u)).
\end{equation}
This method combines the two steps \eqref{eqn:dini-aiming-direction} and \eqref{eqn:dini-aiming-feeback} in one single optimization problem.
In comparison to other techniques, $\kappa_\delta(x)$ computed in \eqref{eqn:optim-based-feedb}, explicitly depends on $\delta$.

\subsection{Inf-convolution-based stabilization} \label{sub:infC-stab}

In this technique, a minimizer of \eqref{eqn:InfC} \ie
\begin{equation} \label{eqn:InfC-minimizer}
	y_\alpha(x) \in \arginf_{y \in \R^n} \left\{ V(y) + \frac{1}{2 \alpha^2} \nrm{y - x}^2 \right\}
\end{equation}
is computed to define a proximal subgradient $\zeta_\alpha (x) := \frac{x - y_\alpha(x)}{\alpha^2}$ \citep{Clarke1997-stabilization}.
A control law is obtained by
\begin{equation} \label{eqn:InfC-ctrl-law}
	\kappa(x) = \arginf_{u \in \U_{\Y}} \scal{\zeta_\alpha(x), f(y_\alpha(x), u)}
\end{equation}
where $\Y$ is a compact set which contains $y_\alpha(x)$.

Numerical studies with the above described methods can be found in \citep{Braun2017-SH-stabilization-Dini-aim,Osinenko2018-practical-SH}.
The next section discusses specifically the use of $F$-disassembled subdifferentials and CLFs for nonsmooth backstepping.

\section{Nonsmooth backstepping} \label{sec:backstepping}

In this section, a variant of nonsmooth backstepping on the example of three-wheel robot with dynamical actuators and dynamical Artstein's circles is presented based on $F$-disassembled subdifferentials and CLFs.

\subsection{Three-wheel robot with dynamical actuators} \label{sub:ENDI}

A three-wheel robot with dynamical control of the driving and steering torques can be described as follows:
\begin{equation} \label{eqn:ENDI}
	\tag {ENDI}
	\begin{aligned}
		\dot x_1 &= \eta_1, & \dot \eta_1 &= u_1,\\
		\dot x_2 &= \eta_2, & \dot \eta_2 &= u_2, \\
		\dot x_3 &= \eta_1 x_2 - x_1 \eta_2. & \\

	\end{aligned}
\end{equation}
The system \eqref{eqn:ENDI}, $\frac{\mathrm d}{\mathrm d t} \begin{bmatrix} x \\ \eta \end{bmatrix} = f_{\text{ENDI}}(x, \eta, u)$, is also called \textit{extended nonholonomic dynamical integrator} (ENDI) \citep{Abbasi2017-backstepping, Sankaranarayanan2009-switched, Pascoal2002-practical}.
It is essentially the Brockett's \textit{nonholonomic integrator} (NI) with additional integrators before the control inputs.
The former reads, accordingly, as
\begin{equation} \label{eqn:NI}
	\tag {NI}

	\dot x =

	f_{\text{NI}}(x, u) =
	\underbrace{\begin{pmatrix}
		1 \\ 0 \\ -x_2
	\end{pmatrix}}_{=: g_1(x)} u_1 + \underbrace{\begin{pmatrix}
		0 \\ 1 \\ x_1
	\end{pmatrix}}_{=: g_2(x)} u_2.

\end{equation}

The following functions are nonsmooth LDGD CLFs for \eqref{eqn:NI} \citep{Braun2017-SH-stabilization-Dini-aim, Clarke2011-discont-stabilization}:
\begin{align*}
	V_1(x) &= x_1^2 + x_2^2 + 2 x_3^2 - 2 \abs{x_3} \sqrt{x_1^2 + x_2^2}, \\
	V_2(x) &= x_1^2 + x_2^2 + 2 x_3^2 + \abs{x_3} (10 - 2 (\abs{x_1} + \abs{x_2})).
\end{align*}
For the ENDI, nonsmooth backstepping based on \citep{Matsumoto2015-position} may be utilized.
To this end, consider the following function \citep{Kimura2015-asymptotic}:
\begin{equation} \label{eqn:NI-F}
	F(x; \theta) := x_1^4 + x_2^4 + \frac{\abs{x_3}^3}{(x_1 \cos (\theta) + x_2 \sin(\theta) + \sqrt{\abs{x_3}})^2}.
\end{equation}

So, $V(x) := \min_{\theta \in \Theta} F(x; \theta)$ is an $F$-disassembled CLF for \eqref{eqn:NI}.
Choose a minimizer $\theta^\star \in \Theta^\star := \argmin_{\theta \in \Theta} F(x; \theta)$ and compute $\zeta(x; \theta) = \nabla_x F(x; \theta)$.
Then, $\zeta(x; \theta^\star) \in \partial_D^F V(x) = \nabla_x F(x; \theta) \mid_{\theta \in \Theta^\star}$ holds.

Use the Sontag's formula \citep{Nakamura2013-asymptotic} to express, formally,

\begin{equation} \label{eqn:feedback-NI}
	\kappa(x; \theta) = -\begin{pmatrix}
		\scal{\zeta(x; \theta), g_1(x)} \\ \scal{\zeta(x; \theta), g_2(x)}
	\end{pmatrix},
\end{equation}

and plugging in a minimizer $\theta^*$ into the above gives a corresponding stabilizing controller for \eqref{eqn:NI}.

We have, for all $x \in \X \sm \{ 0 \}$ and $\theta \in \Theta$,
\begin{equation} \label{eqn:decay-NI-theta}
	\begin{aligned}
		&\scal{\zeta(x; \theta), f_{\text{NI}}(x, u)} \\
		&=\scal{\zeta(x; \theta), g_1(x)} u_1 + \scal{\zeta(x; \theta), g_2(x)} u_2 \\
		&= -\scal{\zeta(x; \theta), g_1(x)}^2 - \scal{\zeta(x; \theta), g_2(x)}^2 < 0.
	\end{aligned}
\end{equation}

Observe, that for all $\theta \in \Theta$, it holds formally that
\begin{equation}
	\scal{\zeta(x; \theta), G(x) \kappa(x; \theta)} \leq 0
\end{equation}
with $G(x) = \begin{bmatrix} g_1(x) & g_2(x)\end{bmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ -x_2 & x_1 \end{pmatrix}$.

Now, the $F$-disassembled CLF \eqref{eqn:NI-F} of \eqref{eqn:NI} is augmented in the spirit of backstepping as follows
\begin{equation} \label{eqn:ENDI-V}
	V_c(x, \eta) = \min_{\theta \in \Theta} \left\{ F(x; \theta) + \frac{1}{2} \nrm{ \eta - \kappa(x; \theta) }^2 \right\}.
\end{equation}

The $F$-disassembled subdifferential for $\theta_c^\star \in \Theta^\star_c$, where $\Theta^\star_c := \argmin_{\theta \in \Theta} \{ F(x; \theta) + \nicefrac 1 2 \nrm{\eta - \kappa(x; \theta)}^2 \}$

is given by
\begin{equation} \label{eqn:disass-grad-ENDI}
	\begin{aligned}
		\partial_D^F V_c(x, \eta) &= \left.\pdiff{F (x; \theta) + \frac{1}{2} \nrm{ \eta - \kappa(x; \theta) }^2 }{\begin{bmatrix} x \\ \eta \end{bmatrix}} \right\vert_{\theta \in \Theta^\star_c} \\
		&= \left. \begin{bmatrix}
		\pdiff{F(x; \theta)}{x} - \pdiff{\kappa(x; \theta)}{x} (\eta - \kappa(x; \theta)) \\
		\eta - \kappa(x; \theta)
		\end{bmatrix} \right \vert_{\theta \in \Theta^\star_c}.
	\end{aligned}
\end{equation}

Define formally (for a generic $\theta$) $\zeta_c(x, \eta; \theta) =$ \\ $ \nabla_{(x,\eta)} F_c(x, \eta; \theta)$ and $z(x;\theta) := \eta - \kappa(x; \theta)$.
Then, express
\begin{equation} \label{eqn:decay-ENDI-theta}
	\begin{aligned}
		&\scal{\zeta_c(x, \eta; \theta), f_{\text{ENDI}}(x, \eta,u)} \\
		&= \scal{
		\begin{bmatrix}
			\zeta(x; \theta) - \nabla_x \kappa(x; \theta) z \\
			z
		\end{bmatrix},
		\begin{bmatrix}
			G(x) \eta \\
			u
		\end{bmatrix}
		} \\
		&= \scal{\zeta(x; \theta) - \nabla_x \kappa(x; \theta) z, G(x) \eta} + \scal{z, u} \\
		&= \scal{\zeta(x; \theta), G(x) \eta} - \scal{\nabla_x \kappa(x; \theta) z, G(x) \eta} + \scal{z, u} \\
		&= \scal{\zeta(x; \theta), G(x) z} + \scal{\zeta(x; \theta), G(x) \kappa(x; \theta)} \\
		& \ \ \ - \scal{\nabla_x \kappa(x; \theta) z, G(x) \eta} + \scal{z, u} =: S(x, z; \theta). \\
	\end{aligned}
\end{equation}

Let us pick a $\theta^\star_c \in \Theta^\star_c$ and analyze $S(x, z; \theta^\star_c)$.
The following cases are possible.

Case 1: $z(x; \theta^\star_c) = 0$.
Since $z = 0$, $S(x, z; \theta^\star_c)$ reduces to
$\scal{\zeta(x; \theta^\star_c), G(x) \kappa(x; \theta^\star_c)}$.
But, in this case, $\Theta^\star = \Theta^\star_c$ and so the control for \eqref{eqn:NI} is resembled, and $S(x, z; \theta^\star_c) < 0$ accordingly.

Case 2: $z(x; \theta^\star_c) \not = 0$.
Since
\begin{equation} \label{eqn:feedback-ENDI}
	\begin{aligned}
		&S(x, z; \theta^\star_c) \\
		&= \scal{z, G(x)^{\mathrm T} \zeta(x; \theta^\star_c)} + \scal{\zeta(x; \theta^\star_c), G(x) \kappa(x; \theta^\star_c)} \\
		& \ \ \ - \scal{z, \nabla_x \kappa(x; \theta^\star_c)^{\mathrm T} G(x) \eta} + \scal{z, u} \\
	\end{aligned}
\end{equation}
holds, a suitable choice for $u$ is \eg
\begin{equation} \label{eqn:ENDI-ctrl}
	\begin{aligned}
		& u = \kappa_c(x, \eta; \theta^\star_c) := \\
		& \nabla_x \kappa(x; \theta^\star_c)^{\mathrm T} G(x) \eta - G(x)^{\mathrm T} \zeta(x; \theta^\star_c) - K z
	\end{aligned}
\end{equation}
with $K  > 0$.
This yields
\begin{equation}
	\begin{aligned}
		\scal{\zeta(x; \theta^\star_c), G(x) \kappa(x; \theta^\star_c)} - K \nrm z^2 < 0
	\end{aligned}
\end{equation}
for all $(x,\eta) \in \R^5 \sm \{ 0 \}$.

Thus, the origin of \eqref{eqn:ENDI} is practically stabilized by $\kappa_c(x, \eta; \theta^\star_c)$.
The overall control thus amounts to finding a minimizer $\theta^\star_c$ and plugging it into the expression of $\kappa_c$.
Notice that $\kappa_c$ contains the expression for $\kappa$ and the corresponding generic form of the subgradient $\zeta$.

As a final note, even a control of the form $u = - K z$ would suffice (in the sense of ultimate upper boundedness of the state) using the fact that the terms in $S$ may be compensated accordingly (details omitted).

\subsection{Dynamical Artstein's circles} \label{sub:Artstein}

The same method as described in the previous section can be applied to the following system, a dynamic extension of the Artstein's Circle:
\begin{equation} \label{eqn:Artsteins-Circle}
	\tag{AC}
	\begin{aligned}
		\dot x_1 &= (-x_1^2 + x_2^2) x_3, \\
		\dot x_2 &= -2 x_1 x_2 x_3, \\
		\dot x_3 &= u,
	\end{aligned}
\end{equation}
which can be written as
\begin{equation}
	\begin{aligned}
		\dot v & = g(v) w \\
		\dot w &= u
	\end{aligned}
\end{equation}
with $v := \begin{pmatrix} x_1 & x_2 \end{pmatrix}^{\mathrm T} \in \R^2$ and $w := x_3 \in \R$.

A nonsmooth LDGD CLF for $\dot v = g(v) u$ is given as \citep{Braun2018-complete}:
\begin{equation} \label{eqn:CLF-LDGD-Artsteins-Circle}
	V(v) = \sqrt{3 x_1^2 + 4 x_2^2} - \abs{x_1}.
\end{equation}

Referring to \eqref{eqn:NI-F}, the following marginal function gives rise to an $F$-disassembled CLF:
\begin{equation} \label{eqn:CLF-disass-Artsteins-Circle}
	F(v; \theta) = \sqrt 3 x_1 \cos(\theta_1) + 2 x_2 \sin(\theta_1) + x_1 \left( \frac{\theta_2}{\pi} - 1\right).
\end{equation}
Obviously, for \eqref{eqn:CLF-LDGD-Artsteins-Circle} and \eqref{eqn:CLF-disass-Artsteins-Circle}, \eqref{eqn:marginal-fcn} holds with $\Theta = [0, 2 \pi] \times [0, 2 \pi]$.

After computing the $F$-disassembled differential $\partial_D^F V(v) = \nabla_v F(v; \theta) \mid_{\theta \in \Theta^\star}$, choose $\zeta(v; \theta^\star) \in \partial_D^F V(v)$, where $\theta^\star \in \Theta^\star := \argmin_{\theta \in \Theta} F(v; \theta)$.
Then, with $\kappa(v; \theta) = - \scal{ \zeta(v; \theta), g(v) }$ as a feedback, it follows
\begin{equation}
	\scal{ \zeta(v; \theta), g(v) \kappa(v; \theta)} = -\scal{ \zeta(v; \theta), g(v) }^2 < 0
\end{equation}
which holds also for all $\theta \in \Theta$ and all $v \in \R^2 \sm \{ 0 \}$, and practically stabilizes $\dot v = g(v)w$.
Now, \eqref{eqn:CLF-disass-Artsteins-Circle} can be extended via backstepping to
\begin{equation}
	V_c(v, w) = \min_{\theta \in \Theta_c} \left\{ F(v; \theta) + \frac 1 2 \nrm{ w - \kappa(v; \theta) }^2 \right\}.
\end{equation}

Define $z := w - \kappa(v; \theta)$.
Since the disassembled subdifferential is given similarly to \eqref{eqn:disass-grad-ENDI}, one can show, that for $\zeta_c(v, w; \theta^\star) \in \partial_D^F V_C(v, w)$ with $\theta^\star \in \Theta^\star_c$, $\Theta^\star_c$ defined as in Section \ref{sub:ENDI}, the following equality holds:
\begin{equation}
	\begin{aligned}
		&\scal{
		\begin{bmatrix}
			\zeta(v; \theta) - \nabla_v \kappa(v; \theta) z \\
			z
		\end{bmatrix},
		\begin{bmatrix}
			g w \\
			u
		\end{bmatrix}
		} \\
		&= \scal{\zeta(v; \theta) - \nabla_v \kappa(v; \theta) z, g w} + \scal{z, u} \\
		&= \scal{\zeta(v; \theta), g w} - \scal{\nabla_v \kappa(v; \theta) z, g w} + \scal{z, u} \\
		&= \scal{\zeta(v; \theta), g z} + \scal{\zeta(v; \theta), g \kappa(v; \theta)} \\
		& \ \ \ - \scal{\nabla_v \kappa(v; \theta) z, g w} + \scal{z, u}. \\
	\end{aligned}
\end{equation}
The corresponding control law can be derived in a similar manner as in Section \ref{sub:ENDI}.

\section{Sliding mode} \label{sec:SMC}

Consider a system
\begin{equation}
\label{eqn:sys-SMC}
\begin{array}{ll}
\dot x & = f(x, u), x \in \R^n, \\
u & = \kappa(x) + s(x) \sigma(x),
\end{array}
\end{equation}
where $\kappa, s$ are continuous and
\begin{equation}
\label{eqn:switch-ctrl}
\sigma(x) = \begin{cases}
1, & \chi(x) > 0, \\
-1, & \chi(x) < 0
\end{cases}
\end{equation}
describes the discontinuous part of the controller around a sliding surface $\Sigma = \{x: \chi(x) = 0 \}$.
To process \eqref{eqn:sys-SMC}, unlike in the usual case, two Lyapunov-like functions $V_1, V_2$ are employed that satisfy, for some $\bar w_1 > 0$,

\begin{align*}
& \hspace{-30pt} \begin{array}{ll}
& V_1 \text{ and } V_1 + V_2 \text{ are proper}, V_1(x) = 0 \iff x \in \Sigma, \\
& \text{$V_1$ is continuously differentiable on } \R^n \setminus \Sigma \text{ and} \\
& \langle \nabla V_1(x) , f(x, \kappa(x) + s(x)\sigma(x))\rangle \le - \bar w_1, x \in \R^n \setminus \Sigma;
\end{array} \\
& \hspace{-30pt} \begin{array}{ll}
& V_2(0) = 0, w_2(0) = 0, \\
& x \in \Sigma \setminus \{0\} \implies V_2(x) > 0, w_2(x) > 0, \\
& \text{$V_2$ is continuously differentiable on } \R^n \setminus \{0\} \text{ and} \\

& \sup_{v \in \partial_L^f(x) } \langle \nabla V_2(x), v \rangle \le -w_2(x), x \in \Sigma \setminus \{0\},
\end{array}
\end{align*}

where $\partial_L^f(x) := \{ \lim_{i \ra \infty} f(x, \kappa(x) + s(x) \sigma_i) : x_i \ra x, \sigma_i \in \sigma(x_i) \}$ describes the limiting behavior of possible system velocities (cf. the construction in Definition \ref{dfn:limit-subdiff}).
Let $f_1(x) := f(x, \kappa(x) + s(x))$ and $f_2(x) := f(x, \kappa(x) - s(x))$, and $\X_1 := \{ \chi(x) > 0 \}, \X_2 := \{ \chi(x) < 0 \}$.
In general, \eqref{eqn:sys-SMC} in SMC mode is usually treated in the sense of a differential inclusion
\citep{Perruquetti2002-slid-mode,Slotine1991-nonlin-ctrl}:
\begin{equation}
\label{eqn:SMC-DI}
\dot x	\in F(x) = \begin{cases}
f_1(x), & x	\in \X_1,\\
(1 - \alpha)f1(x) + \alpha f_2(x), \alpha \in [0, 1], & x \in \Sigma, \\
f_2(x), & x	\in \X_2.
\end{cases}
\end{equation}

This is a particular example of \textit{Filippov regularization} \citep{Cortes2008-discont-dyn-sys} and, in this case, is upper semi-continuous, has closed and convex images and thus admits a Filippov solution \citep{Zabczyk2009-mathematical}, which is an absolutely continuous function whose derivative satisfies the differential inclusion almost everywhere.
A Filippov solution is an idealized construction describing a perfect sliding mode in this context.
In practice, since the control is usually realized digitally, S\&H analysis of SMC comes in handy \citep{Clarke2009-slid-mode-stab}.
The general idea thereof is to chose a sampling time bound small enough that the system approach a specially chosen vicinity of $\Sigma$ where $V_2$ still retains some of it decay rate $w_2$ (by a continuity argument).
The attraction to this vicinity is in turn ensured by retaining some of the decay rate $\bar w_1$.
Combining the two decay properties together ensures practical stabilization of \eqref{eqn:sys-SMC} in the S\&H mode \citep{Clarke2009-slid-mode-stab}.

\section{Computational aspects} \label{sec:comp-aspects}

The described practical stabilization methods of the previous sections rely on computation of primarily two things: \begin{enumerate*}
	\item allowed sampling time for the desired stabilization precision;
	\item control actions.
\end{enumerate*}
General techniques of Section \ref{sec:stabilization} rely heavily on various optimizations.
It was shown that optimization accuracy greatly influences stabilization precision \citep{Osinenko2018-practical-SH}.
Moreover, the involved LDGD CLF must satisfy certain regularity properties.
In particular, the following:
For all compact sets $\Y, \F \subset \R^n$ and for all $\nu, \chi > 0$ there exist $\tilde \Y \subseteq \Y, \mu \ge 0$ such that:
\begin{enumerate}
	\item For each $\tilde y \in  \tilde \Y, \theta \in \F$ and $\forall \mu' \in (0, \mu]$ it holds that
	\begin{equation} \label{eqn-loc-hom-cond}
	\tag{hom}
	\abs{ \frac{ V(\tilde y + \mu' \theta) - V(\tilde y) }{\mu'} - \D_\theta V(\tilde y) } \le \nu;
	\end{equation}
	\item For each $y \in \Y$ there exists $\tilde y \in \tilde  \Y$ such that
	\begin{equation} \label{eqn-neighb-pt-cond}
	\tag{npt}
	\nrm{y - \tilde y} \le \chi.
	\end{equation}
\end{enumerate}
Under these conditions, one can find bounds on the optimization accuracy of the inf-convolution-based method so as to achieve the desired stabilization precision.
Furthermore, it was shown in \citep{Osinenko2018-constr-SMC} what regularity properties of the involved Lyapunov-like functions have to be satisfied to effectively compute required bounds on the sampling time.
The respective machinery was addressed on the example of SMC.
In contrast to \citep{Clarke2009-slid-mode-stab}, an actual numerical example of practical SMC stabilization for vehicle slip control was shown.
The computed sampling time bounds were satisfactory for the considered application.

Here, we give a short case study that demonstrates the effect of the described computational uncertainty.
Namely, for \eqref{eqn:ENDI} and its CLF \eqref{eqn:ENDI-V}, a control law was computed in the S\&H framework using \eqref{eqn:InfC-minimizer} and \eqref{eqn:InfC-ctrl-law} with optimization accuracy $\eps$ and $\gamma$, respectively.

The initial condition is set to $x_0 = \begin{pmatrix} -1 & 0.5 & 0.01 & 0.05 & 0.075 \end{pmatrix}^{\mathrm T}$ and the set of admissible controls is given as $\U = [-3,3]$.

Furthermore, we set $\alpha$ in \eqref{eqn:InfC} to $\alpha = 0.1$ and the sampling time $\delta$ to $\delta = 0.005$. \linebreak

The influence of the optimization accuracy on the state and CLF behavior can be seen in Fig. \ref{fig:gfx_epsilon} for different values of $\eps$ and $\gamma$, namely $\eps = \gamma \in \{ 10^{-2}, 10^{-3}, 10^{-4}, 10^{-6}, 10^{-8} \}$.
It can be observed that insufficient accuracy ($\eps = \gamma = 10^{-3}$ and $\eps = \gamma = 10^{-2}$) leads instability.
Higher accuracies lead to ever smaller vicinities of the origin that the state converges into.
This clearly demonstrates that computational uncertainty must be taken into account in practical stabilization.

\begin{figure}
\centering
  \begin{tikzpicture}
	\coordinate (1) at (0,0);
  	\node [right=of 1, scale = 0.8] {\includegraphics[width = 0.5\textwidth]{gfx/CLF_eps_legend}};
  	\coordinate (2) at (1.875,1.75);
  	\node [right=of 2, scale = 0.8] {\includegraphics[width = 0.5\textwidth]{gfx/Norm_eps}};
	\end{tikzpicture}
	\caption{Practical stabilization of the three-wheel robot with dynamic actuators under different computational accuracy.}
	\label{fig:gfx_epsilon}
\end{figure}

\section{Conclusion}

This work surveyed, in a brief form, some of the key modern nonsmooth stabilization tools and techniques.
These include general practical stabilization, sliding-mode control, nonsmooth backstepping.
Examples were provided.
In addition, this work briefly discussed some computational aspects of practical stabilization one should be concerned about in applications.

\bibliography{bib/constructing-LFs,bib/discont-DE,bib/dyn-sys,bib/MPC,bib/non-smooth-analysis,bib/opt-ctrl,bib/sliding-mode,bib/stability,bib/stabilization,bib/nonlin-ctrl,bib/semiconcave,bib/PDE,bib/Osinenko,bib/Filippov-sol}

\end{document}