# A new count model based on Poisson-Transmuted Geometric convolution

Anupama Nandi, Subrata Chakraborty, Aniket Biswas
Department of Statistics
Dibrugarh University

###### Abstract

A novel over-dispersed discrete distribution, namely the PoiTG distribution is
derived by the convolution of a Poisson variate and an independently
distributed transmuted geometric random variable. This distribution
generalizes the geometric, transmuted geometric, and PoiG distributions.
Various important statistical properties of this count model, such as the
probability generating function, the moment generating function, the moments,
the survival function, and the hazard rate function are investigated.
Stochastic ordering for the proposed model are also studied in details. The
maximum likelihood estimators of the parameters are obtained using general
optimization approach and the EM algorithm approach. It is envisaged that the
proposed distribution may prove to be useful for the practitioners for
modelling over-dispersed count data compared to its closest competitors.

Keywords Transmuted geometric distribution, Poisson distribution, Conway-
Maxwell Poisson distribution; $BerG$ distribution; $BerPoi$ distribution;
Incomplete gamma function; EM algorithm.

MSC 2010 60E05, 62E15.

##  1 Introduction

Over-dispersion is a phenomenon that is significant in numerous modelling
applications and is encountered more frequently than under-dispersion or equi-
dispersion. The literature offers several count models designed specifically
for over-dispersed data. In light of the continuing research interest in this
area, it is important to incorporate a simple yet adequate model. The
literature extensively covers several distributions and regression models that
can be utilized to model over-dispersed count data (Campbell et al. [1999],
Wang et al. [2001], Moghimbeigi et al. [2008], Rodrigues-Motta et al. [2013],
Hassanzadeh and Kazemi [2016], Wongrin and Bodhisuwan [2017], Wang et al.
[2017], Sarvi et al. [2019], Moqaddasi Amiri et al. [2019], Tapak et al.
[2020], Tüzen et al. [2020], Altun [2020], and Bar-Lev and Ridder [2021] ).
The Poisson distribution $Poi(\lambda)$, geometric distribution $Geo(\theta)$,
and negative binomial distribution $NB(r,\theta)$ are among the most commonly
used distributions in count data modelling. Despite its widespread use, the
Poisson distribution has a limitation due to its equi-dispersion
characteristics. Several alternative distributions have been developed and
used to address this issue, notable among these distributions are the hyper-
Poisson of Bardwell and Crow [1964], generalized Poisson of Jain and Consul
[1971], double-Poisson of proposed by Efron [1986], weighted Poisson
distribution of Del Castillo and Pérez-Casany [1998], weighted generalized
Poisson of Chakraborty [2010], the popular COM-Poisson distribution of Sellers
and Shmueli [2010], and the Mittag-Leffler function distribution of
Chakraborty and Ong [2017]. In addition to the aforementioned distributions,
several extensions of the geometric distribution are available in the
literature to address over-dispersed count data modeling (Chakraborty and
Gupta [2015], Jain and Consul [1971], Philippou et al. [1983], Tripathi et al.
[1987], Makcutek [2008], Gómez-Déniz [2010], and Nekoukhou et al. [2012] among
others).

An another important generalisation of geometric distribution is transmuted
geometric distribution $TGD(q,\alpha)$ proposed by Chakraborty and Bhati
[2016], which is developed by applying the quadratic rank transmutation
techniques of Shaw and Buckley [2007]. Consider a random variable $Y$ follows
the $TGD(q,\alpha)$, its corresponding probability mass function (pmf) is
given by

$\displaystyle P(Y=y)=(1-\alpha)q^{y}(1-q)+\alpha(1-q^{2})q^{2y},\hskip
6.0pty=0,1,2...\,\,.$ (1)

where $0<q<1$ and $-1<\alpha<1$. It exhibits under-dispersion for
$\alpha\in(-1,0)$ and over-dispersion for $\alpha\in(0,1)$.

Despite the existence of numerous discrete distributions in the literature,
there remains a significant opportunity for developing new discrete
distributions that are appropriate for over-dispersion situations and possess
simple structures and clear, concise interpretations. The convolution method
is a powerful and relatively simple technique for constructing new probability
distributions. By convolving two known distributions, one can obtain a new
distribution with different properties and characteristics, which can be
useful for modelling various types of data. The $BerG$ and $BerPoi$
distributions are two recent examples of such distributions that have been
proposed in the literature. Bourguignon and Weiß [2017] have introduced the
$BerG$ distribution by using the convolution of a Bernoulli random variable
and a geometric random variable. The Bernoulli distribution demonstrates
under-dispersion as its mean is lower than its variance. Conversely, the
geometric distribution exhibits over-dispersion due to its variance being
greater than its mean. By combining the Bernoulli and geometric distributions,
the $BerG$ distribution can effectively model over-dispersed, under-dispersed,
and equi-dispersed data. Bourguignon et al. [2022] have introduced the
$BerPoi$ distribution from a similar motivation; this is derived using the
convolution of a Bernoulli random variable and a Poisson random variable and
is specially designed to model under-dispersed count data.

In a similar way, Nandi et al. [2023] have introduced a novel over-dispersed
count model, namely $PoiG$ distribution by convolving Poisson and geometric
distributions. Let $X_{1}$ be a random variable that follows a Poisson
distribution with parameter $\lambda$, i.e., $X_{1}\sim P(\lambda)$, where
$\lambda>0$. Let $X_{2}$ be a random variable that follows a geometric
distribution with parameter $q$, i.e., $X_{2}\sim G(q)$, where $0<q<1$. Then
$X=X_{1}+X_{2}$ has the following probability mass function

$P(X=x)=\dfrac{(1-q){q^{x}}}{\Gamma(x+1)}\exp{\left(\dfrac{\lambda(1-q)}{q}\right)}{\Gamma\left(x+1,\dfrac{\lambda}{q}\right)},\hskip
30.0ptx=0,1,2,...\,\,.$ (2)

In this article, we introduce a novel discrete distribution that is
constructed by utilizing the convolution of two independent count random
variables, one following the Poisson distribution and the other following the
transmuted geometric distribution. We refer to this novel discrete
distribution as $PoiTG$ distribution. The $PoiTG$ distribution is simpler to
understand compared to the COM-Poisson distribution, as it doesn’t require a
complex normalizing constant in its probability mass function. It is also
well-suited for regression modeling due to having closed-form expressions for
both its mean and variance. The $PoiTG$ distribution proposed in this study is
an extension of both the Poisson distribution and the transmuted geometric
distribution, as well as the recently proposed $PoiG$ distribution , which is
based on the convolution of Poisson and geometric distributions.

This paper is organized as follows. In section 2, we introduce the $PoiTG$
distribution. In section 3, we describe its important statistical properties
such as recurrence relation, generating functions, moments, skewness,
kurtosis, dispersion index, reliability properties, and stochastic ordering.
In Section 4, we discuss two techniques of parameter estimation: the maximum
likelihood method using numerical optimization techniques and the expectation-
maximization (EM) algorithm. We conclude the article with a few limitations
and future scopes of the current study.

##  2 The $PoiTG$ distribution

This section presents a new discrete distribution that arises from the
analysis of two separate, independent discrete random variables $Y_{1}$ and
$Y_{2}$. Let $Y_{1}$ follows the Poisson distribution $P(\lambda)$ where
$\lambda>0$ and $Y_{2}$ follows the transmuted geometric distribution
$TGD(q,\alpha)$ Chakraborty and Bhati [2016] where $0<q<1$ and $0<\alpha<1$,
respectively. Both the variables are restricted to non-negative integer
values. By performing a convolution operation on these two random variables, a
new variable can be generated, $Y=Y_{1}+Y_{2}$, where $Y\in N_{0}$. Then,

$\displaystyle\Pr(Y=y)$
$\displaystyle=\sum_{i=0}^{y}\Pr(Y_{1}=i)\Pr(Y_{2}=y-i)$
$\displaystyle=\sum_{i=0}^{y}\dfrac{e^{-\lambda}\lambda^{i}}{i!}((1-\alpha)(1-q)q^{y-i}+\alpha(1-q^{2})q^{2(y-i)})$
$\displaystyle=(1-\alpha)(1-q)q^{y}e^{-\lambda}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda}{q}\right)^{i}+\alpha(1-q^{2})q^{2y}e^{-\lambda}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda}{q^{2}}\right)^{i},$
(3)

The distribution in (2) is named as the $PoiTG$ distribution and we write
$Y\sim PoiTG(\lambda,q,\alpha)$. Thus, the probability mass function (pmf) of
$PoiTG(\lambda,q,\alpha)$ can be written as

$\displaystyle p_{Y}(y)$ $\displaystyle=(1-\alpha)w_{1}(y,\lambda,q)+\alpha
w_{2}(y,\lambda,q),\hskip 30.0pty=0,1,2,...\,\,.$ (4)

In the above expression,

$\displaystyle w_{k}(y,\lambda,q)$
$\displaystyle=(1-q^{k})q^{ky}e^{-\lambda}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda}{q^{k}}\right)^{i}$
$\displaystyle=\dfrac{(1-q^{k})q^{ky}}{\Gamma(y+1)}\exp{\left(\dfrac{\lambda}{q^{k}}-\lambda\right)}\Gamma\left(y+1,\dfrac{\lambda}{q^{k}}\right),\hskip
30.0ptk=1,2\,\,.$

The incomplete gamma function (Abramowitz and Stegun [1964]) in (4) is defined
as $\Gamma(n,x)=\int_{x}^{\infty}t^{(n-1)}e^{-t}dt$ and it can also be
rewritten as

$\Gamma(n,x)=(n-1)!\sum_{k=0}^{n-1}\dfrac{e^{-x}x^{k}}{k!},$

which is valid for positive values of $n$ and any value of $x$. Thus the
incomplete gamma function in $w_{k}(y,\lambda,q)$ can be rewritten as

$\displaystyle\Gamma\left(y+1,\dfrac{\lambda}{q^{k}}\right)$
$\displaystyle=y!{\sum_{i=0}^{y}}\dfrac{1}{i!}\exp{\left(-\dfrac{\lambda}{q^{k}}\right)}{\left(\dfrac{\lambda}{q^{k}}\right)^{i}},$

where $y!=\Gamma(y+1)$ and $i!=\Gamma(i+1)$. The nature of the pmf (2) of
$PoiTG(\lambda,q,\alpha)$ is demonstrated by Figure 1 for various combinations
of $(\lambda,q,\alpha)$. These plots clearly depict that the distribution of
$PoiTG(\lambda,q,\alpha)$ is unimodal and has, at most, one exponential tail.

The cumulative distribution function (cdf) of $PoiTG$ distribution is

$\displaystyle F_{Y}(y)$ $\displaystyle=\Pr(Y_{1}+Y_{2}\leq y)$
$\displaystyle=\sum_{y_{1}=0}^{y}F_{TG}(y-y_{1})p_{Y}(y_{1})$
$\displaystyle=\sum_{y_{1}=0}^{y}(1-(1-\alpha)q^{y-y_{1}+1}-\alpha
q^{2(y-y_{1}+1)})p_{Y}(y_{1})$
$\displaystyle=\sum_{y_{1}=0}^{y}\dfrac{e^{-\lambda}\lambda^{y_{1}}}{y_{1}!}-(1-\alpha)W_{1}(y,\lambda,q)-\alpha
W_{2}(y,\lambda,q)$ (5)

An explicit expression of (2) is given by

$\displaystyle F_{Y}(y)$
$\displaystyle=\dfrac{\Gamma(y+1,\lambda)}{\Gamma(y+1)}-(1-\alpha)W_{1}(y,\lambda,q)-\alpha
W_{2}(y,\lambda,q),\hskip 25.0pty=0,1,2,...\,\,.$ (6)

In the above expression,

$\displaystyle W_{k}(y,\lambda,q)$
$\displaystyle=q^{k(y+1)}e^{-\lambda}\sum_{y_{1}=0}^{y}\dfrac{1}{y_{1}!}\left(\dfrac{\lambda}{q^{k}}\right)^{y_{1}}$
$\displaystyle=\dfrac{q^{k(y+1)}}{\Gamma(y+1)}\exp{\left(\dfrac{\lambda(1-q^{k})}{q^{k}}\right)}\Gamma\left(y+1,\dfrac{\lambda}{q^{k}}\right),\hskip
30.0ptk=1,2\,\,.$

The nature of the cdf (6) of $PoiTG(\lambda,q,\alpha)$ is demonstrated by
Figure 2 for various combinations of $(\lambda,q,\alpha)$.

Remark 1: As $\lambda\xrightarrow{}0$, the $PoiTG(\lambda,q,\alpha)$
distribution exhibits behavior similar to the $TGD(q,\alpha)$. Similarly, as
$\alpha\xrightarrow{}0$, it behaves akin to the $PoiG(\lambda,1-q)$
distribution. Furthermore, when both $\lambda$ and $\alpha$ $\xrightarrow{}0$,
it resembles the $G(1-q)$ distribution.

Remark 2: $F_{Y}(0)=p_{Y}(0)=(1-q)(1+\alpha q)e^{-\lambda}$ and
$F_{Y}(1)=p_{Y}(0)+p_{Y}(1)=e^{-\lambda}(1-q)((1+\lambda+q)+\alpha
q(q+q^{2}+\lambda).$

Figure 1: Probability mass function of $PoiTG(\lambda,q,\alpha)$ for different
choices of $\lambda$, $q$, & $\alpha$.  Figure 2: Cumulative distribution
function of $PoiTG(\lambda,q,\alpha)$ for different choices of $\lambda$, $q$,
& $\alpha$.

##  3 Properties of the $PoiG$ distribution

In this section, we analyze various statistical properties of the proposed
$PoiTG(\lambda,q,\alpha)$ distribution. The properties examined here include
recurrence relation, probability generating function (pgf), moment generating
function (mgf), characteristic function (cf), cumulant generating function
(cgf), moments, skewness, kurtosis, index of dispersion(ID) and coefficient of
variation(CV). Additionally, we investigate the reliability properties of the
distribution, such as the survival function, hazard rate function and mean
residual life function. We also analyze the stochastic ordering of the
proposed model.

###  3.1 Recurrence relation

Probability recurrence relation helps in finding the subsequent term using the
preceding term. It usually proves to be advantageous in computing the masses
at different values. Note that,

$\displaystyle p_{Y}(y)$
$\displaystyle=(1-\alpha)(1-q)q^{y}e^{-\lambda}s_{y}^{\prime}+\alpha(1-q^{2})q^{2y}e^{-\lambda}s_{y}^{\prime\prime}$

In the above expression,

$s_{y}^{\prime}=\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda}{q}\right)^{i},$
$s_{y}^{\prime\prime}=\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda}{q^{2}}\right)^{i},$
$s_{y+1}^{\prime}=s_{y}^{\prime}+\dfrac{1}{\Gamma(y+2)}\left(\dfrac{\lambda}{q}\right)^{y+1},$

and

$s_{y+1}^{\prime\prime}=s_{y}^{\prime\prime}+\dfrac{1}{\Gamma(y+2)}\left(\dfrac{\lambda}{q^{2}}\right)^{y+1}.$

Now,

$\displaystyle p_{Y}(y+1)$
$\displaystyle=(1-\alpha)(1-q)q^{y+1}e^{-\lambda}s_{y+1}^{\prime}+\alpha(1-q^{2})q^{2(y+1)}e^{-\lambda}s_{y+1}^{\prime\prime}.$
(7)

This is the recurrence formula of the $PoiTG$ distribution. It is easy to
check that

$\displaystyle\dfrac{p_{Y}(y+1)}{p_{Y}(y)}$
$\displaystyle=q\dfrac{s_{y+1}^{\prime}+\dfrac{\alpha(1+q)}{1-\alpha}q^{y+1}s_{y+1}^{\prime\prime}}{s_{y}^{\prime}+\dfrac{\alpha(1+q)}{1-\alpha}q^{y}s_{y}^{\prime\prime}}$
$\displaystyle=q\left(1+\dfrac{\alpha(1+q)q^{y}s_{y}^{\prime\prime}}{(1-\alpha)s_{y}^{\prime}}\right)^{-1}$
(8)

From (8), it is clear that the behaviour of the tail of the distribution
depends on $q$. When $q\xrightarrow{}0$, the tail of the distribution decays
fast, which implies short tail. When $q\xrightarrow{}1$, the tail of the
distribution decays relatively slowly, which implies long tail. The peak of
the distribution curve depends on the value $\alpha$. As the value of $\alpha$
increases, the shape of the probability density function of the
$PoiTG(\lambda,q,\alpha)$ distribution becomes more peaked. These can easily
be verified from Figure 1.

###  3.2 Generating functions

We use the notation $G$ to denote a pgf and use the notation of the
corresponding random variable in the subscript. For $Y_{1}\sim P(\lambda)$ and
$Y_{2}\sim TGD(\alpha,q)$,

$\displaystyle G_{Y_{1}}(s)$ $\displaystyle=e^{\lambda(s-1)}\hskip
10.0pt\textrm{and}\hskip 10.0ptG_{Y_{2}}(s)=\dfrac{(1-q)(1+\alpha
q(1-s)-q^{2}s)}{(1-qs)(1-q^{2}s)},\hskip 5.0pt|q^{2}s|<1.$

Now by using the convolution property of probability generating function we
obtain the pgf of $PoiTG(\lambda,q,\alpha)$ as

$\displaystyle G_{Y}(s)$ $\displaystyle=\dfrac{e^{\lambda(s-1)}(1-q)(1+\alpha
q(1-s)-q^{2}s)}{(1-qs)(1-q^{2}s)}.$ (9)

Similar methods are used to obtain the other generating functions, including
the mgf $M_{Y}(t)$, cf $\phi_{Y}(t)$ and cgf $K_{Y}(t)$. These are given below
as

$\displaystyle M_{Y}(t)=\dfrac{e^{\lambda(e^{t}-1)}(1-q)(1+\alpha
q(1-e^{t})-q^{2}e^{t})}{(1-qe^{t})(1-q^{2}e^{t})},$ (10)
$\displaystyle\phi_{Y}(t)$
$\displaystyle=\dfrac{e^{\lambda(e^{it}-1)}(1-q)(1+\alpha
q(1-e^{it})-q^{2}e^{it})}{(1-qe^{it})(1-q^{2}e^{it})},$ (11)

and

$\displaystyle K_{Y}(t)$
$\displaystyle=\lambda(e^{t}-1)+log\left\\{\dfrac{(1-q)(1+\alpha
q(1-e^{t})-q^{2}e^{t})}{(1-qe^{t})(1-q^{2}e^{t})}\right\\}.$ (12)

###  3.3 Moments and related concepts

By taking the $r^{th}$ derivative of the moment-generating function (mgf) in
equation (10) with respect to $t$ and evaluating it at $t=0$, we can compute
the $r^{th}$ order raw moment. That is,

$\displaystyle E(Y^{r})=M_{Y}^{(r)}(0)=\dfrac{d^{r}}{dt^{r}}[M_{Y}(t)]_{t=0}.$

Let $\mu_{r}^{\prime}$ denote the $r^{th}$ order raw moment, that is
$\mu_{r}^{\prime}=E(Y^{r})$. The explicit expressions of the first four
moments are

$\displaystyle\mu_{1}^{\prime}$
$\displaystyle=\lambda+\dfrac{q(1-\alpha)+q^{2}}{1-q^{2}},$ (13)
$\displaystyle\mu_{2}^{\prime}$
$\displaystyle=\dfrac{1}{(1-q^{2})^{2}}[\lambda(1+\lambda)+q(1-\alpha)(1+2\lambda)+q^{2}(3-2\alpha-2\lambda^{2})$
$\displaystyle\hskip
56.9055pt+q^{3}(1-\alpha)(3-2\lambda)+q^{4}(1-\lambda+\lambda^{2})],$ (14)
$\displaystyle\mu_{3}^{\prime}$
$\displaystyle=\dfrac{1}{(1-q^{2})^{3}}[\lambda(\lambda^{2}+3\lambda+1)+q(1-\alpha)(3\lambda^{2}+6\lambda+1)$
$\displaystyle\hskip
51.21504pt-q^{2}(3\lambda^{3}+6\lambda^{2}-9\lambda+6\alpha(\lambda+1)-7)-2q^{3}(1-\alpha)(3\lambda^{2}-8)$
$\displaystyle\hskip
51.21504pt+q^{4}(3\lambda^{3}+3\lambda^{2}-9\lambda+6\alpha(\lambda-2)+16)+q^{5}(1-\alpha(3\lambda^{2}-6\lambda+7)$
$\displaystyle\hskip 51.21504pt-q^{6}(\lambda^{3}+\lambda-1)],$ (15)
$\displaystyle\text{and}\hskip 6.0pt\mu_{4}^{\prime}$
$\displaystyle=\dfrac{1}{(1-q^{2})^{4}}[\lambda(\lambda^{3}+6\lambda^{2}+7\lambda+1)+q(1-\alpha)(4\lambda^{3}+18\lambda^{2}+14\lambda+1)$
$\displaystyle\hskip
51.21504pt-q^{2}(4\lambda^{4}+20\lambda^{3}-2\lambda^{2}-46\lambda+2\alpha(7+18\lambda+6\lambda^{2})-15)$
$\displaystyle\hskip
51.21504pt-q^{3}(1-\alpha)(12\lambda^{3}+30\lambda^{2}-54\lambda-61)$
$\displaystyle\hskip
51.21504pt+q^{4}(6\lambda^{4}+24\lambda^{3}-24\lambda^{2}+8\alpha(3\lambda^{2}-13)+115)$
$\displaystyle\hskip
51.21504pt+q^{5}(1-\alpha)(12\lambda^{3}+30\lambda^{2}-54\lambda+115)$
$\displaystyle\hskip
51.21504pt-q^{6}(4\lambda^{4}+12\lambda^{3}-14\lambda^{2}+46\lambda+2\alpha(25-18\lambda+6\lambda^{2})-61)$
$\displaystyle\hskip
51.21504pt-q^{7}(1-\alpha)(-15+14\lambda-6\lambda^{2}+4\lambda^{3})+q^{8}(\lambda^{4}+2\lambda^{3}+\lambda^{2}-\lambda+1)].$
(16)

By utilizing the above expressions for the raw moments, we can derive the
explicit formulae for the first four central moments of the random variable
$Y$. The $r^{th}$ central moment, denoted by $\mu_{r}$, is defined as
$\mu_{r}=E(Y-\mu_{1}^{\prime})^{r}$. The explicit expressions for the first
four central moments of the $PoiTG(\lambda,q,\alpha)$ distribution are

$\displaystyle\mu_{1}$ $\displaystyle=0,$ (17) $\displaystyle\mu_{2}$
$\displaystyle=\lambda+\dfrac{q(1-\alpha+q(2+q(1-\alpha)-\alpha^{2}))}{(1-q^{2})^{2}},$
(18) $\displaystyle\mu_{3}$
$\displaystyle=\dfrac{1}{(1-q^{2})^{3}}(\lambda+q(1-\alpha)-q^{2}(3\alpha^{2}+3\lambda-4)-2q^{3}(\alpha^{3}+2\alpha-3)$
$\displaystyle\hskip
133.72786pt-q^{4}(3\alpha^{2}-3\lambda-4)+q^{5}(1-\alpha)-q^{6}\lambda),$ (19)
$\displaystyle\text{and}\hskip 6.0pt\mu_{4}$
$\displaystyle=\dfrac{1}{(1-q^{2})^{4}}((1+q^{8})\lambda(3\lambda+1)+(q+q^{7})(1-\alpha)(6\lambda+1)-(q^{2}+q^{6})$
$\displaystyle\hskip
71.13188pt(\alpha^{2}(6\lambda+4)+6\alpha+12\lambda^{2}-8\lambda-11)+(q^{3}+q^{5})(1-\alpha)$
$\displaystyle\hskip
28.45274pt(6\alpha^{2}+12\alpha-6\lambda+35)+q^{4}(3\alpha^{4}-4\alpha^{2}(3\lambda-7)+12\alpha-2(9\lambda^{2}-9\lambda+25))).$
(20)

The mean and variance of the $PoiTG(\lambda,q,\alpha)$ distribution correspond
to its first raw moment and second central moment, respectively. However, due
to computational complexity and the involvement of large equations, the
expressions for skewness and kurtosis are not provided. Instead, 3-D surface
plots for these measures are presented in Figure 3 and 4.

Figure 3: skewness of $PoiTG(\lambda,q,\alpha)$ for various combination of
$\lambda$, q, $\alpha$.

Figure 3 displays three rows of surface plots representing skewness. Each row
corresponds to a different combination of fixed parameters: $\lambda$, $q$,
and $\alpha$. The first row shows the surface plot of skewness with respect to
$q$ and $\alpha$, while keeping $\lambda$ constant. The second row depicts the
surface plot of skewness in relation to $\lambda$ and $\alpha$, with $q$ held
constant. The third row exhibits the surface plot of skewness with respect to
$\lambda$ and $q$, while maintaining $\alpha$ at a fixed value. These surface
plots visually illustrate how skewness varies when specific parameters are
held constant, offering insights into the relationship between the parameters
and the resulting skewness in the $PoiTG(\lambda,q,\alpha)$ distribution.
Notably, all the surface plots indicate positive skewness for any combination
of $\lambda$, $q$, and $\alpha$.

Regarding the first row’s plot, it demonstrates that as the value of $\lambda$
increases, the skewness of the distribution decreases. In other words, higher
values of $\lambda$ correspond to a reduction in positive skewness. The second
row’s plot shows that the distribution becomes more skewed as the value of $q$
increases. This suggests that higher values of $q$ result in a greater degree
of positive skewness. However, the skewness of the distribution does not seem
to be influenced by the parameter $\alpha$. This indicates that varying the
value of $\alpha$ does not noticeably affect the skewness of the
$PoiTG(\lambda,q,\alpha)$ distribution.

Similarly, surface plots for kurtosis were obtained, with each row using the
same parameter combination as in Figure 3. In the kurtosis plots, the
horizontal surface of each plot is drawn at a value of 3, which never
intersects the kurtosis surface. This suggests a leptokurtic nature of the
$PoiTG(\lambda,q,\alpha)$ distribution. The second row of the kurtosis plot
indicates that the parameter $q$ has a significant influence on the kurtosis.
As $q$ increases, the $PoiTG(\lambda,q,\alpha)$ distribution becomes more
leptokurtic. From the first and third rows, it can be observed that for large
values of $\lambda$ and $\alpha$, the distribution does not show significant
differences in the peakedness of the curve.

Figure 4: kurtosis of $PoiTG(\lambda,q,\alpha)$ for various combination of
$\lambda$, q, $\alpha$.

###  3.4 Index of dispersion(ID) and coefficient of variation(CV)

The index of dispersion (Hoel [1943]), denoted by $I_{Y}$, provides a measure
of the level of dispersion of a distribution. It determines whether a
distribution is suitable for modeling an over-dispersed, under-dispersed or
equi-dispersed dataset. If $I_{Y}$ is greater than one, the distribution of
$Y$ can accommodate over-dispersion, and if $I_{Y}$ is less than one, the
distribution can accommodate under-dispersion. The distribution is said to be
equi-dispersed when $I_{Y}=1$. The dispersion index is given by

$\displaystyle I_{Y}$
$\displaystyle=\dfrac{{\sigma}^{2}}{\bar{x}}=\dfrac{\lambda(1-q^{2})^{2}+q(1-\alpha+q(2+q(1-\alpha)-\alpha^{2}))}{(1-q^{2})(\lambda(1-q^{2})+q(1-\alpha)+q^{2})}$

We present a 3-D Plot 5 of the index of dispersion for the
$PoiTG(\lambda,q,\alpha)$ distribution for various combinations of
$(\lambda,q,\alpha)$ to illustrate the nature of the dispersion. The plot in
the first row shows that the index of dispersion ($I_{Y}$) is slightly
influenced by the choice of $\lambda$ values. As the value of $\lambda$
increases, the corresponding distribution becomes less over-dispersed. The
values of $I_{Y}$ for $\lambda=1.5$ suggest a lower level of over-dispersion
compared to the case where $\lambda=0.5$ and $\lambda=1$.

Figure 5: Index of dispersion of $PoiTG(\lambda,q,\alpha)$ for various
combination of $\lambda$, q, $\alpha$.

The plot in the second row shows that $I_{Y}$ is greatly influenced by the
choice of $q$ values. As the value of $q$ increases, the corresponding
distribution becomes more over-dispersed. The values of $I_{Y}$ for $q=0.6$
and $q=0.8$ suggest a higher level of over-dispersion compared to the case
where $q=0.2$.

The plot in the third row shows that $I_{Y}$ is influenced by the choice of
$\alpha$ values. As the value of $\alpha$ increases, the corresponding
distribution becomes less over-dispersed. The values of $I_{Y}$ decrease as
$\alpha$ increases in the direction of 0.2, 0.6, and 0.8.

The coefficient of variation (CV) is a statistical measure that indicates the
degree of variation in a dataset relative to its mean. A higher value of the
CV indicates that the data has a higher degree of variability. In the context
of probability distributions, a distribution with a higher CV is generally
considered to be better at modeling data with higher variability.

$\displaystyle
CV(Y)=\dfrac{\sqrt{\lambda(1-q^{2})^{2}+q(1-\alpha+q(2+q(1-\alpha)-\alpha^{2}))}}{(\lambda(1-q^{2})+q(1-\alpha)+q^{2})}\times
100\%.$

###  3.5 Reliability properties

The survival function of a discrete random variable $Y$ at $y$ is defined as
the probability of $Y$ assuming values greater than or equal to $y$. The
survival function is a common tool in reliability analysis and is used to
estimate the probability that a system will survive beyond a certain time or
number of events. It is also known as the reliability function because it
provides information about the reliability of a system. The survival function
(sf) of $Y\sim PoiTG(\lambda,q,\alpha)$ is

$\displaystyle S_{Y}(y)$ $\displaystyle=P(Y\geq
y)=1-\dfrac{\Gamma(y,\lambda)}{\Gamma y}+(1-\alpha)W_{1}(y-1,\lambda,q)+\alpha
W_{2}(y-1,\lambda,q).$ (21)

In the above expression,

$\displaystyle W_{k}(y-1,\lambda,q)$ $\displaystyle=\dfrac{q^{ky}}{\Gamma
y}\exp{\left(\dfrac{\lambda(1-q^{k})}{q^{k}}\right)}\Gamma\left(y,\dfrac{\lambda}{q^{k}}\right),\hskip
30.0ptk=1,2\,.$

The hazard rate or failure rate of a discrete random variable $T$ at time
point $t$ is defined as the conditional probability of experiencing a failure
at time $t$, given that the system has survived up to time $t$. The hazard
rate function (hrf) of $Y\sim PoiG(\lambda,q,\alpha)$ is given as

$\displaystyle h_{Y}(y)$ $\displaystyle=\dfrac{P(Y=y)}{S_{Y}(y)}$
$\displaystyle=\frac{(1-q)q^{y}\left(\alpha(1+q)e^{\frac{\lambda}{q^{2}}}q^{y}\Gamma\left(y+1,\frac{\lambda}{q^{2}}\right)+(1-\alpha)e^{\frac{\lambda}{q}}\Gamma\left(y+1,\frac{\lambda}{q}\right)\right)}{ye^{\lambda}(\Gamma
y-\Gamma(y,\lambda))+q^{y}y\left(\alpha
e^{\frac{\lambda}{q^{2}}}q^{y}\Gamma\left(y,\frac{\lambda}{q^{2}}\right)+(1-\alpha)e^{\frac{\lambda}{q}}\Gamma\left(y,\frac{\lambda}{q}\right)\right)}$
(22)

Figure 6 displays the hazard rate function (hrf) for various parameter
settings. In the case of the $PoiTG$ distribution, the behavior of the hrf
varies depending on the values of the parameters. When the value of $\lambda$
is small, the $PoiTG$ distribution demonstrates a constant failure rate. As
$\lambda$ increases, the failure rate shows a slight increase up to a specific
time period, after which it returns to a constant rate. On the other hand, as
the parameter $\alpha$ increases, the failure rate decreases after a certain
time period and remains constant thereafter.

Figure 6: Hazard rate function of $PoiTG(\lambda,q,\alpha)$ for various
combination of $\lambda$, q, $\alpha$.

In the context of reliability studies, the mean residual life refers to the
expected remaining lifetime of a component given that it has already survived
up to a fixed time. The mean residual life (mrl) of a variable $Y\sim
PoiTG(\lambda,q,\alpha)$ is given as

$\displaystyle\mu_{Y}(y)$ $\displaystyle=E(Y-y|Y\geq y)$
$\displaystyle=\frac{\sum_{y=t}^{\infty}\bar{F}(y)}{\bar{F}(t-1)}$
$\displaystyle=\dfrac{\sum_{y=t}^{\infty}\left(1-\dfrac{\Gamma(y,\lambda)}{\Gamma
y}+(1-\alpha)W_{1}(y-1,\lambda,q)+\alpha
W_{2}(y-1,\lambda,q)\right)}{1-\dfrac{\Gamma(t-1,\lambda)}{(\Gamma
t-1)}+(1-\alpha)W_{1}(t-2,\lambda,q)+\alpha W_{2}(t-2,\lambda,q)}$ (23)

In (23), $\bar{F}(y)$ is defined as $1-F(y-1)$, and $\bar{F}(t-1)$ is defined
as $1-F(t-2)$. Similarly, The functions $W_{k}(y-1,\lambda,q)$ and
$W_{k}(t-2,\lambda,q)$ are defined in equation (21), where in the latter case,
$y-1$ is replaced with $t-2$.

###  3.6 Stochastic ordering

The likelihood ratio order is a concept to compare the relative likelihoods of
two random variables. We write $Y\geq_{lr}X$ to mean that $X$ is likelihood
ratio smaller than or equal to $Y$. This is an important statistical property
used to compare the behaviour of two different random variables (Bakouch et
al. [2014]). Establishing the likelihood ratio order between two random
variables can have several practical utilities and implications. Likelihood
ratio order can be used to assess and compare the risk associated with
different events or outcomes. By establishing the order between random
variables representing different risks, we can prioritize interventions or
actions based on their relative likelihoods. In decision theory, establishing
the likelihood ratio order allows us to rank uncertain outcomes based on their
relative likelihoods. Likelihood ratio order can also be used in utility
theory to compare and rank different random variables based on their expected
utilities.

Let $X\sim PoiTG(\lambda_{1},q_{1},\alpha)$ and $Y\sim
PoiTG(\lambda_{2},q_{2},\alpha)$, then $Y$ is said to be smaller than $X$ in
the usual likelihood ratio order, that is $Y\leq_{lr}X$ if
$L(x)=p_{X}(x)/p_{Y}(x)$ is an increasing function in $x$, that is $L(x)\leq
L(x+1)$, which is possible only when $\lambda_{2}\leq\lambda_{1}$ and
$q_{1}\leq q_{2}$ for all $0<\alpha<1$. Note that, for each $x=0,1,2,...$.

$\displaystyle p_{X}(x)$
$\displaystyle=(1-\alpha)(1-q_{1})q_{1}^{y}e^{-\lambda_{1}}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda_{1}}{q_{1}}\right)^{i}+\alpha(1-q_{1}^{2})q_{1}^{2y}e^{-\lambda_{1}}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda_{1}}{q_{1}^{2}}\right)^{i}$
(24) $\displaystyle p_{Y}(x)$
$\displaystyle=(1-\alpha)(1-q_{2})q_{2}^{y}e^{-\lambda_{2}}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda_{2}}{q_{2}}\right)^{i}+\alpha(1-q_{2}^{2})q_{2}^{2y}e^{-\lambda_{2}}\sum_{i=0}^{y}\dfrac{1}{i!}\left(\dfrac{\lambda_{2}}{q_{2}^{2}}\right)^{i}$
(25)

From (24) and (25), we can see that $L(x)\leq L(x+1)$ for all $x=0,1,2..$ when
$\lambda_{2}<\lambda_{1}$ and $q_{1}\leq q_{2}$ for all $0<\alpha<1$. Let
$Y\leq_{st}X$ denote $P(Y\geq x)\leq P(X\geq x)$ for all $x$. This is the
notion of stochastic ordering. Similarly, the hazard rate order $Y\leq_{hr}X$
implies

$\dfrac{p_{X}(x)}{P(X\geq x)}\leq\dfrac{p_{Y}(x)}{P(Y\geq x)},$

for all $x$. The reversed hazard rate order $Y\leq_{rh}X$ implies

$\dfrac{p_{Y}(x)}{P(Y\leq x)}\leq\dfrac{p_{X}(x)}{P(X\leq x)},$

for all $x$. It’s important to note that the likelihood ratio order is a
partial order and not a total order. This means that it does not provide a
strict ranking between all pairs of random variables. Instead, it allows us to
compare the likelihoods of two random variables relative to each other.

From the likelihood ratio order of $X$ and $Y$, the conclusions $Y\leq_{st}X$,
$Y\leq_{hr}X$ and $Y\leq_{rh}X$ are immediate Bakouch et al. [2014].

##  4 Parameter Estimation

Let $\mathbf{Y}=(Y_{1},Y_{2},\ldots,Y_{n})$ denote a random sample of size $n$
drawn from the $PoiTG(\lambda,q,\alpha)$ distribution, and
$\mathbf{y}=(y_{1},y_{2},\ldots,y_{n})$ denote a realization on $\mathbf{Y}$.
The objective of this section is estimate the parameters $\lambda$, $q$ and
$\alpha$ based on the available data $\mathbf{y}$. We present two different
methods of estimation. These methods include maximum likelihood (ML) estimator
obtained through numerical methods and the one obtained via the EM Algorithm.
Additionally, we establish asymptotic confidence intervals for all the
parameters using the maximum likelihood estimates.

###  4.1 Maximum Likelihood Estimator

Using the pmf of $Y\sim PoiTG(\lambda,q,\alpha)$ in (4), the log-likelihood
function of the parameters $\lambda$, $q$ and $\alpha$ can easily be found as

$\displaystyle l(\lambda,q,\alpha;\bm{y})=\log
L(\lambda,q,\alpha;\bm{y})=\sum_{i=1}^{n}\log((1-\alpha)w_{i1}+\alpha
w_{i2}).$ (26)

In the above expression,

$\displaystyle w_{i1}$
$\displaystyle=(1-q)e^{\lambda(1-q)/q}\dfrac{q^{y_{i}}}{y_{i}!}\Gamma\left(y_{i}+1,\dfrac{\lambda}{q}\right),\,\,\,\text{and}$
$\displaystyle w_{i2}$
$\displaystyle=(1-q^{2})e^{\lambda(1-q^{2})/q^{2}}\dfrac{q^{2y_{i}}}{y_{i}!}\Gamma\left(y_{i}+1,\dfrac{\lambda}{q^{2}}\right).$

Here are some notations that will be used in the score functions :

$\displaystyle u_{ij}$
$\displaystyle=\dfrac{e^{-\lambda}q^{j(y_{i}-1)}}{y_{i}!}\left(\dfrac{\lambda}{q^{j}}\right)^{y_{i}}(1-q^{j}),$
$\displaystyle v_{ij}$
$\displaystyle=\dfrac{e^{\lambda(1-q^{j})/q^{j}}q^{jy_{i}}}{y_{i}!}\Gamma\left(y_{i}+1,\dfrac{\lambda}{q^{j}}\right),\,\,\,\text{and}$
$\displaystyle q_{jk}$ $\displaystyle=\dfrac{1-q^{j}}{q^{k}},\hskip
42.67912ptj=1,2\hskip 4.0pt\&\hskip 2.0ptk=0,1,2,3.$

To obtain the score functions, we differentiate (26) with respect to the
parameters $\lambda$, $q$, and $\alpha$.

$\displaystyle s_{1}(\lambda,q,\alpha;\bm{y})$
$\displaystyle=\dfrac{\partial}{\partial\lambda}l(\lambda,q,\alpha;\bm{y})$
$\displaystyle=\sum_{i=1}^{n}\left(\dfrac{\alpha(u_{i1}-u_{i2}+q^{2}q_{22}^{2}v_{i2}-qq_{11}^{2}v_{i1})+u_{i1}+qq_{11}^{2}v_{i1}}{\alpha
q_{20}v_{i2}+(1-\alpha)q_{10}v_{i1}}\right),$ (27) $\displaystyle
s_{2}(\lambda,q,\alpha;\bm{y})=\dfrac{\partial}{\partial
q}l(\lambda,q,\alpha;\bm{y})$
$\displaystyle=\sum_{i=1}^{n}\left(\dfrac{2\alpha(\dfrac{\lambda
u_{i2}}{q}-qv_{i2}-\lambda
q_{23}v_{i2}+q_{21}v_{i2}y_{i})+(1-\alpha)(\dfrac{\lambda
u_{i1}}{q}-v_{i1}-\lambda q_{12}v_{i1}+q_{11}v_{i1}y_{i})}{\alpha
q_{20}v_{i2}+(1-\alpha)q_{10}v_{i1}}\right),$ (28) and $\displaystyle
s_{3}(\lambda,q,\alpha;\bm{y})=\dfrac{\partial}{\partial\alpha}l(\lambda,q,\alpha;\bm{y})=\sum_{i=1}^{n}\left(\dfrac{q_{20}v_{i2}-q_{10}v_{i1}}{\alpha
q_{20}v_{i2}+(1-\alpha)q_{10}v_{i1}}\right).$ (29)

The maximum likelihood estimators are ideally obtained by simultaneously
solving the equations obtained by setting the right hand sides of (27), (28),
and (29) equal to zero. However, in this case, the structural complexity of
the equations prevents obtaining explicit expressions for the maximum
likelihood estimators. Therefore, we resort to numerical optimization
techniques to directly maximize the log-likelihood function with respect to
the parameters. To achieve this, we use constrOptim function of R, which
implements Nelder-Mead method very effectively. In this R function, we put
appropriate constraint matrix and constraint vector as arguments for the
parameter vector $(\lambda,q,\alpha)$. Let $\hat{\lambda}_{ML}$,
$\hat{q}_{ML}$, and $\hat{\alpha}_{ML}$ denote the maximum likelihood
estimates (MLE) of $\lambda$, $q$, and $\alpha$, respectively.

To obtain the information matrix, we need to compute the second-order partial
derivatives of the log-likelihood function with respect to the parameters
$\lambda$, $q$, and $\alpha$. However, the exact expressions for all the
second-order partial derivatives can be cumbersome and difficult to obtain.

In section 1, we discussed the $PoiG(\lambda,q)$ distribution, which was
recently introduced by Nandi et al. [2023]. This distribution is a special
case of the $PoiTG(\lambda,q,\alpha)$ distribution when the parameter $\alpha$
is set to zero. Suppose $x_{1},x_{2},...x_{n}$ is a random sample of size $n$
drawn from the $PoiG(\lambda,q)$ distribution. In that case, we can express
the pmf of $PoiTG(\lambda,q,\alpha)$ as a combination of the pmf of
$PoiG(\lambda,q)$, as shown in (2).

$\displaystyle P_{Y}(y)$ $\displaystyle=(1-\alpha)X_{1}+\alpha X_{2}$

Here, $X_{1}$ is the pmf of $PoiG(\lambda,q)$ and $X_{2}$ is the
$PoiG(\lambda,q^{2})$. The likelihood and log-likelihood function of
$PoiG(\lambda,q)$ distribution is

$\displaystyle L^{\prime}(\lambda,q;y_{i})$
$\displaystyle=(1-q)q^{y_{i}}{y_{i}!}\exp{\left(\dfrac{\lambda(1-q)}{q}\right)}\dfrac{\Gamma\left(y_{i}+1,\dfrac{\lambda}{q}\right)}{\Gamma(y_{i}+1)},$
(30) $\displaystyle\text{and}\,\,\,l^{\prime}(\lambda,q;\bm{y})$
$\displaystyle=n\log(1-q)+n\overline{y}\log
q+\dfrac{n\lambda(1-q)}{q}+\sum_{i=0}^{n}\log\left(\dfrac{\Gamma\left(y_{i}+1,\dfrac{\lambda}{q}\right)}{\Gamma(y_{i}+1)}\right).$
(31)

The score functions are

$\displaystyle
s_{1}^{\prime}(\lambda,q;\bm{y})=\dfrac{\partial}{\partial\lambda}l^{\prime}(\lambda,q;\bm{y})$
$\displaystyle=\frac{n(1-q)}{q}-\sum_{i=1}^{n}\alpha_{1}(y_{i})\beta^{y_{i}},$
$\displaystyle\text{and}\,\,\,\,\,s_{2}^{\prime}(\lambda,q;\bm{y})=\dfrac{\partial}{\partial
q}l^{\prime}(\lambda,q;\bm{y})$
$\displaystyle=-\frac{n}{1-q}-\frac{n(\lambda-\bar{y})}{q}-\frac{n\lambda(1-q)}{q^{2}}+\sum_{i=1}^{n}\lambda\alpha_{2}(y_{i})\beta^{y_{i}}.$

In the above expression, $\beta=\dfrac{\lambda}{q}$ and
$\alpha_{j}(y_{i})=\dfrac{e^{-\beta}}{\Gamma\left(y_{i}+1,\beta\right)}\dfrac{1}{q^{j}}.$
for $j=1,2,3,...\,\,$.

The second-order partial derivative of the log-likelihood of the
$PoiG(\lambda,q)$ are given below.

$\displaystyle\dfrac{\partial^{2}l^{\prime}(\lambda,q;\bm{y})}{\partial\lambda^{2}}$
$\displaystyle=\sum_{i=1}^{n}\left[(\beta^{y_{i}}-y_{i}\beta^{y_{i}-1})\alpha_{2}(y_{i})-\beta^{2y_{i}}\alpha_{1}(y_{i})^{2}\right],$
$\displaystyle\dfrac{\partial^{2}l^{\prime}(\lambda,q;\bm{y})}{\partial\lambda\partial
q}$
$\displaystyle=-\dfrac{n}{q^{2}}-\sum_{i=1}^{n}\left[\lambda(\beta^{y_{i}}-y_{i}\beta^{y_{i}-1})\alpha_{3}(y_{i})-\beta^{y_{i}}\alpha_{2}(y_{i})-\dfrac{\lambda}{q}\beta^{2y_{i}}\alpha_{1}(y_{i})^{2}\right],$
$\displaystyle\text{and}\,\,\,\,\dfrac{\partial^{2}l^{\prime}(\lambda,q;\bm{y})}{\partial
q^{2}}$
$\displaystyle=\dfrac{2n\lambda-n\bar{y}q}{q^{3}}-\dfrac{n}{(1-q)^{2}}+$
$\displaystyle\hskip 12.0pt\sum_{i=1}^{n}\left[((\lambda^{2}-2\lambda
q)\beta^{y_{i}}-\lambda^{2}y_{i}\beta^{y_{i}-1})\alpha_{4}(y_{i})-\lambda^{2}\beta^{2y_{i}}\alpha_{2}(y_{i})^{2}\right].$

The above listed second-order partial derivatives of the log-likelihood
functions for $PoiG(\lambda,q)$ can be utilized to derive the expression for
the second-order partial derivatives of the log-likelihood function for
$PoiTG(\lambda,q,\alpha)$.

$\displaystyle l(\lambda,q,\alpha;\bm{y})$ $\displaystyle=\sum_{i=1}^{n}\log
L(\lambda,q,\alpha;y_{i})=\sum_{i=1}^{n}\log((1-\alpha)P(\lambda,q;y_{i})+\alpha
P(\lambda,q_{*};y_{i})).$

Let consider,

$\displaystyle A_{i1}$
$\displaystyle=(1-\alpha)L^{\prime}(\lambda,q;y_{i})\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial\lambda}+\alpha
L^{\prime}(\lambda,q_{*};y_{i})\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda},$ $\displaystyle A_{i2}$
$\displaystyle=(1-\alpha)L^{\prime}(\lambda,q;y_{i})\bigg{(}\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial\lambda}\bigg{)}^{2}+\dfrac{\partial^{2}l^{\prime}(\lambda,q;y_{i})}{\partial\lambda^{2}}\bigg{)},$
$\displaystyle A_{i3}$ $\displaystyle=\alpha
L^{\prime}(\lambda,q_{*};y_{i})\bigg{(}\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda}\bigg{)}^{2}+\dfrac{\partial^{2}l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda^{2}}\bigg{)},$
$\displaystyle B_{i1}$
$\displaystyle=(1-\alpha)L^{\prime}(\lambda,q;y_{i})\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial q}+2\alpha
qL^{\prime}(\lambda,q_{*};y_{i})\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial q^{*}},$ $\displaystyle B_{i2}$
$\displaystyle=(1-\alpha)L^{\prime}(\lambda,q;y_{i})\bigg{(}\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial
q}\bigg{)}^{2}+\dfrac{\partial^{2}l^{\prime}(\lambda,q;y_{i})}{\partial
q^{2}}\bigg{)},$ $\displaystyle B_{i3}$ $\displaystyle=2\alpha
L^{\prime}(\lambda,q_{*};y_{i})\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial q_{*}}+2q^{2}\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial
q_{*}}\bigg{)}^{2}+2q^{2}\dfrac{\partial^{2}l^{\prime}(\lambda,q_{*};y_{i})}{\partial
q_{*}^{2}}\bigg{)},$ $\displaystyle C_{i1}$
$\displaystyle=(1-\alpha)L^{\prime}(\lambda,q;y_{i})\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q;y_{i}}{\partial\lambda}\hskip 1.0pt\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial
q}+\dfrac{\partial^{2}l^{\prime}(\lambda,q;y_{i})}{\partial\lambda\partial
q}\bigg{)},$ $\displaystyle C_{i2}$ $\displaystyle=2\alpha
qL^{\prime}(\lambda,q_{*};y_{i})\bigg{(}\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda}\hskip 1.0pt\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial
q_{*}}+\dfrac{\partial^{2}l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda\partial
q_{*}}\bigg{)},$ $\displaystyle\text{and}\hskip 28.45274ptC_{i3}$
$\displaystyle=2qL^{\prime}(\lambda,q_{*};y_{i})\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial
q_{*}}-L^{\prime}(\lambda,q;y_{i})\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial q}.$

The second-order partial derivatives of the log-likelihood function for
$PoiTG(\lambda,q,\alpha)$ are given as

$\displaystyle\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda^{2}}$
$\displaystyle=\sum_{i=1}^{n}\dfrac{L(\lambda,q,\alpha;y_{i})(A_{i2}+A_{i3})-A_{i1}^{2}}{L(\lambda,q,\alpha;y_{i})^{2}},$
(32) $\displaystyle\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q^{2}}$
$\displaystyle=\sum_{i=1}^{n}\dfrac{L(\lambda,q,\alpha;y_{i})(B_{i2}+B_{i3})-B_{i1}^{2}}{L(\lambda,q,\alpha;y_{i})^{2}},$
(33)
$\displaystyle\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial
q}$ $\displaystyle=\sum_{i=1}^{n}\dfrac{L(\lambda,q,\alpha
y_{i})(C_{i1}+C_{i2})-A_{i1}B_{i1}}{L(\lambda,q,\alpha;y_{i})^{2}},$ (34)
$\displaystyle\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\alpha^{2}}$
$\displaystyle=-\sum_{i=1}^{n}\dfrac{\left(L^{\prime}(\lambda,q_{*};y_{i})-L^{\prime}(\lambda,q;y_{i})\right)^{2}}{L(\lambda,q,\alpha;)^{2}},$
(35)
$\displaystyle\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial\alpha}$
$\displaystyle=\sum_{i=1}^{n}\dfrac{\left(L^{\prime}(\lambda,q_{*};y_{i})\dfrac{\partial
l^{\prime}(\lambda,q_{*};y_{i})}{\partial\lambda}-L^{\prime}(\lambda,q;y_{i})\dfrac{\partial
l^{\prime}(\lambda,q;y_{i})}{\partial\lambda}\right)}{L(\lambda,q,\alpha;y_{i})^{2}},$
(36) $\displaystyle\text{and}\hskip
28.45274pt\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q\partial\alpha}$
$\displaystyle=\sum_{i=1}^{n}\dfrac{L(\lambda,q,\alpha;y_{i})C_{i3}}{L(\lambda,q,\alpha;y_{i})^{2}}-\sum_{i=1}^{n}\dfrac{\left(L^{\prime}(\lambda,q_{*};y_{i})-L^{\prime}(\lambda,q;y_{i})\right)B_{i1}}{L(\lambda,q,\alpha;y_{i})^{2}}.$
(37)

The Fisher’s information matrix for $(\lambda,q,\alpha)$ is

$\displaystyle
I_{Y}(\lambda,q,\alpha)=\begin{pmatrix}-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda^{2}}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial
q}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial\alpha}\right)\\\
&\\\
-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial
q}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q^{2}}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q\partial\alpha}\right)\\\ &\\\
-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial\alpha}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q\partial\alpha}\right)&-E\left(\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\alpha^{2}}\right)\end{pmatrix}.$

This can be approximated by

$\displaystyle\widehat{I}_{Y}(\lambda,q,\alpha)\approx\begin{pmatrix}-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda^{2}}&\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial
q}&-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial\alpha}\\\
&\\\ -\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial
q}&-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q^{2}}&-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q\partial\alpha}\\\ &\\\
-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\lambda\partial\alpha}&-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial
q\partial\alpha}&-\dfrac{\partial^{2}l(\lambda,q,\alpha;\bm{y})}{\partial\alpha^{2}}\end{pmatrix}_{(\lambda,q,\alpha)=(\hat{\lambda}_{ML},\hat{q}_{ML},\hat{\alpha}_{ML})}.$

For large value of $n$, the maximum likelihood estimators
$\hat{\lambda}_{ML}$, $\hat{q}_{ML}$, and $\hat{\alpha}_{ML}$ are consistent
and asymptotically normal with the mean vector $(0,0,0)$ and the dispersion
matrix $\hat{I}^{-1}=[d_{ij}]_{3\times 3}$, where
$d_{ij}=d(\hat{\lambda}_{ML},\hat{q}_{ML},\hat{\alpha}_{ML};y)$. The
dispersion matrix $\hat{I}^{-1}$ includes the variances of
$\hat{\lambda}_{ML}$, $\hat{q}_{ML}$, and $\hat{\alpha}_{ML}$, which are
represented by $d_{11}$, $d_{22}$, and $d_{33}$, respectively. Consequently,
the asymptotic $(1-\alpha)\times 100\%$ confidence interval for $\lambda$, $q$
and $\alpha$ are given respectively by

$\left(\hat{\lambda}_{ML}-z_{\alpha/2}\sqrt{d_{11}}\hskip 5.0pt,\hskip
5.0pt\hat{\lambda}_{ML}+z_{\alpha/2}\sqrt{d_{11}}\right)\text{,}\left(\hat{q}_{ML}-z_{\alpha/2}\sqrt{d_{22}}\hskip
5.0pt,\hskip 5.0pt\hat{q}_{ML}+z_{\alpha/2}\sqrt{d_{22}}\right)$

and

$\left(\hat{\alpha}_{ML}-z_{\alpha/2}\sqrt{d_{33}}\hskip 5.0pt,\hskip
5.0pt\hat{q}_{ML}+z_{\alpha/2}\sqrt{d_{33}}\right).$

###  4.2 EM algorithm

The expectation-maximization (EM) algorithm is a useful iterative procedure to
compute maximum likelihood estimates of parameters of a mixture distribution.
First, the mixture distribution is represented in such a way that the mixing
coefficient is the parameter of a latent Bernoulli random variable. As pointed
out earlier, $Y=(1-\alpha)X_{1}+\alpha X_{2}$, Where $X_{1}\sim
PoiG(\lambda,q)$ and $X_{2}\sim PoiG(\lambda,q^{2})$. We consider $Z$ to be
another random variable such that, $Z$ is independent of $X_{1}$ and $X_{2}$,
and $Y$ is expressed as the following linear combination of $X_{1}$ and
$X_{2}$.

$\displaystyle Y$ $\displaystyle=(1-Z)X_{1}+ZX_{2}$ (38)

Clearly, a reasonable assumption is $Z\sim Bernoulli(\alpha)$, $0<\alpha<1$.
The estimates obtained by EM algorithm are consistent and unique (Dempster et
al. [1977], Redner and Walker [1984]). Each iteration of EM algorithm consists
of two steps, which are the expectation or E-step and maximization or M-step.
We have $n$ iid copies of $Y$, namely ${Y_{1},Y_{2},...Y_{n}}$. However, in
the incomplete-data framework of EM algorithm, we find that observations of
$Z$ are not available. Thus, the hypothetical complete dataset is
${(Y_{i},Z_{i}):i=1,2,...,n}$. Under the formulation, E-step of each EM
iteration requires the expectation of $(Z|Y;\Psi^{(k)})$, where
$\Psi^{(k)}=(\lambda^{(k)},q^{(k)},\alpha^{(k)})$ is the estimate of
$\Psi=(\lambda,q,\alpha)$ in the $k^{th}$ iteration. Since,

$Z_{i}|Y_{i},\Psi^{(k)}\sim Bernoulli(\alpha_{i}^{(k)}),$
$\displaystyle\alpha_{i}^{(k+1)}=\dfrac{\alpha^{(k)}(1-q^{(k)^{2}})2^{2y_{i}}\hskip
2.0pte^{\dfrac{\lambda^{(k)}(1-q^{(k)^{2}})}{q^{(k)^{2}}}}\Gamma\left(y_{i}+1,\dfrac{\lambda^{(k)}}{q^{(k)^{2}}}\right)}{\alpha^{(k)}(1-q^{(k)^{2}})2^{2y_{i}}\hskip
2.0pte^{\dfrac{\lambda^{(k)}(1-q^{(k)^{2}})}{q^{(k)^{2}}}}\Gamma\left(y_{i}+1,\dfrac{\lambda^{(k)}}{q^{(k)^{2}}}\right)+(1-\alpha^{(k)})(1-q^{k})q^{y_{i}}\hskip
2.0pte^{\dfrac{\lambda^{(k)}(1-q^{(k)})}{q^{(k)}}}\Gamma\left(y_{i}+1,\dfrac{\lambda^{(k)}}{q^{(k)}}\right)}.$
(39)

Now, it is clear that ,

$E[Z_{i}|Y_{i},\Psi^{(k)}]=\alpha_{i}^{k+1}\text{and
}V[Z_{i}|Y_{i},\Psi^{(k)}]=\alpha_{i}^{k+1}(1-\alpha_{i}^{k+1}).$ (40)

We estimate $\alpha$ for the $(k+1)^{th}$ iteration by

$\alpha^{k+1}=\dfrac{1}{n}\sum_{i=1}^{n}\alpha_{i}^{k+1}.$ (41)

Then, we proceed to M-step, where we maximize the log-likelihood of observed
$(y_{1},y_{2},...,y_{n})$ on $(Y_{1},Y_{2},...,Y_{n})$ with respect to
$\lambda$ and $q$ for given $\alpha^{(k+1)}$ and obtain $\lambda^{(k+1)}$ and
$q^{(k+1)}$. That is,

$(\lambda^{(k+1)},q^{(k+1)})=\text{arg}\hskip
3.0pt\underset{\lambda,q}{max}\hskip
3.0ptl(\lambda,q,\alpha^{(k+1)};\underaccent{\tilde}{y}).$

We continue iterating the successive E-step and M-step until
$|\lambda^{(k+1)}-\lambda^{(k)}|<\epsilon$, $|q^{(k+1)}-q^{(k)}|<\epsilon$,
and $|\alpha^{(k+1)}-\alpha^{(k)}|<\epsilon$, simultaneously. We set
$\alpha_{0}=0.5$ and $\epsilon=0.0001$ in the EM algorithm.

##  5 Discussion

This article introduces a novel discrete distribution with three parameters,
which is thoroughly studied. Although the primary focus of this study lies in
theoretical development, its applied aspect is also important. From the
application point of view, the proposed model is easy to use for modeling
over-dispersed data. Despite the existence of alternative count models for
over-dispersion, the interpretability of the proposed model’s parameters could
lead to widespread adoption. The parameters $\lambda$, $q$, and $\alpha$
controls the distribution’s tail behavior, while $q$ and $\alpha$ account for
the over-dispersion in the dataset. Their combined influence provides
flexibility in shaping the distribution. In the case of large $\lambda$, the
distribution maintains a bell-shaped mass distribution. In the complete
version of the paper, comparative real life data analysis and a simulation
experiment to examine the performance of point estimators obtained using the
two above-mentioned estimation methods will be reported.

## References

  * Abramowitz and Stegun [1964] M Abramowitz and I A Stegun.  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , volume 55.  US Government printing office, 1964.
  * Altun [2020] E Altun.  A new generalization of geometric distribution with properties and applications.  _Communications in Statistics-Simulation and Computation_ , 49(3):793–807, 2020.
  * Bakouch et al. [2014] H S Bakouch, M A Jazi, and S Nadarajah.  A new discrete distribution.  _Statistics_ , 48(1):200–240, 2014.
  * Bar-Lev and Ridder [2021] S K. Bar-Lev and Ad Ridder.  Exponential dispersion models for overdispersed zero-inflated count data.  _Communications in Statistics-Simulation and Computation_ , pages 1–19, 2021.
  * Bardwell and Crow [1964] G E Bardwell and E L Crow.  A two parameter family of hyper-Poisson distributions.  _Journal of the American Statistical Association_ , 59:133–141, 1964.
  * Bourguignon and Weiß [2017] M Bourguignon and Christian H Weiß.  An INAR (1) process for modeling count time series with equidispersion, underdispersion and overdispersion.  _Test_ , 26(4):847–868, 2017.
  * Bourguignon et al. [2022] M Bourguignon, D I Gallardo, and R MR Medeiros.  A simple and useful regression model for underdispersed count data based on Bernoulli–Poisson convolution.  _Statistical Papers_ , 63(3):821–848, 2022.
  * Campbell et al. [1999] N L Campbell, L J Young, and G A Capuano.  Analyzing over-dispersed count data in two-way cross-classification problems using generalized linear models.  _Journal of Statistical Computation and Simulation_ , 63(3):263–281, 1999.
  * Chakraborty [2010] S Chakraborty.  On some distributional properties of the family of weighted generalized poisson distribution.  _Communications in Statistics—Theory and Methods_ , 39(15):2767–2788, 2010.
  * Chakraborty and Bhati [2016] S Chakraborty and D Bhati.  Transmuted geometric distribution with applications in modeling and regression analysis of count data.  _SORT-Statistics and Operations Research Transactions_ , pages 153–176, 2016.
  * Chakraborty and Gupta [2015] S Chakraborty and R D Gupta.  Exponentiated geometric distribution: another generalization of geometric distribution.  _Communications in Statistics-Theory and Methods_ , 44(6):1143–1157, 2015.
  * Chakraborty and Ong [2017] S Chakraborty and S H Ong.  Mittag-leffler function distribution-a new generalization of hyper-Poisson distribution.  _Journal of Statistical distributions and applications_ , 4(1):1–17, 2017.
  * Del Castillo and Pérez-Casany [1998] J Del Castillo and M Pérez-Casany.  Weighted poisson distributions for overdispersion and underdispersion situations.  _Annals of the Institute of Statistical Mathematics_ , 50(3):567–585, 1998.
  * Dempster et al. [1977] A P Dempster, N M Laird, and D B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the royal statistical society: series B (methodological)_ , 39(1):1–22, 1977.
  * Efron [1986] B Efron.  Double exponential-families and their use in generalized linear-regression.  _Journal of the American Statistical Association_ , 81:709–721, 1986.
  * Gómez-Déniz [2010] E. Gómez-Déniz.  Another generalization of the geometric distribution.  _Test_ , 19(2):399–415, 2010.  doi: 10.1007/s11749-009-0169-3.
  * Hassanzadeh and Kazemi [2016] F Hassanzadeh and I Kazemi.  Analysis of over-dispersed count data with extra zeros using the Poisson log-skew-normal distribution.  _Journal of Statistical Computation and Simulation_ , 86(13):2644–2662, 2016.
  * Hoel [1943] P G Hoel.  On indices of dispersion.  _The Annals of Mathematical Statistics_ , 14(2):155–162, 1943.
  * Jain and Consul [1971] G C. Jain and P C. Consul.  A generalized negative binomial distribution.  _SIAM Journal on Applied Mathematics_ , 21(4):501–513, 1971.  doi: 10.1137/0121056.
  * Makcutek [2008] J. Makcutek.  A generalization of the geometric distribution and its application in quantitative linguistics.  _Romanian Reports in Physics_ , 60(3):501–509, 2008.
  * Moghimbeigi et al. [2008] A Moghimbeigi, M R Eshraghian, K Mohammad, and B Mcardle.  Multilevel zero-inflated negative binomial regression modeling for over-dispersed count data with extra zeros.  _Journal of Applied Statistics_ , 35(10):1193–1202, 2008.
  * Moqaddasi Amiri et al. [2019] M Moqaddasi Amiri, L Tapak, and J Faradmal.  A mixed-effects least square support vector regression model for three-level count data.  _Journal of Statistical Computation and Simulation_ , 89(15):2801–2812, 2019.
  * Nandi et al. [2023] A Nandi, S Chakraborty, and A Biswas.  A new over-dispersed count model.  _Under Revision_ , 2023.
  * Nekoukhou et al. [2012] V. Nekoukhou, M H. Alamatsaz, and H. Bidram.  A discrete analogue of the generalized exponential distribution.  _Communications in Statistics - Theory and Methods_ , 41(11):2000–2013, 2012.  doi: 10.1080/03610926.2011.555044.
  * Philippou et al. [1983] A N. Philippou, C. Georghiou, and GN. Philippou.  A generalized geometric distribution and some of its properties.  _Statistics and Probability Letters_ , 1(4):171–175, 1983.  doi: 10.1016/0167-7152(83)90025-1.
  * Redner and Walker [1984] R A Redner and H F Walker.  Mixture densities, maximum likelihood and the em algorithm.  _SIAM review_ , 26(2):195–239, 1984.
  * Rodrigues-Motta et al. [2013] M Rodrigues-Motta, H P Pinheiro, E G Martins, Márcio S Araújo, and Sérgio F dos Reis.  Multivariate models for correlated count data.  _Journal of Applied Statistics_ , 40(7):1586–1596, 2013.
  * Sarvi et al. [2019] F Sarvi, A Moghimbeigi, and H Mahjub.  GEE-based zero-inflated generalized Poisson model for clustered over or under-dispersed count data.  _Journal of Statistical Computation and Simulation_ , 89(14):2711–2732, 2019.
  * Sellers and Shmueli [2010] K F Sellers and G Shmueli.  A flexible regression model for count data.  _The Annals of Applied Statistics_ , pages 943–961, 2010.
  * Shaw and Buckley [2007] W T Shaw and I R Buckley.  The alchemy of probability distributions: Beyond gram-charlier & cornish-fisher expansions, and skew-normal or kurtotic-normal distributions.  _Submitted, Feb_ , 7:64, 2007.
  * Tapak et al. [2020] L Tapak, O Hamidi, P Amini, and G Verbeke.  Random effect exponentiated-exponential geometric model for clustered/longitudinal zero-inflated count data.  _Journal of Applied Statistics_ , 47(12):2272–2288, 2020.
  * Tripathi et al. [1987] R C. Tripathi, R C. Gupta, and T J White.  Some generalizations of the geometric distribution.  _Sankhya Ser. B_ , 49(3):218–223, 1987.
  * Tüzen et al. [2020] F Tüzen, S Erbaş, and H Olmuş.  A simulation study for count data models under varying degrees of outliers and zeros.  _Communications in Statistics-Simulation and Computation_ , 49(4):1078–1088, 2020.
  * Wang et al. [2017] S Wang, N G Cadigan, and H P Benoît.  Inference about regression parameters using highly stratified survey count data with over-dispersion and repeated measurements.  _Journal of Applied Statistics_ , 44(6):1013–1030, 2017.
  * Wang et al. [2001] Y Wang, L J Young, and D E Johnson.  A UMPU test for comparing means of two negative binomial distributions.  _Communications in Statistics-Simulation and Computation_ , 30(4):1053–1075, 2001.
  * Wongrin and Bodhisuwan [2017] W Wongrin and W Bodhisuwan.  Generalized Poisson–Lindley linear model for count data.  _Journal of Applied Statistics_ , 44(15):2659–2671, 2017.


# Robust nonlinear state-feedback control of second-order systems

Michael Ruderman and Denis Efimov M. Ruderman is with University of Agder.
Grimstad, Norway. Email<EMAIL_ADDRESS>Efimov is with Inria, Univ.
Lille, CNRS, CRIStAL. Lille, France.This work was supported by RCN grant
number 340782.

###### Abstract

This note proposes a novel nonlinear state feedback controller for perturbed
second-order systems. In analogy to a linear proportional-derivative (PD)
output feedback control, the proposed nonlinear scheme uses the output state
of interest and its time derivative for a robust finite-time regulation. The
control has only one free design parameter, and the closed-loop system is
shown to be uniformly asymptotically stable in the presence of matched
disturbances. We derive a strict Lyapunov function for the closed control loop
with a bounded exogenous perturbation, and use it for both the control tuning
and analysis of the finite-time convergence. Apart from the numerical results,
a revealing experimental example is also shown in favor of the proposed
control and in comparison with PD and sub-optimal nonlinear damping
regulators.

##  I Introduction

Convergence properties of a transient response and the ability to compensate
for unknown (external) disturbances are often the trade-offs when designing a
feedback control at large. As well known, a PD regulator, which is equivalent
to the state-feedback control for a second-order system, is sufficient for the
unperturbed case provided that an exponential asymptotic convergence is
acceptable. When the system behavior is, however, subject to the static
errors, an integral control action can be added, thus resulting in a standard
PID (or an integral-plus-state feedback) scheme in case of the linear feedback
regulation [1]. From another side, an integral feedback action increases
inherently the system order and can also lead to the so-called wind-up effects
and even destabilization [2]. If the disturbances are not static and possess
certain regularity only (e.g., Lipschitz continuity or/and boundedness), then
the robust feedback regulators like, for instance, high-gain (see, e.g., [3]
and references therein) or finite-time (see [4] and a recent text [5])
controllers can be required. Also the input-to-state stability (ISS) theory
[6] becomes necessary for studying the stability and robustness in presence of
external disturbances, once the overall closed-loop system is nonlinear [7].

Second-order systems capture a relatively large class of the (perturbed)
dynamic processes, especially in mechanical, electro-magnetic, and fluid-flow
domains. In that way, a perturbed double-integrator represents an often
encountered system process. A comparative overview of different feedback
control principles applied to a classical double-integrator can be found in
[8]. Also the sliding-mode control methods, see [9, 10], become relatively
popular for compensating the unknown matched perturbations. Indeed, a feedback
control action proportional to the sign of the control error has already been
recognized in earlier works as particularly efficient for a time-optimal
stabilization [11, 12], and disturbance compensation for perturbed systems
[13]. However, a discontinuous control signal, especially a continuously
switching one, can be undesirable for several controlled processes, especially
with regard to the wear effects, overloading of actuators, noise, and energy
consumption. Therefore, a feedback control law that would allow both, a sign
conditioned robustness to the unknown but bounded perturbations, and such
continuous damping action which would prevent steady-state oscillations, can
be desirable. The idea of an inverse output-dependent damping in the second-
order systems was originally proposed in the so-called sub-optimal nonlinear
damping control [14]. Recently, the sub-optimal nonlinear damping control was
augmented by a linear integral term in [15], with additional analysis, so as
to cope with the matched constant disturbances, while the ISS property was
established for time-varying bounded perturbations.

The present work builds up on the idea proposed in [14] and in the sliding
mode control theory introducing the novel nonlinear state feedback control for
the perturbed second-order systems. The dynamic perturbations are assumed to
be upper bounded, and the controller has only one free design parameter,
similarly to [14].

The note is organized as follows. The proposed nonlinear state feedback
control is introduced in Section II, also demonstrating qualitatively the
phase portrait and shape of the converging state trajectory and control value,
and that in comparison with a linear state feedback control (equivalent to PD)
and sub-optimal nonlinear damping control [14]. The main part of the control
system analysis is delivered in Section III. First, the system is shown to be
uniformly asymptotically stable, then the finite-time convergence is formally
proved. Section IV reports the experimental evaluation of the proposed control
and compares it with an equivalent PD and sub-optimal nonlinear damping
controllers. The results reported in this communication are briefly summarized
in Section V.

##  II Novel nonlinear feedback control

In this section, first, we rigorously describe the studied control problem;
second, the control algorithm is presented; third, we finish with the
illustrations of the behavior for the closed-loop dynamics.

###  II-A Problem statement

We consider a sufficiently general class of perturbed second-order systems:

$\frac{\mathrm{d}^{2}y(t)}{\mathrm{d}t^{2}}=u(t)+d(t),\quad t\geq 0$ (1)

with the control signal $u(t)\in\mathbb{R}$ and the matched unknown
disturbance quantity $d(t)\in\mathbb{R},$ where
$\|d\|_{\infty}=\text{ess}\sup_{t\geq 0}|d(t)|\leq D$ with the upper bound
$D>0$ assumed to be known. Both dynamic system states,
$x(t)=\left(x_{1}(t),x_{2}(t)\right)^{\top}\equiv\bigl{(}y(t),\dot{y}(t)\bigr{)}^{\top}\in\mathbb{R}^{2}$
(including output derivative111Later, we will use the doted symbol, i.e.
$\dot{a}\equiv\mathrm{d}a/\mathrm{d}t$, for denoting the time derivative of a
variable $a$, equally as for the first order differential equations.), are
available for a feedback design, while the robust finite-time stabilization at
$x=0$ is our primary goal.

More precisely, it is required to design a state feedback
$u(t)=u(x_{1}(t),x_{2}(t))$ guaranteeing that there exist
$\sigma,T\in\mathcal{K}$ such that $|x(t)|\leq\sigma(|x(0)|)$ for all $t\geq
0$ and $|x(t)|=0$ for all $t\geq T(|x(0)|)$ for any
$x(0)\in\mathbb{X}=\\{x\in\mathbb{R}^{2}\,:\,x_{1}\neq 0\\}\cup\\{0\\}$ and
$\|d\|_{\infty}\leq D$. So, we assume that at the initial time instant, the
position $x_{1}(0)$ is not zero, and it is required to bring $x(t)$ to zero in
a finite time independently in the presence of a properly bounded disturbance
$d$.

Note that in case of an additional (inertial) scaling factor in the left hand-
side of (1), like in the experimental study of a moving lumped mass shown in
Section IV, this should be additionally incorporated as a normalization factor
of the control gains, cf. [16].

###  II-B Control law

The proposed second-order nonlinear state-feedback control

$u(t)=\begin{cases}-|x_{1}(t)|^{-1}\bigl{(}\gamma
x_{1}(t)+|x_{2}(t)|x_{2}(t)\bigr{)}&x_{1}(t)\neq 0\\\
-\gamma\,\text{sign}\bigl{(}x_{1}(t)\bigr{)}&x_{1}(t)=0\end{cases},$ (2)

where $\gamma>D$ is the single control design parameter, yields for
$x(t)\in\mathbb{X}\setminus\\{0\\}$ the closed-loop system dynamics to

$\displaystyle\frac{\mathrm{d}x_{1}(t)}{\mathrm{d}t}=x_{2}(t),$ (3)
$\displaystyle\frac{\mathrm{d}x_{2}(t)}{\mathrm{d}t}=\qquad\;\;$
$\displaystyle\begin{cases}-|x_{1}(t)|^{-1}\bigl{(}\gamma
x_{1}(t)+|x_{2}(t)|x_{2}(t)\bigr{)}+d(t)&x_{1}(t)\neq 0\\\
-\gamma|x_{1}(t)|^{-1}x_{1}(t)+d(t)&x_{1}(t)=0\end{cases}.$

As we can conclude, the proposed control is discontinuous at $x_{1}=0$, and in
the set $\mathbb{X}$ the irregularity happens at the origin only, which is the
equilibrium of the system. Moreover, there are two discontinuous terms: the
first one is the conventional sign function
$\text{sign}(x_{1}(t))=x_{1}(t)|x_{1}(t)|^{-1}$ widely used in the sliding
mode control theory, the second discontinuity comes from the fraction
$|x_{2}(t)|x_{2}(t)|x_{1}(t)|^{-1}$. For definition of solutions in the face
of the first kind of discontinuity, the Filippov’s theory can be used with
convex embedding [9, 10]. However, the second term becomes unbounded at
$x_{1}=0$ and another approach should be utilized. Below, in order to
substantiate existence and uniqueness of solutions for (3), we will
demonstrate similarity of the system (3) via a time change to another system
with well defined solutions, as we have already done for the sub-optimal
nonlinear damping control in [15]. Next, the stability and performance of the
closed-loop system can be analyzed. Before doing it, we proceed with some
illustrative simulation results showing the intuition behind the proposed
design and qualitative properties of the state trajectories.

###  II-C Illustrations

An exemplary phase-portrait of the closed-loop control system (3) is
illustrated in Fig. 1. Here the control gain factor $\gamma=100$ is assigned.
Due to a symmetric behavior with respect to the origin, the trajectories for
different initial conditions are shown in the 1st and 4th quadrants only.

Figure 1:  Phase portrait of the closed-loop control system (3).

The trajectories of $x_{2}(t)$ and the control inputs $u(t)$ are compared
(qualitatively) in Fig. 2 between the proposed control (3), sub-optimal
nonlinear damping control [14], and critically damped linear state feedback
control.

Figure 2:  Response of the dynamic state $x_{2}(t)$ in (a) and control signal
$u$ in (b) for the proposed control (3), sub-optimal nonlinear damping control
[14], and critically damped linear state feedback control.

All three controllers are assigned with the same output feedback gain factor
$\gamma=100$, while the response of $x_{2}(t)$ is depicted on the logarithmic
scale for the sake of a better comparison in Fig. 2 (a). The corresponding
control signals $u(t)$ are shown in Fig. 2 (b). All three feedback controllers
disclose the same control value at the beginning, this is due to the equal
gaining factor $\gamma$ and initial value ($x_{1}(0)=-1$), then one can
recognize that the control (3) uses the range of the admissible values more
efficiently during the convergence phase.

##  III Analysis

We start this section with proof of existence and uniqueness of solutions in
(3). Next, the uniform asymptotic stability at the origin will be
substantiated. The rate of convergence is evaluated at the end.

###  III-A Existence of solutions

In the Appendix, for the system (8) operating in time $\tau\geq 0$ the
following properties are demonstrated for any $x(0)\in\mathbb{R}^{2}$ and the
inputs with $\|d\|_{\infty}\leq D$:

  1. 1.

global existence and uniqueness of solutions for all $t\geq 0$,

  2. 2.

invariance of the sets $\mathcal{X}=\\{x\in\mathbb{R}^{2}:x_{1}=0\\}$ and
$\mathbb{X}=\left(\mathbb{R}^{2}\setminus\mathcal{X}\right)\cup\\{0\\}$,

  3. 3.

existence of the equilibrium at the origin, which is uniformly globally
asymptotically stable.

Introducing a new (scaled) time argument

$t=\int\limits_{0}^{\tau}\bigl{|}x_{1}(s)\bigr{|}ds,$

which is well-defined on the non-negative real numbers, i.e.,
$t\in\mathbb{R}_{+}$ for $\tau\in\mathbb{R}_{+}$ (indeed, $t=0$ for $\tau=0$,
and
$\lim_{\tau\to+\infty}\int\limits_{0}^{\tau}\bigl{|}x_{1}(s)\bigr{|}ds\leq+\infty$
since $\bigl{|}x_{1}(\tau)\bigr{|}\to 0$ as $\tau\to+\infty$), and that
implies

$\mathrm{d}t=\bigl{|}x_{1}(\tau)\bigr{|}\mathrm{d}\tau$

transforming (8) to the closed-loop system (3). Therefore, the systems (8) and
(3) have the same solutions, where the only difference is that each solution
becomes scaled along the independent time variable based on the norm of
output.

Since (8) has a unique solution asymptotically approaching the origin for any
initial condition $x(0)\in\mathbb{X}$ and input with $\|d\|_{\infty}\leq D$,
(3) also admits this solution while demonstrating the same behavior in finite
or infinite time. In the former case, since the origin is the equilibrium of
(3), the zero solution can be extended for all $t\geq 0$. Therefore, the
following result has been proven:

######  Proposition 1.

The system (3) admits an unique solution for any $x(0)\in\mathbb{X}$ and
$\|d\|_{\infty}\leq D$ well defined with $x(t)\in\mathbb{X}$ for all $t\geq
0$.

In fact, via the similarity of (8) and (3), a stronger property can be proven:
uniform asymptotic stability of the origin for (3) in $\mathbb{X}$. However,
below we would like to apply the Lyapunov function method directly to (3), so
as to establish this property with demonstration of the finite-time
convergence rate.

###  III-B Uniform asymptotic stability

We consider first an unperturbed closed-loop dynamics (3), i.e., with $d=0$.
For that system, first, we suggest the non-strict Lyapunov function candidate

$E(x)=\gamma|x_{1}|+\frac{1}{2}\,x_{2}^{2},$ (4)

which is positive definite for all $x\in\mathbb{R}^{2}$. Evaluating the time
derivative of (4) along the solutions of (3) we obtain

$\dot{E}=\gamma\,\mathrm{sign}(x_{1})\dot{x}_{1}+x_{2}\dot{x}_{2},$

and substituting the right hand-side of (3) into $\dot{E}$ results in

$\dot{E}=-|x_{2}|x_{2}^{2}\leq 0,$

which is negative semi-definite only. It is enough to guarantee boundedness of
the state trajectories $x(t)$ in the closed-loop system for all $t\geq 0$. In
order to ensure the convergence of the unperturbed model (3) to the
equilibrium, we may combine the non-strict Lyapunov function (4) with the
LaSalle invariance principle. For $x\in\mathbb{X}$ with $x_{2}=0$ (the case
when $\dot{E}=0$), the dynamics of trajectories become

$\dot{x}_{1}(t)=0,\quad\dot{x}_{2}(t)=-\gamma\,\text{sign}\left(x_{1}(t)\right)\neq
0,$ (5)

that is always true as long as $x_{1}\neq 0$. Then, the largest invariant set
of (3) with $d=0$ where $\dot{E}=0$ is $\\{0\\}$. This may show that the zero
equilibrium of an unperturbed system (3) is asymptotically stable.

Now, we are in the position to prove the uniform asymptotic stability of the
origin for the perturbed system (3) with $d\neq 0$.

######  Theorem 2.

The origin for the system (3) is uniformly asymptotically stable for any
$x(0)\in\mathbb{X}$ and $\|d\|_{\infty}\leq D$ if

$\gamma>D^{1.5}+D+\frac{1}{2}.$ (6)

###### Proof:

As in the Appendix for the system (8), consider for (3) a strict Lyapunov
function candidate

$V(x)=\frac{1}{2}\,\left(\begin{array}[]{c}z_{1}\\\
x_{2}\end{array}\right)^{\top}\,\left(\begin{array}[]{cc}2\,\gamma&\varepsilon\\\
\varepsilon&1\end{array}\right)\left(\begin{array}[]{c}z_{1}\\\
x_{2}\end{array}\right)$

with $z_{1}=\sqrt{|x_{1}|}\,\mathrm{sign}(x_{1})$, which is positive definite
for $\sqrt{2\gamma}>\varepsilon$. This Lyapunov function is continuously
differentiable for $x\in\mathbb{X}\setminus\\{0\\}$. Taking the time
derivative of $V$ one obtains

$\displaystyle\dot{V}$ $\displaystyle=$ $\displaystyle
x_{2}\,d-x_{2}^{2}\frac{|x_{2}|}{|x_{1}|}-\varepsilon\bigl{(}\gamma-\mathrm{sign}(x_{1})\,d\bigr{)}\,\sqrt{|x_{1}|}$
$\displaystyle+\,\varepsilon\frac{x_{2}^{2}}{\sqrt{|x_{1}|}}\,\Bigl{(}0.5-\mathrm{sign}(x_{1})\mathrm{sign}(x_{2})\Bigr{)}.$

Using the _Young’s inequality_ [17] which postulates

$ab\leq\frac{1}{p}\,a^{p}+\frac{p-1}{p}\,b^{p/(p-1)}$

for any $a,b\in\mathbb{R}_{+}$ and $p>1$, one can show that

$\dot{V}\leq-\biggl{(}\frac{2}{3}-\varepsilon\biggr{)}\,x_{2}^{2}\frac{|x_{2}|}{|x_{1}|}-\biggl{[}\varepsilon\Bigl{(}\gamma-\frac{1}{2}-D\Bigr{)}-\frac{2}{3}D^{1.5}\biggr{]}\,\sqrt{|x_{1}|}$
(7)

for any $x\in\mathbb{X}$ and $|d|\leq D$. It can be seen, that taking

$\frac{2}{3}\frac{D^{1.5}}{\gamma-\frac{1}{2}-D}<\varepsilon<\frac{2}{3}$

both expressions in the brackets on the right hand-side of $\dot{V}$ are
positive under (6), hence, $\dot{V}<0$ that completes the proof. ∎

The low bound of the control gain $\gamma$, which is satisfying the parametric
condition (6), is shown as a function of the disturbance value’s magnitude,
i.e., $D$, in Fig. 3.

Figure 3:  Low bound of control gain depending on perturbation bound.

###  III-C Finite-time convergence

In order to prove the finite-time convergence (see, e.g., [4] for details) of
the system (3), which can be supposed from the numerical results given in Fig.
2 (a), we make further use of the strict Lyapunov function $V$.

######  Theorem 3.

The system (3) is finite-time convergent to the origin for any
$x(0)\in\mathbb{X}$ and $\|d\|_{\infty}\leq D$ provided that (6) is satisfied.

###### Proof:

The estimate (7) implies boundedness of solutions of (3) in $\mathbb{X}$ for
the inputs with magnitude smaller than $D$, with their asymptotic convergence
to the origin. Let us consider an auxiliary variable
$\zeta(t)=x_{2}(t)x_{1}(t)^{-1}$, which is well defined on the solutions of
the system (3) (in the set $\mathbb{X}\setminus\\{0\\}$), and it has the
dynamics:

$\displaystyle\dot{\zeta}=\frac{d\,x_{1}-\bigl{(}\gamma|x_{1}|+\text{sign}(x_{1})|x_{2}|x_{2}\bigr{)}-x_{2}^{2}}{x_{1}^{2}}$
$\displaystyle=\frac{d\,\text{sign}(x_{1})-\gamma}{|x_{1}|}-|\zeta|\zeta-\zeta^{2}$
$\displaystyle=\frac{d\,\text{sign}(x_{1})-\gamma}{|x_{1}|}-\begin{cases}2\zeta^{2}&\zeta>0\\\
0&\zeta\leq 0\end{cases}.$

Note that $(d\,\text{sign}(x_{1})-\gamma)|x_{1}|^{-1}<0$ due to (6), which
implies that for any $x(0)\in\mathbb{X}$ and $\|d\|_{\infty}\leq D$ there is
$t_{1}\geq 0$ such that $\zeta(t)<0$ for all $t\geq t_{1}$, i.e., the
trajectories enter in the second or the fourth quadrants after finite
transients. Next, consider another auxiliary variable
$z(t)=x_{2}^{2}(t)|x_{1}(t)|^{-1}$, which admits the dynamics:

$\displaystyle\dot{z}=\frac{|x_{2}|}{|x_{1}|}\,\Bigl{(}2\text{sign}(x_{2})\bigl{(}d-\gamma\,\text{sign}(x_{1})\bigr{)}$
$\displaystyle-\bigl{(}2+\text{sign}(x_{1})\text{sign}(x_{2})\bigr{)}z\Bigr{)}.$

For $t\geq t_{1}$ and taking into account (6),
$\text{sign}(x_{2})(d-\gamma\text{sign}(x_{1}))>0$ out the origin (the lines
$x_{1}=0$ and $x_{2}=0$ can be reached only there). Hence, for any
$x(0)\in\mathbb{X}$ and $\|d\|_{\infty}\leq D$, there is $t_{2}\geq 0$ such
that $z(t)\geq\theta>0$ for all $t\geq t_{2}$ (while the trajectories do not
converge to zero), where $\theta$ depends only on the choice of $\gamma$. So
there are two possibilities, either the trajectories approach the origin in a
finite time, or the variable $z(t)$ becomes strictly positive and separated
with zero, which returning to (7) implies for $t\geq t_{2}$:

$\displaystyle\dot{V}\leq-\biggl{(}\frac{2}{3}-\varepsilon\biggr{)}\,\theta|x_{2}|-\biggl{[}\varepsilon\Bigl{(}\gamma-\frac{1}{2}-D\Bigr{)}-\frac{2}{3}D^{1.5}\biggr{]}\,\sqrt{|x_{1}|}$
$\displaystyle\leq-\kappa\sqrt{V}$

for some $\kappa>0$. The latter estimate implies a finite-time convergence of
$V$ to zero, according to [4]. ∎

The controlled output behavior of (3) is visualized in Fig. 4 (a). Here a
periodic disturbance $d(t)=50\sin(20\pi t)$ is used with the amplitude equal
to half of the control gain, i.e., $\gamma=100$. Both the disturbance signal
and the control value are shown opposite to each other in Fig. 4 (b).

Figure 4:  Output response of the perturbed control system (3) in (a), and the
control signal versus disturbance value in (b).

##  IV Experimental control results

The proposed control (3) is experimentally evaluated on the electro-magneto-
mechanical actuator system with one translational degree of freedom, see Fig.
5.

Figure 5:  Experimental setup of the actuator system.

Technical details, including the identified system model and most relevant
hardware parameters of the setup, can be found in [16]. The sampling rate is
$10$ kHz. We note that the overall upper-bounded perturbation $d$ is due to
the constant gravity term, but also due to the time-varying friction and
cogging force by-effects of the electro-magneto-mechanical actuator.

The step response of the control (3) (denoted further as NLG owing to the
nonlinear gains) is compared once with a standard proportional-derivative (PD)
control

$u(t)_{\mathrm{PD}}=-\gamma x_{1}-\sigma x_{2},$

and once with the sub-optimal nonlinear damping control [14] (denoted further
as SOND). The PD control is first designed so that the unperturbed closed-loop
control system has two critically damped real poles (denoted further as
$\mathrm{PD}_{1}$), i.e., with an optimal value of $\sigma$ that corresponds
to the unity damping ratio.

TABLE I:  Evaluated controllers param. \ cont. |  $\mathrm{PD}_{1}$ |  $\mathrm{PD}_{2}$ |  $\mathrm{SOND}$ |  $\mathrm{NLG}$
---|---|---|---|---
$x_{1}$-gain $\gamma$ |  $500$ |  $750$ |  $750$ |  $25$
$x_{2}$-gain $\sigma$ |  $2$ |  $4$ |  $-$ |  $-$

Then, the output control gain $\gamma$ is enhanced (denoted further as
$\mathrm{PD}_{2}$), while allowing also for the conjugate complex pole pair
with a transient overshoot. The parameters of all four evaluated controllers
are listed in Table I. For a fair comparison, the SOND control is assigned
with the same gaining factor $\gamma=750$ as the faster PD-control
$\mathrm{PD}_{2}$. The feedback gain of the NLG-control $\gamma=25$ is
intentionally set an order of magnitude lower than the feedback gains of the
PD and SOND controllers. This is for demonstrating that the NLG-control
behaves less aggressive during the transient response and, at the same time,
achieves a much higher accuracy at steady-state. All experimentally measured
control responses, to the same reference step of $y_{\mathrm{ref}}=8$ mm, are
shown opposite each other in Fig. 6.

Figure 6:  Measured step response of PD, SOND, and NLG controllers.

From the response of NLG-control, one can also recognize a nearly uniform
behavior of $x_{2}(t)$ during the whole transient phase, cf. Fig. 2 (a). Note
that the corresponding inclination of the $x_{1}(t)$ slope depends on the
assigned $\gamma$ parameter. It is also worth noting that a residual non-zero,
however minor, steady-state error of the NLG-control is due to the sensor’s
bias and noise, while the latter is additionally impacting the used $x_{2}(t)$
signal, which is obtained via differentiation and filtering of the measured
$x_{1}(t)$.

##  V Summary

A novel robust nonlinear state-feedback control for second-order systems was
introduced. The control law combines the ideas of the state-dependent
nonlinear damping introduced in [14] and the sliding mode control theory,
which equips the overall feedback law with an output-sign-dependent term able
to compensate for bounded perturbations. The resulted control has only one
tuning parameter, the overall output feedback gain, whose selection should be
based on the upper bound of the perturbation magnitude, and this closed
analytic dependence was shown in the paper. We provided a detailed analysis of
uniform asymptotic stability property of the control system and finite-time
convergence based on the proposed strict Lyapunov function. The control signal
appears to be continuous during the transients, apart from the convergence
instant, where the state value reaches the origin (similarly to some high
order sliding mode control algorithms). Remarkable features of the proposed
control is that it: (i) is equivalent (from the structural viewpoint) to a
standard linear PD feedback control and, at the same time, (ii) can guarantee
for full compensation of the bounded matched disturbances in the second-order
systems, while (iii) providing a global finite-time convergence. The proposed
control was already evaluated experimentally on an electro-magneto-mechanical
actuator system with considerable bounded disturbances.

## References

  * [1] K. J. Åström and T. Hägglund, _Advanced PID control_. ISA - Instrumentation, Systems and Automation Society, 2006.
  * [2] P. Hippe, _Windup in control: its effects and their prevention_. Springer, 2006.
  * [3] A. Ilchmann and E. P. Ryan, “High-gain control without identification: a survey,” _GAMM-Mitteilungen_ , vol. 31, no. 1, pp. 115–125, 2008.
  * [4] S. P. Bhat and D. S. Bernstein, “Finite-time stability of continuous autonomous systems,” _SIAM Journal on Control and optimization_ , vol. 38, no. 3, pp. 751–766, 2000.
  * [5] D. Efimov and A. Polyakov, “Finite-time stability tools for control and estimation,” _Foundations and Trends in Systems and Control_ , vol. 9, no. 2-3, pp. 171–364, 2021.
  * [6] E. Sontag, “Input to state stability: Basic concepts and results,” _Nonlinear and optimal control theory: lectures given at the CIME summer school held in Cetraro, Italy June 19–29, 2004_ , pp. 163–220, 2008.
  * [7] H. Khalil, _Nonlinear Systems_ , 3rd ed. Prentice Hall, 2002.
  * [8] V. G. Rao and D. S. Bernstein, “Naive control of the double integrator,” _IEEE Control Systems Magazine_ , vol. 21, no. 5, pp. 86–97, 2001.
  * [9] Y. Shtessel, C. Edwards, L. Fridman, and A. Levant, _Sliding mode control and observation_. Springer, 2014.
  * [10] V. Utkin, A. Poznyak, Y. V. Orlov, and A. Polyakov, _Road map for sliding mode control design_. Springer, 2020.
  * [11] A. Fuller, “Relay control systems optimized for various performance criteria,” _IFAC Proceedings Volumes_ , vol. 1, no. 1, pp. 520–529, 1960\.
  * [12] V. T. Haimo, “Finite time controllers,” _SIAM Journal on Control and Optimization_ , vol. 24, no. 4, pp. 760–770, 1986.
  * [13] Y. Z. Tsypkin, _Relay control systems_. Cambridge University Press, 1984.
  * [14] M. Ruderman, “Optimal nonlinear damping control of second-order systems,” _Journal of the Franklin Institute_ , vol. 358, no. 8, pp. 4292–4302, 2021\.
  * [15] D. Efimov and M. Ruderman, “Stability analysis and integral extension of sub-optimal nonlinear damping control,” _International Journal of Robust and Nonlinear Control_ , vol. in press, p. pp, 2024.
  * [16] M. Ruderman, “Motion control with optimal nonlinear damping: from theory to experiment,” _Control Engineering Practice_ , vol. 127, p. 105310, 2022.
  * [17] W. H. Young, “On classes of summable functions and their Fourier series,” _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , vol. 87, no. 594, pp. 225–229, 1912\.

Consider the second-order system

$\displaystyle\frac{\mathrm{d}x_{1}(\tau)}{\mathrm{d}\tau}$ $\displaystyle=$
$\displaystyle\bigl{|}x_{1}(\tau)\bigr{|}x_{2}(\tau),$ (8)
$\displaystyle\frac{\mathrm{d}x_{2}(\tau)}{\mathrm{d}\tau}$ $\displaystyle=$
$\displaystyle-\gamma
x_{1}(\tau)-\bigl{|}x_{2}(\tau)\bigr{|}x_{2}(\tau)+\bigl{|}x_{1}(\tau)\bigr{|}d(\tau),$

in the time $\tau\geq 0$, with the state
$x(\tau)=(x_{1}(\tau)\;x_{2}(\tau))^{\top}\in\mathbb{R}^{2}$ and a bounded
input $d(\tau)\in\mathbb{R}$ with $\|d\|_{\infty}\leq D$ for a given $D>0$,
and with a gain $\gamma>D$. The right hand-side of (8) is locally Lipschitz
continuous in $x(\tau)$ and $d(\tau)$, hence, the solutions are uniquely
defined for any initial condition $x(0)\in\mathbb{R}^{2}$ and a measurable
essentially bounded input $d$ at least locally in time. Moreover, it is easy
to check that the origin is the equilibrium of the system and the submanifold
$\mathcal{X}=\\{x\in\mathbb{R}^{2}:x_{1}=0\\}$ is invariant uniformly in $d$.

######  Lemma 4.

The origin for (8) is uniformly globally asymptotically stable if

$\gamma>D^{1.5}+D+\frac{1}{2}.$ (9)

###### Proof:

Let us consider the following Lyapunov function candidate

$\displaystyle V(x)$
$\displaystyle=\frac{1}{2}\xi^{\top}(x)\left[\begin{array}[]{cc}2\gamma&\varepsilon\\\
\varepsilon&1\end{array}\right]\xi(x),$ $\displaystyle\xi(x)$
$\displaystyle=\left[\begin{array}[]{c}\sqrt{|x_{1}|}\text{sign}(x_{1})\\\
x_{2}\end{array}\right],$

which is positive definite providing that the parameter
$\varepsilon\in(0,\sqrt{2\gamma})$. The function $V$ is continuously
differentiable for all $x\in\mathbb{R}^{2}\setminus\mathcal{X}$. Note that
$V(x)=\frac{x_{2}^{2}}{2}$ for $x\in\mathcal{X}$ (it is an invariant subspace
for the system, and the trajectories of (8) stay on it while being defined)
and it is straightforward to check that (we denote $V(\tau)=V(x(\tau))$)

$\frac{\mathrm{d}V(\tau)}{\mathrm{d}\tau}=-\bigl{|}x_{2}(\tau)\bigr{|}^{3}<0,$

implying convergence of $x_{2}(\tau)$ to zero. For $x(0)\notin\mathcal{X}$
(note that $\left(\mathbb{R}^{2}\setminus\mathcal{X}\right)\cup\\{0\\}$ is
also invariant for (8)), we obtain

$\displaystyle\frac{\mathrm{d}V(\tau)}{\mathrm{d}\tau}=-(\gamma-\text{sign}(x_{1}(\tau))d(\tau))\varepsilon\sqrt{|x_{1}(\tau)|}|x_{1}(\tau)|-\bigl{|}x_{2}(\tau)\bigr{|}^{3}$
$\displaystyle+\varepsilon\sqrt{|x_{1}(\tau)|}x_{2}^{2}(\tau)\left(\frac{1}{2}-\text{sign}(x_{1}(\tau))\text{sign}(x_{2}(\tau))\right)$
$\displaystyle+\bigl{|}x_{1}(\tau)\bigr{|}x_{2}(\tau)d(\tau)$
$\displaystyle\leq-(\gamma-D)\varepsilon\sqrt{|x_{1}(\tau)|}|x_{1}(\tau)|-\bigl{|}x_{2}(\tau)\bigr{|}^{3}$
$\displaystyle+\varepsilon\sqrt{|x_{1}(\tau)|}x_{2}^{2}(\tau)\left(\frac{1}{2}-\text{sign}(x_{1}(\tau))\text{sign}(x_{2}(\tau))\right)$
$\displaystyle+D|x_{1}(\tau)||x_{2}(\tau)|$
$\displaystyle\leq-[(\gamma-\frac{1}{2}-D)\varepsilon-\frac{2}{3}D^{1.5}]\sqrt{|x_{1}(\tau)|}|x_{1}(\tau)|$
$\displaystyle-(\frac{2}{3}-\varepsilon)\bigl{|}x_{2}(\tau)\bigr{|}^{3},$

where on the last step the following inequalities have been utilized (that are
obtained from Young’ s inequality [17]):

$\displaystyle\sqrt{|x_{1}(\tau)|}x_{2}^{2}(\tau)\leq\frac{1}{3}\sqrt{|x_{1}(\tau)|}|x_{1}(\tau)|+\frac{2}{3}\bigl{|}x_{2}(\tau)\bigr{|}^{3},$
$\displaystyle
D|x_{1}(\tau)||x_{2}(\tau)|\leq\frac{2}{3}D^{1.5}\sqrt{|x_{1}(\tau)|}|x_{1}(\tau)|+\frac{1}{3}|x_{2}(\tau)|^{3}.$

Hence, for $\gamma>D+\frac{1}{2}$ and

$\frac{2}{3}\frac{D^{1.5}}{\gamma-\frac{1}{2}-D}<\varepsilon<\min\\{\frac{2}{3},\sqrt{2\gamma}\\}$

the desired property

$\frac{\mathrm{d}V(\tau)}{\mathrm{d}\tau}<0$

while $V(\tau)\neq 0$ is obtained. The interval for the values of
$\varepsilon$ is not empty provided that the condition (9) is verified.
Indeed, this condition implies that $\gamma>0.5$ for any $D>0$, hence,
$\min\\{\frac{2}{3},\sqrt{2\gamma}\\}=\frac{2}{3}$ and
$\frac{2}{3}\frac{D^{1.5}}{\gamma-\frac{1}{2}-D}<\frac{2}{3}$ is equivalent to
(9). Combining the derived estimates for
$\frac{\mathrm{d}V(\tau)}{\mathrm{d}\tau}$ on the submanifold of $\mathcal{X}$
and outside, the uniform global asymptotic stability property of the origin
for (8) follows the conventional results [7]. ∎