# Rainfall-runoff prediction using a Gustafson-Kessel clustering based Takagi-
Sugeno Fuzzy model

Subhrasankha Dey*1, Tanmoy Dam2 *Corresponding author1Subhrasankha Dey, Dept.
of Infrastructure Engineering, The University of Melbourne, Parkville,
Australia<EMAIL_ADDRESS>Tanmoy Dam, School of Engineering and
Information Technology, University of New South Wales Canberra, Australia
<EMAIL_ADDRESS>
###### Abstract

A rainfall-runoff model predicts surface runoff either using a physically-
based approach or using a systems-based approach. Takagi-Sugeno (TS) Fuzzy
models are systems-based approaches and a popular modeling choice for
hydrologists in recent decades due to several advantages and improved accuracy
in prediction over other existing models. In this paper, we propose a new
rainfall-runoff model developed using Gustafson-Kessel (GK) clustering-based
TS Fuzzy model. We present comparative performance measures of GK algorithms
with two other clustering algorithms: (i) Fuzzy C-Means (FCM), and (ii)
Subtractive Clustering (SC). Our proposed TS Fuzzy model predicts surface
runoff using: (i) observed rainfall in a drainage basin and (ii) previously
observed precipitation flow in the basin outlet. The proposed model is
validated using the rainfall-runoff data collected from the sensors installed
on the campus of the Indian Institute of Technology, Kharagpur. The optimal
number of rules of the proposed model is obtained by different validation
indices. A comparative study of four performance criteria: Root Mean Square
Error (RMSE), Coefficient of Efficiency (CE), Volumetric Error (VE), and
Correlation Coefficient of Determination (R) have been quantitatively
demonstrated for each clustering algorithm.

###### Index Terms:

Rainfall-runoff prediction, Takagi Sugeno Fuzzy inference, Gustafson-Kessel
clustering, Comparative error measure

##  I Introduction

A rainfall-runoff is an overland water flow due to precipitation happening
over a drainage basin. A rainfall-runoff model predicts accumulated overland
flow at an outlet of a drainage basin due to excess rainfall. A rainfall-
runoff model can be developed using either a physically-based model or a data-
driven black-box model. A physically-based model is constructed by the mass,
momentum, and energy transformation equations such that the parameters of the
model are directly related to the characteristics of a drainage basin (or a
catchment area) [1]. However, such models ignore frequent topographical
changes of a catchment due to urbanization/human interventions and demand
additional data (e.g. initial soil moisture, land use, evaporation and
infiltration data, distribution, and rainfall duration) to improve model
accuracy [2]. In contrast to physically-based approaches, data-driven models
can handle the topographic variation of a catchment and can produce output
even with the absence of additional parameters and data, those are needed by a
physically-based model [3].

Recently, Fuzzy rule-based models are becoming a popular data-driven approach
for rainfall-runoff modelling [2, 4, 5, 6, 3]. One of the important Fuzzy
rule-based rainfall-runoff models is the Neuro-Fuzzy system, known as Adaptive
Network-based Fuzzy Inference System (ANFIS) [7]. ANFIS has been shown
superiority when compared with other data-driven models such as Artificial
Neural Network (ANN), Auto-Regressive Moving Average (ARMA), and Auto-
Regressive with exogenous inputs (ARX) models. However, ANFIS is reported to
be computationally expensive and tends to be over-fitting [8, 3].

Performance of the existing Fuzzy rule-based models varies with the
topographic nature of a catchment [3]. Further, the models often suffer from
uncertainty in the selection of lag-time, variability of the training data in
the presence of many influencing parameters of a hydrological process. Hence,
a comparative performance measure is adapted by researchers in order to
achieve higher accuracy to predict runoff due to rainfall [2, 9]. Also,
Takagi–Sugeno (TS) fuzzy systems remains a popular choice of capturing the
non-linear relationships between rainfall and runoff in a catchment [2, 9].

The main contributions of this paper are as follows,

  * •

Gustafson-Kessel’s (GK) clustering method has been adapted to learn the fuzzy
structure, mainly premise and consequence parameters in TS fuzzy model
identification. To obtain the optimal number of rules directly from the data,
we used different cluster validating indices.

  * •

To validate the GK-based TS fuzzy model using the data obtained at a catchment
area of Indian Institute of Technology, Kharagpur (an academic campus) [10].
We evaluate our proposed model’s performance in the multi-step-ahead runoff
prediction scheme. We have also compared the performance with other existing
clustering-based Fuzzy methods.

##  II Literature review

In earlier research, linear regression-based runoff prediction models are
found [11]. However, nonlinear dynamics of runoff due to rainfall were not
considered and hence the model performance is not adequate in linear
regression. In recent time, the following nonlinear data-driven approaches are
becoming popular (i) Artificial Neural Network (ANN) [12, 13, 11] and (ii)
Fuzzy rule-based Neuro-Fuzzy Model [4, 5, 6, 14].

The Fuzzy rule-based models are popularly based on the trial and error method
to find out the optimal number of rules [11, 15]. In a Fuzzy rule-based
approach, membership functions and Fuzzy rules are determined with some degree
of uncertainty which is present in a hydrological process [3]. An improvement
of a Fuzzy rule-based system is found in Neuro-Fuzzy models where a neural
network is used to train the model’s membership functions and Fuzzy rules [16,
3]. Neuro-Fuzzy systems are often based on Takagi–Sugeno-Kang structure [3].
ANFIS is a widely popular Neuro-Fuzzy system and is used in a number of data-
driven rainfall-runoff modelling despite the overfitting issues of finding
suitable membership functions [16, 17, 18, 19]. Other studies which have
employed Neuro-Fuzzy Models (NFMs) for flow forecasting include DENFIS
(Dynamic Evolving Neural-Fuzzy Inference System) models [20]. Wavelet-based
approaches are also integrated with the ANFIS model to handle the uncertainty
in the observed data for removing the noise samples [21]. The effect of lag
time on rainfall-runoff modeling by ANFIS has been studied earlier [22].
Additionally, an input selection method based on correlation and mutual
information analysis is able to identify an optimum set of rainfall inputs for
rainfall-runoff modeling by ANFIS [22]. Sometimes, non-sequential rainfall
antecedents can produce better results compared to sequential rainfall inputs
[22]. Hence, selection of the type and number of inputs are also important in
a TS Fuzzy rule-based rainfall-runoff model [23]. Rainfalls are considered by
some researchers as the only inputs of the data-driven model [24, 25]. When
the flow measurement is unavailable for sub-catchments in a catchment,
rainfall data can be chosen as inputs to identify the more contributing sub-
catchments [8]. However, the performance of ANFIS could be affected due to
unnecessary complexity when so many inputs are involved [22]. The
computational time can be significantly reduced using a fewer number of
inputs.

A combination of rainfall and runoff are also used as inputs to the Fuzzy
models [15, 26] with increased computation complexity and overfitting errors.
Some studies have suggested pruning of the unnecessary inputs by selecting a
narrower time window around the most correlated rainfall antecedent with
runoff [27, 15]. Hence, in this paper, we develop a simple four-input-based
Takagi-Sugeno Fuzzy model to reduce computational complexity and reducing
overfitting errors. We present a comparative study of three different Fuzzy
clustering algorithms based the proposed TS Fuzzy model using real data
collected at a pilot study area of a drainage basin.

##  III General Structure of a TS Fuzzy model

In TS Fuzzy model, any nonlinear function is approximated by a set of linear
models where each linear model is described as a rule [28]. The structure of
each rule is represented by premise and consequent variables. An input
variable transformed into Fuzzy Membership Function (MF) is known as premise
variable, whereas the output of TS Fuzzy model is known as consequent
variable. The consequent output is described by a linear functional relation
between premise variables. In a TS Fuzzy model , each rule is described as
follows,

$\displaystyle{P^{i}}:\,IF\,$
$\displaystyle{x_{k1}}\,is\,{A_{1}}^{i}\,{\rm{AND}}\,IF\,{x_{k2}}\,is\,{A_{2}}^{i}...\,{\rm{AND}}\,IF\,{x_{kn}}\,is\,{A_{n}}^{i}$
$\displaystyle THEN\,{y_{k}}^{i}$
$\displaystyle={a_{0}}^{i}+{x_{k}}_{1}.{a_{1}}^{i}+...+{x_{kn}}.{a_{n}}^{i}$
$\displaystyle={a_{0}}^{i}+\sum\limits_{j=1}^{n}{{x_{kj}}.{a_{j}}^{i}},\,\,j=1,2,..,n,\,i=1,2,..,C$
$\displaystyle=[1\,\,{x_{k}}^{T}]{\theta^{i}}\,=\bar{x}{\theta^{i}}$ (1)

where $P^{i}$ is the structure of the rule $(i=1,2,...,Crule)$ and $n$ is the
number of input variable. ${{\bf x}_{k}}=[{x_{k1}},{x_{k2}},...{x_{kn}}]$ be
the input variables or premise variables of the model and $y_{k}$ is the
consequent’s output of $i^{th}$ rule.
${\theta^{i}}=[{p_{0}}^{i},{p_{1}}^{i},...,{p_{n}}^{i}]\in{R^{(n+1)}}$ is the
coefficients of consequent variable. The dumbbell shape MF (i.e. Gaussian MF)
for each premise variable is considered here as,

$\displaystyle{A_{j}}^{i}({x_{kj}})=\exp[-\frac{{{{({x_{kj}}-{\vartheta_{j}}^{i})}^{2}}}}{{{{({\sigma_{j}}^{i})}^{2}}}}]\in[{\rm{0}},{\rm{1}}],$
$\displaystyle\textrm{where, }(i=1,2,..,C;j=1,2,..,n)$ (2)

where
$V=[{\vartheta_{1}}^{i},{\vartheta_{2}}^{i},...,{\vartheta_{n}}^{i}]\in{R^{n}}$
is the vector of mean of the GMF. The width is represented by
$[{\sigma_{1}}^{i},{\sigma_{2}}^{i},...,{\sigma_{n}}^{i}]\in{R^{n}}$.

The final output of the TS Fuzzy model is represented by weighted average of
each individual firing rule. Let, the system input variables and output
variables are,
${x_{k}}=[{x_{k}}_{1},{x_{k}}_{2},...,{x_{k}}_{n}]{\,{}^{T}}(k=1,2,...,N^{th}$
observation$)$ and ${y_{k}}(k=1,...,N)$ respectively. Therefore, the output of
TS Fuzzy model is given by,

${{{\hat{y}}_{k}}=\frac{{\sum\limits_{i=1}^{C}{{w_{k}}^{i}.{y_{k}}^{i}}}}{{\sum\limits_{i=1}^{C}{{w_{k}}^{i}}}}}$
(3)

where, the membership value (${A_{j}}^{i}({x_{kj}})$)of $j^{th}$ premise
variable of $k^{th}$ sample is defined in (2). The overall truth value
(${{w_{k}}^{i}}$) of each rule can be calculated by minimum operator or logic
AND operation. The truth value of $k^{th}$ sample is expressed by,

$\displaystyle{w_{k}}^{i}=\mathop{\min}\limits_{j}\\{{A_{j}}^{i}({x_{kj}})\\},$
$\displaystyle{(i=1,2,...,{\mkern 1.0mu}C{\mkern 1.0mu};{\mkern 1.0mu}{\mkern
1.0mu}j=1,2,...,{\mkern 1.0mu}n\,;\,k=1,2,...,N)}$ (4)

##  IV Fuzzy clustering based TS Fuzzy model

###  IV-A Fuzzy partition matrix from the model data

TS Fuzzy model can be designed from the input-output data by Fuzzy partition
matrix. Let, the data vector is represented as,
$D=\\{{Z_{k}}:{Z_{k}}={[{x_{k}}^{T},{y_{k}}]^{T}},{x_{k}}={[{x_{k1}},{x_{k2}},...,{x_{kn}}]^{T}},k=1,2,...,N\\}$.
Hence, the Fuzzy partition matrix in a data set for premise variables is
represented by clusters with certain degree of membership values. Each
membership value for each cluster is a part of the Fuzzy partition matrix.

Definition Let, $\Im=\\{{Z_{1}},{Z_{2}},...,{Z_{N}}\\}$ be a finite set and
$2\leq C<N$ be number of clusters, then Fuzzy partition matrix for set $\Im$,
is defined as

${U_{fc}}=\\{U\in{R^{C\times
N}}|{\mu_{ik}}\in[0,1];\,\,\sum\limits_{i=1}^{C}{{\mu_{ik}}=1;}\,\,0\prec\sum\limits_{k=1}^{N}{{\mu_{ik}}<N,\,\,\forall}i\\}$

where, ${\mu_{ik}}(i=1,2,...,C\,and\,k=1,2,...,N)$ be Fuzzy membership value
of ${Z_{k}}$ sample is belonging to $\emph{i}^{th}$ cluster.

###  IV-B GK clustering algorithm based Fuzzy partition matrix

The whole input-output space is divided by Fuzzy clustering algorithm to
achieve the Fuzzy partition matrix. Therefore, many clustering algorithms
based TS Fuzzy models are found in [29], [30], [31]. FCM clustering algorithm
is popular of its kind [32]. The modified version of FCM algorithm is known as
GK (Gustafson-Kessel) Algorithm [33]. Cluster covariance values of each class
is updated by the norm inducing distance metric. The geometrical shape
provided by GK clustering algorithm is elliptical shape in nature whereas the
FCM gives hyper-spherical shape. The objective function of GK algorithm is
given by,

${\begin{array}[]{l}{J_{GK}}(D;U,{\rm{\theta}},{\rm{v,}}\,{M_{i}})=\sum\limits_{i=1}^{C}{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}}}{G_{ik}}^{2}\\\
subject\,to,\,\,\sum\limits_{i}^{C}{{\mu_{ik}}}=1\end{array}}$ (5)

where the distance metric is
${G_{ik}}^{2}={[{Z_{k}}-{\textbf{v}^{i}}]^{T}}{M_{i}}[{Z_{k}}-{\textbf{v}^{i}}]$,
${{\bf{v}}^{i}}\in{R^{(n+1)}}$ is cluster prototype in the input-output
product space and ${M_{i}}$ is symmetric positive definite matrix, represents
as volume size of $i^{th}$ cluster. Equation (5) is a constrained optimization
problem. Applying the Lagrangian multipliers ($\lambda$) to convert it into an
unconstrained optimization form,

${\begin{array}[]{l}L(U,\lambda,{\rm{\theta}},{\rm{v}})={J_{GK}}(D;U,{\rm{\theta}},{\rm{v,}}\,{M_{i}})-\sum\limits_{i=1}^{C}{{\beta_{i}}}(|{M_{i}}|-\rho)\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\sum\limits_{k=1}^{N}{{\lambda_{k}}(\sum\limits_{i=1}^{C}{{\mu_{ik}}-1})}\end{array}}$
(6)

For minimization of (6) with the respect to ${\mu_{ik}}$, assuming
${G_{ik}}^{2}\succ 0$, the final expression is as follows,

${{\mu_{ik}}=\frac{1}{{\sum\limits_{q=1}^{C}{{{\left({\frac{{{G_{ik}}^{2}}}{{{G_{qk}}^{2}}}}\right)}^{\frac{1}{{m-1}}}}}}}}$
(7)

Similarly, the cluster center ($v^{i}$) is obtained by taking derivative of
(5) w.r.t. ${\textbf{v}^{i}}$ as,

${\begin{array}[]{l}\frac{{\partial
L(...)}}{{\partial{{\rm{v}}^{i}}}}=\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}[{x_{k}}-{{\rm{v}}^{i}}]}=0\\\
\Rightarrow\,\,\,{{\rm{v}}^{i}}=\frac{{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}{Z_{k}}}}}{{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}}}}\end{array}}$
(8)

Taking partial derivative (6) w.r.t $M^{i}$,

${\begin{array}[]{l}\frac{{\partial
L(...)}}{{\partial{M_{i}}}}=\sum\limits_{k=1}^{N}{{\mu_{i}}{{{}_{k}}^{m}}[{{\rm{Z}}_{k}}-{{\rm{v}}^{i}}]}{[{{\rm{Z}}_{k}}-{{\rm{v}}^{i}}]^{T}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-{\beta_{i}}\left|{{M_{i}}}\right|{M_{i}}^{-1}=0\\\
\Rightarrow{M_{i}}^{-1}=\frac{{\sum\limits_{k=1}^{N}{{\mu_{i}}{{{}_{k}}^{m}}[{{\rm{Z}}_{k}}-{{\rm{v}}^{i}}]}{{[{{\rm{Z}}_{k}}-{{\rm{v}}^{i}}]}^{T}}}}{{{\beta_{i}}\left|{{M_{i}}}\right|}}\\\
\,\,\,\,\,\,\,\,=\sqrt[n]{{(\frac{1}{{\rho\left|{{F_{i}}}\right|}})}}\,\,\,{F_{i}}.\end{array}}$
(9)

where $F_{i}$ is the covariance matrix. $\rho$ is a constant.

Inputs:

Fuzziness parameter ($m=2$), Iteration number ( $l=0$).

Termination constant ($\xi=0.001$)

${G_{ik}}^{2}$: Calculate the distance metric,
${G_{ik}}^{2}={[{Z_{k}}-{{\bf{v}}^{i}}]^{T}}({\rho}(det{({F_{i}})^{{\raise
2.1097pt\hbox{$1$}\\!\mathord{\left/{\vphantom{1n}}\right.\kern-1.2pt}\\!\lower
2.1097pt\hbox{$n$}}}}){F_{i}}^{-1})[{Z_{k}}-{{\bf{v}}^{i}}]$ .

${U^{0}}\in{\mathbb{R}^{C\times N}}$: Fuzzy partition matrix for first

iteration.

Output:

${U^{l}}$: Fuzzy partition matrix for $l^{th}$ iteration.

Steps:

1. Calculate $v^{i}$ by using (8).

2. Calculate covariance matrix $F_{i}=\frac{{\sum\limits_{k=1}^{N}{\mu_{ik}^{m}\left({Z_{k}-v^{i}}\right)\left({Z_{k}-v^{i}}\right)^{T}}}}{{\sum\limits_{k=1}^{N}{\mu_{ik}^{m}}}}$.

3. Update ${U^{l}}$ with distance metric ${G_{ik}}^{2}$.

${{U^{l}}=\left\\{{\begin{array}[]{*{20}{c}}{\frac{1}{{\sum\limits_{q=1}^{C}{{{\left({\frac{{{G_{ik}}^{2}({\theta^{i}})}}{{{G_{qk}}^{2}({\theta^{q}})}}}\right)}^{\frac{1}{{m-1}}}}}}}}&{{\rm{If}}{\mkern
1.0mu}{G_{ik}}({\theta^{i}})>{\mkern 1.0mu}0}\\\
0&{{\rm{otherwise}}}\end{array}}\right.}$ (10)

until ${\left\|{{U^{l}}-{U^{l-1}}}\right\|_{2}}\leq\xi$ then stop;

otherwise, $l=l+1$ and go to step 1.

Algorithm 1 GK algorithm for a TS Fuzzy model

###  IV-C Cluster validity index

Cluster validity index provides a clear idea about the optimal number of
partitions in the data space. The partition in the data space is obtained by
clustering algorithms viz. FCM, GK and SC [34]. The optimal number of
partitions (subspace) in a data space are obtained by varying the clustering
numbers. However, varying the clustering numbers only may not be sufficient
enough to provide the optimal numbers due to its dependency in the cluster
shape. Therefore, a suitable validity index with proper partitioning algorithm
is required to obtain the optimal numbers. Hence, six different validity
measures are applied to find out the optimal number of clusters in the data
space as follows:

####  IV-C1 Partition Coefficient (PC)

The Partition Coefficient (PC) [35]index measures of overlapping between
clusters. It is defined by,

${{{V_{PC}}=\frac{{\sum\limits_{i=1}^{C}{\sum\limits_{k=1}^{N}{{{\left({{\mu_{ik}}}\right)}^{2}}}}}}{N}}}$
(11)

where ${\mu_{ik}}$ is the membership value of $k^{th}$ data point belonging to
$i^{th}$ cluster. The PC method does not hold any connection between shape of
data. It is concerned only with the partition matrix ($U$) i.e. it does not
have any relation with data space. Therefore, this validity index may not
favorable for highly complex dataset. The index provides an optimal number of
clusters at a point that gives maximum value by varying the cluster number
from $C=2$ to $C_{max}$.

####  IV-C2 Partition Entropy (PE)

The Partition Entropy (PE) index [35] defines the fuzziness between partitions
and is given by as,

${V_{PE}}=-\frac{{\sum\limits_{k=1}^{N}{\sum\limits_{i=1}^{C}{{\mu_{ik}}\log({\mu_{ik}})}}}}{N}$
(12)

where the minimization of ($V_{PE}$) provides the optimal number by varying
cluster number from $C=2$ to $C_{max}$. Similarly, for the $V_{PE}$ index,
that may not provide the exact number of partitions in a dataset because it
does not hold any relation between the partitions shape.

####  IV-C3 Modified Partition Coefficient (MPC)

Both $V_{PC}$ and $V_{PE}$ indices monotonic decreasing tendency with varying
$C$. The modified $PC$ (MPC) [36] can remove the monotonic decreasing tendency
and is given by as,

${{{\rm{V}}_{{\rm{MPC}}}}{\rm{=1-}}\frac{C}{{C-1}}{\rm{(1-}}{{\rm{V}}_{{\rm{PC}}}}{\rm{)}}}$
(13)

The optimal number of clusters is obtained by maximizing $V_{MPC}$ over the
cluster number, varying from $C=2$ to $C_{max}$.

####  IV-C4 Partition Index (SC)

This index [37] is obtained by ratio of the compactness and separation of a
cluster. The summation of the individual cluster validity measure is
normalized by dividing cardinality of each cluster. The index is given by,

${{V_{SC}}=\sum\limits_{i=1}^{C}{\frac{{\sum\limits_{k=1}^{N}{{{\left({{\mu_{ik}}}\right)}^{2}}{{\left\|{{Z_{k}}-{v^{i}}}\right\|}^{2}}}}}{{{N_{i}}\sum\limits_{j=1}^{C}{{{\left\|{{v^{j}}-{v^{i}}}\right\|}^{2}}}}}}}$
(14)

A lower value of the index indicates the optimal numbers.

####  IV-C5 Separation Index (S)

This index [37] is generally do the opposite effect of SC index. Minimum
separation provides the optimal numbers of clusters. It is given by as,

${{V_{S}}=\frac{{\sum\limits_{i=1}^{C}{\sum\limits_{k=1}^{N}{{{\left({{\mu_{ik}}}\right)}^{2}}{{\left\|{{Z_{k}}-{v^{i}}}\right\|}^{2}}}}}}{{{N_{i}}\mathop{\min}\limits_{i,j}{{\left\|{{v^{j}}-{v^{i}}}\right\|}^{2}}}}}$
(15)

####  IV-C6 Xie-Beni (XB)

The Xie-Beni index [34] is based on the compactness of the cluster and minimum
separation variation between the clusters. The optimal value is obtained from
the minimum values of the two ratios by varying the cluster number from $C=2$
to $C_{max}$. The mathematical expression for $V_{XB}$ is given by,

${{V_{XB}}=\frac{{\sum\limits_{i=1}^{C}{\sum\limits_{k=1}^{N}{{{\left({{\mu_{ik}}}\right)}^{2}}{{\left\|{{Z_{k}}-{v^{i}}}\right\|}^{2}}}}}}{{N\mathop{\min}\limits_{i,j}{{\left\|{{v^{j}}-{v^{i}}}\right\|}^{2}}}}}$
(16)

###  IV-D Parameter estimation

In this section, optimal number of rule-base Fuzzy partition matrix obtained
by GK algorithm is used to estimate the premise parameters. Here, the Gaussian
Membership Function (GMF) for premise variables is considered. The mean and
the width (standard deviation) [38, 31, 32, 39, 40] are obtained as follows,

${{\vartheta_{\rm{j}}}^{i}=\frac{{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}{{\rm{x}}_{kj}}}}}{{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}}}}}$
(17)
${{\sigma_{j}}^{i}=\sqrt{\frac{{2*[\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}{{({x_{kj}}-{\vartheta_{j}}^{i})}^{2}}]}}}{{\sum\limits_{k=1}^{N}{{\mu_{ik}}^{m}}}}}\,\,}$
(18)

The premise model parameters are identified from the Fuzzy partition matrix
$(U)$ and model data $(Z)$. The coefficients of consequent part are
identified. The obtained premise parameters are applied to make a global
matrix form of consequent coefficients. The matrix form of coefficients is
expressed as,

${Y=\pi\zeta+\varepsilon}$ (19)

where $\pi$ is the regression vector that is derived from the input vector and
truth value $(w)$. $Y=[{y_{1}},{y_{2}},...,{y_{N}}]$ is the observed data
vector and $\varepsilon$ is the modeling error for $N^{th}$ number of
observations. The global form of the consequent coefficients is
$\zeta=[{p_{0}}^{1},{p_{1}}^{1},...,{p_{n}}^{1},{p_{0}}^{2},{p_{1}}^{2},...,{p_{n}}^{2},...,{p_{0}}^{C},{p_{1}}^{C},...,{p_{n}}^{C}]$.
The regression vector is represented as,

${\pi({x_{k}})=[{\nabla_{1k}}\,\,{\nabla_{1k}}{x_{k1}}\,\,{\nabla_{1k}}{x_{k2}}...\,{\nabla_{1k}}{x_{kn}}\,\,\,...{\nabla_{Ck}}\,...\,{\nabla_{{}_{Ck}}}{x_{kn}}]}$
(20)

where, the truth value is obtained as follows,

${{\nabla_{ik}}\,=\frac{{{w_{k}}^{i}}}{{\sum\limits_{i=1}^{C}{{w_{k}}^{i}}}}}$
(21)

The matrix inversion of Equation (19) may not feasible for all the time to
determine the coefficients of $\zeta$ [31]. Therefore, the Orthogonal Least
Square (OLS) has been applied to the Equation (19) to identify the
coefficients of consequent parameters [32].

##  V Proposed TS Fuzzy model

The selection of input and output variables is the first step to develop a TS-
Fuzzy rainfall-runoff model. Runoff at the outlet of a drainage basin is
directly related to previous values of runoff and rainfall at different
locations. However, the accumulated water due to rainfall takes a variable
time to reach the exit of the basin depending on the catchment topography. The
time interval from the center of mass of rainfall excess to the peak of the
resulting hydrograph is known as lag-time. The effect of lag-time has been
removed from the output data-set by adding the same delay in all the input
variables. We define our proposed model after (removing the lag-time in the
training dataset) as $M_{4}$ model (with four input variable), having the
following data structure:

${{M_{4}}:\,\,\,\,{y_{k}}=f({y_{k-1}},x_{k}^{1},x_{k}^{2},x_{k}^{3})}$ (22)

where $y_{k}$ is the observed runoff and $x_{k}^{1}$, $x_{k}^{2}$ and
$x_{k}^{3}$ are the rainfall collected at different locations at the $k^{th}$
time interval of a time horizon ${0,T}$ hours. Let $\delta^{0}_{pred}$ = time
difference between $k^{th}$ and ${k-1}^{th}$ data point. We define the set
$\delta_{pred}={\delta^{0}_{pred},...,\delta^{i_{N}}_{pred},...}$ such that
$\delta^{i_{N}-1}_{pred}=i_{N}\cdot\delta^{0}_{pred}$ and so on, where $i_{N}$
is an integer. A subset of $\delta_{pred}$ is called as a time ahead runoff
prediction scheme.

###  V-A Flow chart of a GK algorithm based TS Fuzzy model

Figure 1 describes the general architecture of the proposed GK algorithm based
TS Fuzzy model. The model architecture is broadly classified into the
following:

##### Finding optimal rule-base: The GK clustering algorithm (Algorithm1)
along with validation index has been applied to input-output data set to
create optimal Fuzzy subspaces (Fuzzy rules). Each optimal Fuzzy subspace is
associated to a Fuzzy rule i.e. the output of optimal Fuzzy ubspaces is
described by Fuzzy partition matrix (${U}$). Once the optimal subspace is
obtained from the both algorithms, then obtained $U$ has been used to estimate
the model parameters.

##### Estimating TS Fuzzy model parameters: Gaussian membership function
parameters (mean and width) are obtained by using eq.(17) and eq.(18)
respectively. We have formulated a global matrix form of regression vector to
identify the consequent parameters. OLS is applied to eq.(19) estimate the
consequent parameters.

Figure 1: Architecture of the GK algorithm based TS Fuzzy model

##  VI Data and validation technique

A part of the campus of IIT Kharagpur is used as a pilot test-bed for the data
acquisition process as shown in Figure 2. The validation data are obtained
from an real-time web-based hydro-meteorological sensor network platform that
exists in the study area [10]. In Figure 2, the red markers are indicating
locations of tipping bucket rain-gauges for on-line measurement of rainfall,
and the blue up arrow indicates outlet of the drainage basin, where a
rectangular weir with pressure type digital water depth sensor is installed to
calculate pressure head (or water level). Pressure head is converted to
surface runoff flow using Equation 1, 2 and 3 of [10].

$\CIRCLE$ \- Rainfall measuring nodes; $\Uparrow$ \- Water-level measuring
node

Figure 2: Actual locations of spatially distributed sensing nodes deployed in
the pilot study area

A rainfall-runoff data set collected from the storm event in the $05^{th}$
May, $2015$ between $21:12$ Hours to $23:52$ Hours, comprising of $302$
observations (Data Set 1) is used to train the model parameters. Collected
pressure head data from the storm event in the $22^{nd}$ September, $2015$
between $21:40$ Hours to $23:40$ Hours, comprising of $273$ observations (Data
Set 2) are used to validate the models. Each observation used for training and
validation purpose is comprises of rainfall data obtained from three spatially
distributed location and pressure head data obtained at the exit of the study
area. In Data set 1 (training data) and Data set 2 (validation data),
collected data are taken at $30$ seconds interval. Hence, the minimum time
difference between each observation $\delta^{0}_{pred}$ is $30$ seconds in
Data Set 1 and Data Set 2. The proposed model is trained with Data Set 1 and
validated using Data Set 2 without changing the dimension of each data point.
Next, we normalise the input output dataset to develop a dimensionless
training and validation data set. The models are labelled with a parenthesis
$(N)$ when dimensionless input and output data set used to calculate the
performance measure. The model performance is validated through statistical
and hydrological performance measures such as Root Mean Square Error (RMSE),
Coefficient of Efficiency (CE), Volumetric Error (VE) and correlation
coefficient of determination (R) between observed data and simulated data. The
model performance criteria are given by,

${\begin{array}[]{l}RMSE=\sqrt[2]{{\frac{{\sum\limits_{k=1}^{N}{{{({y_{k}}-{{\hat{y}}_{k}})}^{2}}}}}{N}}}\\\
{CE}=\left({1-\frac{F}{{{F_{0}}}}}\right);\\\
VE=\left({\frac{{\sum{{y_{k}}-\sum{{{\hat{y}}_{k}}}}}}{{\sum{{y_{k}}}}}}\right)\times
100\\\
R=\frac{{\sum\limits_{k=1}^{N}{({y_{k}}-\bar{y})({{\hat{y}}_{k}}-{{\bar{y}}_{p}})}}}{{\sqrt[2]{{\sum\limits_{k=1}^{N}{{{({y_{k}}-\bar{y})}^{2}}}}}\sqrt[2]{{\sum\limits_{k=1}^{N}{{{({{\hat{y}}_{k}}-{{\bar{y}}_{p}})}^{2}}}}}}};{{\bar{y}}_{p}}=\frac{1}{N}\sum\limits_{k=1}^{N}{{{\hat{y}}_{k}}}\end{array}}$
(23)

where ${y_{k}}$ is the observed pressure head (in mm), ${{\hat{y}}_{k}}$ is
the predicted pressure head from the model. ${\bar{y}}$ and ${{\bar{y}}_{p}}$
are the mean of observed and predicted data respectively. The RMSE performance
measure shows the algorithm capability for predicting the observed data. Lower
values of RMSE are indicating better fit with the observed data. The $CE$
performance measure indicates the quality of the model verification data taken
for different lengths over different time interval. The correlation
coefficient ($R$) is used to check the goodness-of-fit of the model.

##  VII Results and Discussion

We discuss two types of errors: (i) prediction error i.e. error observed in
the next immediate observation, and (ii) time ahead runoff prediction error
i.e. error observed in a particular time ahead runoff prediction scheme.

###  VII-A Prediction error comparison

TABLE I: Prediction error comparison for all the different algorithm based $M_{4}$ model in both dimensional and dimensionless input-output setting | Training | Validation
---|---|---
Algorithm | RMSE |  VE |  CE |  R | RMSE |  VE |  CE |  R
$GK$ | 0.00 |  0.00 |  1.0 |  1.00 | 0.00 |  0.00 |  1.0 |  1.0
$GK(N)$ | 0.00 |  0.002 |  0.99 |  0.99 | 0.00 |  0.002 |  0.99 |  1
$FCM$ | 0.04 |  0.013 |  1 |  1 | 0.49 |  0.05 |  1 |  1
$FCM(N)$ | 0.02 |  1.09 |  0.99 |  1 | 0.02 |  2.29 |  0.99 |  1
$SC$ | 0.03 |  0.01 |  1 |  1 | 0.06 |  0.01 |  1 |  1
$SC(N)$ | 0.01 |  1.12 |  1 |  1 | 0.04 |  2.96 |  0.98 |  0.99

GK, FCM and SC algorithms based $M_{4}$ TS Fuzzy models are constructed for
both training (Data Set 1) and validation (Data Set 2). Table I, specifies the
performances of $M_{4}$. The performance of the $M_{4}$ model is listed in the
Table I for both dimensional and dimensionless input of the Data Set 2. The
performances of six different cluster validation indices as presented in
Section IV-C have given minimum values at $3$ clusters. Hence, there are $3$
optimal number of clusters (rules) present in the model data set of $M_{4}$.
Error measures for the six $M4$ models are presented in Table I for both the
training and validation data. In Table I, the row that contains dimensionless
input are identifiable by the ’(N)’ beside the algorithm (e.g. $GK(N)$).

The FCM algorithm based $M_{4}$ model is generated by MATLAB "genfis3"
function. The optimal number of rules in the FCM based Fuzzy structure are
similar to the optimal number of rules, obtained by GK validation indices.
Table I indicates that FCM algorithm based model is not as good as GK
algorithm while estimating the runoff data. The SC algorithm based $M_{4}$
model is generated by MATLAB "genfis2" function where number of rules in the
model structure is similar to the optimal number of rules (provided by GK
validation indices). Table I indicates that SC based TS Fuzzy model has a
better performance than that of FCM algorithm but it is not as good as the GK
algorithm based $M_{4}$ model.

Figure 3: Observed vs estimated pressure head for GK based $M_{4}$ model
Figure 4: Observed vs estimated pressure head for GK based $M_{4}$ model
Figure 5: Observed vs estimated pressure head for FCM based $M_{4}$ model
Figure 6: Observed vs estimated pressure head for FCM based $M_{4}$ model
Figure 7: Observed vs estimated pressure head for SC based $M_{4}$ model
Figure 8: Observed vs estimated pressure head for SC based $M_{4}$ model

The GK based $M_{4}$ model performances for the training dataset and
validation dataset are shown in the Figure 3 and Figure 4. The training and
testing dataset performances by SC based $M_{4}$ model are shown in the Figure
7 and Figure 8 respectively. The FCM based $M_{4}$ model performances for the
training and testing dataset are shown in the Figure 5 and Figure 6
respectively.

###  VII-B Time ahead runoff prediction error comparison

This section presents an error performance comparison of the algorithmic
performances of the $M_{4}$ model when the data is predicted according to
three different prediction ahead time i.e. three different schemes selected
from the set $\delta_{pred}$. We demonstrate a comparison between
$\delta^{1}_{pred}=60$ seconds, $\delta^{4}_{pred}=150$ seconds and
$\delta^{9}_{pred}=300$ seconds in Table II. We should note that Table I
represents the errors obtained in the time ahead prediction scheme
$\delta^{0}_{pred}=30$ seconds.

TABLE II: Error comparison for three different prediction ahead time $\delta^{1}_{pred}$ | Training | Validation
---|---|---
| RMSE |  VE |  CE |  R | RMSE |  VE |  CE |  R
$GK$ | 0.00 |  0.00 |  1 |  1 | 8.99 |  0.3 |  0.99 |  0.99
$FCM$ | 0.93 |  0.35 |  0.99 |  0.99 | 9.42 |  1.07 |  0.99 |  0.99
$SC$ | 0.33 |  0.05 |  1 |  1 | 9.25 |  0.31 |  0.99 |  0.99
$\delta^{4}_{pred}$ | Training | Validation
| RMSE |  VE |  CE |  R | RMSE |  VE |  CE |  R
$GK$ | 0.00 |  0.00 |  1 |  1 | 12 |  0.61 |  0.99 |  0.99
$FCM$ | 0.26 |  0.11 |  1 |  1 | 13.01 |  0.88 |  0.98 |  0.99
$SC$ | 0.07 |  0.02 |  1 |  1 | 12.94 |  0.63 |  0.98 |  0.99
$\delta^{9}_{pred}$ | Training | Validation
| RMSE |  VE |  CE |  R | RMSE |  VE |  CE |  R
$GK$ | 0.00 |  0.00 |  1 |  1 | 16.05 |  1.12 |  0.97 |  0.99
$FCM$ | 0.014 |  0.004 |  1 |  1 | 18.06 |  1.23 |  0.96 |  0.99
$SC$ | 0.006 |  0.002 |  1 |  1 | 18.05 |  1.21 |  0.96 |  0.99

Table II shows that the RMSE error of the $M_{4}$ model for GK algorithm has a
$\pm 16\%$ tolerance while predicting the runoff 10 times (300s/30s) earlier.
The RMSE error of the $M_{4}$ model for FCM and SC algorithms have a $\pm
18\%$ tolerance for the same prediction time interval. Figure 9 is showing the
estimated pressure head before $5$ minutes of observed runoff for FCM based
$M_{4}$ model that has the highest RMSE error among rest of the models.

Figure 9: Estimated pressure head before $5$ minutes of observed runoff for
FCM based $M_{4}$ model

##  VIII Conclusion and Future Work

This paper presents a GK clustering algorithm based TS fuzzy model to predict
runoff due to rainfall using read world observations. Further, a comparative
study between FCM and SC clustering algorithms is developed with various error
measures. The proposed model is trained and validated with a different set of
observations recorded on different days. Validation results show that the GK
algorithm based $M_{4}$ model performs significantly better than FCM and SC
algorithms. GK algorithm also performs better in the time ahead runoff
prediction error that increases with $i_{N}$ (the superscript of
$\delta^{i_{N}}_{pred}$).

Proposed models have been trained and validated over a limited number of
rainfall-runoff data due to the catchment geography and climate
characteristics. As future work, the methodologies described in this paper may
be applied to those catchments where rainfall-runoff data is available for a
longer duration of time to predict the surface runoff due to rainfall. Other
hydrological inputs (e.g. Antecedent moisture conditions) and/or an increasing
number of spatially distributed rainfall stations may be integrated with the
input data set for a better prediction of the output.

## References

  * [1] S. N. Kuiry, D. Sen, and P. D. Bates, “Coupled 1d–quasi-2d flood inundation model with unstructured grids,” _Journal of Hydraulic Engineering_ , vol. 136, no. 8, pp. 493–506, 2010.
  * [2] A. P. Jacquin and A. Y. Shamseldin, “Development of rainfall–runoff models using takagi–sugeno fuzzy inference systems,” _Journal of Hydrology_ , vol. 329, no. 1-2, pp. 154–173, 2006.
  * [3] Y. Morales, M. Querales, H. Rosas, H. Allende-Cid, and R. Salas, “A self-identification neuro-fuzzy inference framework for modeling rainfall-runoff in a chilean watershed,” _Journal of Hydrology_ , vol. 594, p. 125910, 2021.
  * [4] A. K. Lohani, N. Goel, and K. Bhatia, “Development of fuzzy logic based real time flood forecasting system for river narmada in central india,” in _International conference on innovation advances and implementation of flood forecasting technology_ , vol. 17, 2005, pp. 1–10.
  * [5] A. Lohani, N. Goel, and K. Bhatia, “Takagi–sugeno fuzzy inference system for modeling stage–discharge relationship,” _Journal of Hydrology_ , vol. 331, no. 1, pp. 146–160, 2006.
  * [6] Y. Hundecha, A. Bardossy, and H.-W. WERNER, “Development of a fuzzy logic-based rainfall-runoff model,” _Hydrological Sciences Journal_ , vol. 46, no. 3, pp. 363–376, 2001.
  * [7] J. shing Roger Jang, “Anfis: Adaptive-network-based fuzzy inference system,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. 23, pp. 665–685, 1993.
  * [8] T. K. Chang, A. Talei, S. Alaghmand, and M. P.-L. Ooi, “Choice of rainfall inputs for event-based rainfall-runoff modeling in a catchment with multiple rainfall stations using data-driven techniques,” _Journal of Hydrology_ , vol. 545, pp. 100–108, 2017.
  * [9] R. Tabbussum and A. Q. Dar, “Performance evaluation of artificial intelligence paradigms—artificial neural networks, fuzzy logic, and adaptive neuro-fuzzy inference system for flood prediction,” _Environmental Science and Pollution Research_ , vol. 28, no. 20, pp. 25 265–25 282, 2021.
  * [10] S. Dey, N. Bhattacharya, S. Chakrabarti, and D. Sen, “Real-time ogc compliant online data monitoring and acquisition network for management of hydro-meteorological hazards,” _ISH Journal of Hydraulic Engineering_ , pp. 1–10, 2016.
  * [11] A. K. Lohani, N. Goel, and K. Bhatia, “Comparative study of neural network, fuzzy logic and linear transfer function techniques in daily rainfall-runoff modelling under different input domains,” _Hydrological Processes_ , vol. 25, no. 2, pp. 175–193, 2011.
  * [12] C. W. Dawson and R. Wilby, “An artificial neural network approach to rainfall-runoff modelling,” _Hydrological Sciences Journal_ , vol. 43, no. 1, pp. 47–66, 1998.
  * [13] G. Tayfur and V. P. Singh, “Ann and fuzzy logic models for simulating event-based rainfall-runoff,” _Journal of hydraulic engineering_ , vol. 132, no. 12, pp. 1321–1330, 2006.
  * [14] F.-J. Chang and Y.-C. Chen, “A counterpropagation fuzzy-neural network modeling approach to real time streamflow prediction,” _J. Hydrol._ , vol. 245, no. 1-4, pp. 153–164, 2001.
  * [15] P. Nayak, K. Sudheer, and K. Ramasastri, “Fuzzy computing based rainfall–runoff model for real time flood forecasting,” _Hydrological Processes_ , vol. 19, no. 4, pp. 955–968, 2005.
  * [16] P. C. Nayak, K. Sudheer, D. Rangan, and K. Ramasastri, “A neuro-fuzzy computing technique for modeling hydrological time series,” _Journal of Hydrology_ , vol. 291, no. 1, pp. 52–66, 2004.
  * [17] P. Nayak, K. Sudheer, D. Rangan, and K. Ramasastri, “Short-term flood forecasting with a neurofuzzy model,” _Water Resources Research_ , vol. 41, no. 4, 2005.
  * [18] A. Mukerji, C. Chatterjee, and N. S. Raghuwanshi, “Flood forecasting using ann, neuro-fuzzy, and neuro-ga models,” _Journal of Hydrologic Engineering_ , vol. 14, no. 6, pp. 647–652, 2009.
  * [19] R. Remesan, M. A. Shamim, D. Han, and J. Mathew, “Runoff prediction using an integrated hybrid modelling scheme,” _Journal of Hydrology_ , vol. 372, no. 1, pp. 48–60, 2009.
  * [20] N. K. Kasabov and Q. Song, “Denfis: dynamic evolving neural-fuzzy inference system and its application for time-series prediction,” _IEEE transactions on Fuzzy Systems_ , vol. 10, no. 2, pp. 144–154, 2002.
  * [21] S. A. Akrami, V. Nourani, and S. Hakim, “Development of nonlinear model based on wavelet-anfis for rainfall forecasting at klang gates dam,” _Water resources management_ , vol. 28, no. 10, pp. 2999–3018, 2014.
  * [22] A. Talei and L. H. Chua, “Influence of lag time on event-based rainfall–runoff modeling using the data driven approach,” _Journal of hydrology_ , vol. 438, pp. 223–233, 2012.
  * [23] R. Govindraju, “Artificial neural networks in hydrology: Ii. hydrological applications,” _J. Hydrol. Eng_ , vol. 5, no. 2, pp. 124–137, 2000.
  * [24] L. H. Chua, T. S. Wong, and L. Sriramula, “Comparison between kinematic wave and artificial neural network models in event-based runoff simulation for an overland plane,” _Journal of hydrology_ , vol. 357, no. 3, pp. 337–348, 2008\.
  * [25] A. Talei, L. H. C. Chua, and C. Quek, “A novel application of a neuro-fuzzy computational technique in event-based rainfall–runoff modeling,” _Expert Systems with Applications_ , vol. 37, no. 12, pp. 7456–7468, 2010\.
  * [26] A. Talei, L. H. C. Chua, C. Quek, and P.-E. Jansson, “Runoff forecasting using a takagi–sugeno neuro-fuzzy model with online learning,” _Journal of hydrology_ , vol. 488, pp. 17–32, 2013.
  * [27] P. Nayak, K. Sudheer, and S. Jain, “Rainfall-runoff modeling through hybrid intelligent system,” _Water Resources Research_ , vol. 43, no. 7, 2007.
  * [28] T. Takagi and M. Sugeno, “Fuzzy identification of systems and its applications to modeling and control,” _IEEE transactions on systems, man, and cybernetics_ , no. 1, pp. 116–132, 1985.
  * [29] S. L. Chiu, “Fuzzy model identification based on cluster estimation,” _Journal of Intelligent & fuzzy systems_, vol. 2, no. 3, pp. 267–278, 1994\.
  * [30] J. Abonyi, R. Babuska, and F. Szeifert, “Modified gath-geva fuzzy clustering for identification of takagi-sugeno fuzzy models,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 32, no. 5, pp. 612–621, 2002.
  * [31] T. Dam and A. K. Deb, “Interval type-2 modified fuzzy c-regression model clustering algorithm in ts fuzzy model identification,” in _Fuzzy Systems (FUZZ-IEEE), 2016 IEEE International Conference on_. IEEE, 2016, pp. 1671–1676.
  * [32] ——, “Ts fuzzy model identification by a novel objective function based fuzzy clustering algorithm,” in _Computational Intelligence in Ensemble Learning (CIEL), 2014 IEEE Symposium on_. IEEE, 2014, pp. 1–7.
  * [33] R. Babuka, P. Van der Veen, and U. Kaymak, “Improved covariance estimation for gustafson-kessel clustering,” in _Fuzzy Systems, 2002. FUZZ-IEEE’02. Proceedings of the 2002 IEEE International Conference on_ , vol. 2. IEEE, 2002, pp. 1081–1085.
  * [34] X. L. Xie and G. Beni, “A validity measure for fuzzy clustering,” _IEEE Transactions on pattern analysis and machine intelligence_ , vol. 13, no. 8, pp. 841–847, 1991.
  * [35] J. C. Bezdek, “Numerical taxonomy with fuzzy sets,” _Journal of Mathematical Biology_ , vol. 1, no. 1, pp. 57–71, 1974.
  * [36] R. N. Dave, “Validating fuzzy partitions obtained through c-shells clustering,” _Pattern Recognition Letters_ , vol. 17, no. 6, pp. 613–623, 1996.
  * [37] A. M. Bensaid, L. O. Hall, J. C. Bezdek, L. P. Clarke, M. L. Silbiger, J. A. Arrington, and R. F. Murtagh, “Validity-guided (re) clustering with applications to image segmentation,” _IEEE Transactions on Fuzzy Systems_ , vol. 4, no. 2, pp. 112–123, 1996.
  * [38] T. Dam and A. K. Deb, “A clustering algorithm based ts fuzzy model for tracking dynamical system data,” _Journal of the Franklin Institute_ , vol. 354, no. 13, pp. 5617–5645, 2017.
  * [39] T. Dam and A. Deb, “Block sparse representations in modified fuzzy c-regression model clustering algorithm for ts fuzzy model identification,” in _2015 IEEE Symposium Series on Computational Intelligence_. IEEE, 2015, pp. 1687–1694.
  * [40] T. Dam and A. K. Deb, “Interval type-2 recursive fuzzy c-means clustering algorithm in the ts fuzzy model identification,” in _2015 IEEE Symposium Series on Computational Intelligence_. IEEE, 2015, pp. 22–29.


# A specific model of Hilbert geometry on the unit disc

Charalampos Charitos, Ioannis Papadoperakis  and Georgios
Tsapogas111Corresponding author
Agricultural University of Athens

###### Abstract

A new metric on the open 2-dimensional unit disk is defined making it a
geodesically complete metric space whose geodesic lines are precisely the
Euclidean straight lines. Moreover, it is shown that the unit disk with this
new metric is not isometric to any hyperbolic model of constant negative
curvature, nor to any convex domain in $\mathbb{R}^{2}$ equipped with its
Hilbert metric.
2020 Mathematics Subject Classification: 51F99

##  1 Introduction

On the 8th August 1900, at the Second International Congress of Mathematics
held in Paris, David Hilbert delivered a lecture entitled “The future problems
of mathematics”, in which he presented a collection of open problems. The
fourth problem of the list can be stated as follows: If $\Omega$ is a convex
subset of a Euclidean space, find a characterization of all metrics on
$\Omega$ for which the Euclidean lines are geodesics. We can put additional
conditions on these metrics on $\Omega$ by requiring geodesic completeness and
Euclidean lines being the unique geodesics. These geometries with the extra
requirements are of particular interest and they have been studied
extensively.

Before Hilbert, Beltrami in [1] had already shown that the unit disc in the
plane, with the Euclidean chords taken as geodesics of infinite length, is a
model of the hyperbolic geometry. However, Beltrami did not give a formula for
this distance, and this led Klein in [4] to express the distance in the unit
disc in terms of the cross radio. Hilbert’s fourth problem became a very
active research area and it was gradually realized that the discovery of all
metrics satisfying Hilbert’s problem was not plausible. Consequently, each
metric resolving Hilbert’s problem defines a new geometry worth to be studied.
A very important class of such metrics, defined by means of the cross ratio,
are referred to as Hilbert metrics and play a central role in this research
area.

Among the prominent mathematicians worked on the Hilbert’s fourth problem, it
is worthy to mention Busemann and Pogorelov, see for instance [2], [7], [8].
The ideas of the latter to solve Hilbert’s fourth problem came from Busemann,
who introduced integral geometry techniques to approach Hilbert’s problem.
Busemann’s idea was to consider for every two points $x$ and $y$ in a convex
subset $\Omega$ of the real projective space $RP^{n},$ the unique geodesic
segment $[x,y]$ joining these points, and the subset of hyperplanes of
$RP^{n}$ intersecting $[x,y]$ equipped with a non negative measure having
specific properties. In dimension $2,$ Pogorelov’s solution consisted in
showing that every distance between $x$ and $y$ satisfying Hilbert’s problem
is given by a metric $d(x,y)$ constructed with the help of the measure
constructed by Busemann on the subset of hyperplanes mentioned above. There
are generalizations of Pogorelov’s theorem in greater dimensions and one may
see in [6] a detailed discussion on Hilbert’s fourth problem.

However, Pogorelov’s approach is very general and does not allow for further
study of the geometry of these metrics. On the contrary, in the present work a
concrete new metric satisfying Hilbert’s problem is defined without the use of
cross ratio and its geometry is studied. More precisely, it is shown that this
metric makes the open unit disk a geodesically complete metric space whose
geodesics are of infinite length and are precisely the Euclidean lines.
Moreover it is shown that this metric space is not isometric to any hyperbolic
model of constant curvature nor to any convex domain in the plane equipped
with its Hilbert metric. Finally, the natural Euclidean boundary of the unit
disk is shown to coincide with the visual boundary with respect to the new
metric, namely, with the set of equivalence classes of asymptotic geodesic
rays.

##  2 Definitions

Let $f:\left(-1,1\right)\longrightarrow\mathbb{R}$ be the function

$f(t)=\frac{t}{1-|t|}$

$\scriptstyle 1$$\scriptstyle 0$$\scriptstyle-1$

and define a metric $d_{I}$ on the interval $\left(-1,1\right)$ by

$d_{I}\left(s,t\right)=\left|f\left(s\right)-f\left(t\right)\right|.$

Clearly,

$d_{I}\left(t,0\right)=d_{I}\left(0,t\right)=\left|\frac{t}{1-|t|}\right|=\frac{1}{1-|t|}-1,\mathrm{\
\ }\forall t\in(-1,1)$

and

$d_{I}\left(s,t\right)=f\left(|s|\right)+f\left(|t|\right)\mathrm{\ if\
}st<0.$

Moreover, this metric has the following two properties:

  1. P1:

For $s,t\in(-1,1),$ $d_{I}\left(s,t\right)\longrightarrow\infty$ as
$t\longrightarrow-1$ or $1.$

  2. P2:

For $q,s,t\in(-1,1)$ with $q<s<t$ we have
$d_{I}\left(q,s\right)+d_{I}\left(s,t\right)=d_{I}\left(q,t\right).$

We now define a metric $D$ on the open 2-dimensional unit disk
$\mathbb{D}^{2}$ as follows: consider the unit disk in the $xy-$plane and
identify each ray with an angle $\theta\in\left[0,2\pi\right].$ For each such
ray (eg. angle $\theta$) denote by $\Delta_{\theta}$ the diameter determined
by that ray. Observe that the Euclidean length of $\Delta_{\theta}$ is $2.$
For any point $z$ in the disk, denote by $z_{\theta}$ its projection to the
diameter $\Delta_{\theta}.$ Let $x,y$ be two points in the disk.

For each $\theta\in\left[0,2\pi\right],$ set
$d_{\theta}\left(x,y\right):=d_{I}\left(x_{\theta},y_{\theta}\right)$ where
the projection points $x_{\theta},y_{\theta}$ in $\Delta\theta$ are identified
with the corresponding points in the interval $(-1,1).$ Define

$D\left(x,y\right):=\frac{1}{2}\int_{0}^{2\pi}d_{\theta}\left(x,y\right)d\theta.$
(1)

Since, for every $\theta\in[0,\pi]$ the diameters $\Delta_{\theta}$ and
$\Delta_{\theta+\pi}$ coincide, we have
$d_{\theta}\left(x,y\right)=d_{\pi+\theta}\left(x,y\right)$ for any two points
$x,y.$ It follows that

$D\left(x,y\right)=\frac{1}{2}\int_{0}^{2\pi}d_{\theta}\left(x,y\right)d\theta=\int_{0}^{\pi}d_{\theta}\left(x,y\right)d\theta$
(2)

######  Lemma 1 (Triangle Inequality).

Let $x,y,z$ be three points in the disk.
(a) If the Euclidean segment $[x,z]$ in the unit disk contains the point $y$
then

$D(x,z)=D(x,y)+D(y,z).$

(b) If $y$ is a point not contained in $[x,z]$ then

$D(x,z)<D(x,y)+D(y,z).$

###### Proof.

(a) The projection $y_{\theta}$ will be in the interior of the segment
$\left[x_{\theta},z_{\theta}\right]\subset\Delta_{\theta}$ for all directions
$\theta$ except in the case the direction $\theta$ is perpendicular to the
segment $[x,z]$ (in which case $x_{\theta}\equiv z_{\theta}$). By Property P2

$d_{\theta}\left(x,y\right)+d_{\theta}\left(y,z\right)=d_{\theta}\left(x,z\right)$

which shows the desired equality.

(b) Set $\theta_{yx}$ (resp. $\theta_{yz}$) to be the angle which, viewed as a
ray, is perpendicular to the Euclidean segment $[y,x]$ (resp. $[y,z]$). We may
assume that $\theta_{yx}<\theta_{yz}.$ Then for every
$\theta\notin[\theta_{yx},\theta_{yz}]$ the point $y_{\theta}$ is contained in
the segment $[x_{\theta},z_{\theta}]$ so that by Property P2 we have
$d_{\theta}(x,z)=d_{\theta}(x,y)+d_{\theta}(y,z).$
$\scriptstyle\theta{yx}$$\scriptstyle\theta{yz}$$\color[rgb]{1,0,0}\scriptstyle\Delta_{\theta}$$\color[rgb]{1,0,0}z_{\theta}$$\color[rgb]{1,0,0}x_{\theta}$$\color[rgb]{1,0,0}y_{\theta}$$z$$x$$\scriptstyle
y$
For $\theta\in[\theta_{yx},\theta_{yz}]$ the point $y_{\theta}$ is not
contained in the segment $[x_{\theta},z_{\theta}]$ which implies that

$d_{\theta}(y,z)=d_{\theta}(y,x)+d_{\theta}(x,z)\textrm{\ \ or,\ \
}d_{\theta}(y,x)=d_{\theta}(y,z)+d_{\theta}(z,x)$

depending on whether $x_{\theta}\in[y_{\theta},z_{\theta}]$ or,
$z_{\theta}\in[y_{\theta},x_{\theta}].$ In the first case we have

$d_{\theta}(x,z)<d_{\theta}(y,z)\Longrightarrow d_{\theta}(x,z)\lneqq
d_{\theta}(y,z)+d_{\theta}(x,y)$

and the same strict inequality follows in the second case. This completes the
proof of part (b). ∎

######  Lemma 2.

The function
$D\left(x,y\right)=\int_{0}^{\pi}d_{\theta}\left(x,y\right)d\theta$ defines a
metric on the open unit disk.

###### Proof.

To complete the proof that $D$ is a metric we only need to show that the
integral

$\int_{0}^{\pi}d_{\theta}\left(x,y\right)d\theta$

is finite. By Lemma 1, it suffices to show that for any $x$ in the unit disk
the integral $\int_{0}^{\pi}d_{\theta}\left(O,x\right)d\theta$ is finite.
Recall that the point $O$ in the diameter $\Delta_{\theta}$ is identified with
$0\in(-1,1)$ and observe that by definition of the function $f$ the
orientation of the diameter $\Delta_{\theta}$ is irrelevant. Let $\|\cdot\|$
denote Euclidean length. Clearly,

$d_{I}(0,x_{\theta})\leq d_{I}(0,\|Ox\|)$

We then have

$D(O,x)=\int_{0}^{\pi}d_{\theta}(O,x)d\theta=\int_{0}^{\pi}d_{I}(0,x_{\theta})d\theta\leq\int_{0}^{\pi}d_{I}(0,\|Ox\|)d\theta=\frac{\|Ox\|}{1-\|Ox\|}\pi.$

∎

It is well known that a curve with endpoints $x,z$ is a geodesic segment with
respect to a metric $d$ if and only if for every $y$ in the curve we have
$d(x,y)+d(y,z)=d(x,z).$ It follows, by part (a) of Lemma 1 above, that
Euclidean lines in the unit disk are geodesics with respect to the metric $D$
and part (b) shows that only the Euclidean lines are geodesics with respect to
$D.$ Hence, we have the following

######  Proposition 3.

The metric space $\left(\mathbb{D}^{2},D\right)$ is a geodesic metric space
whose geodesics are precisely the Euclidean lines in $\mathbb{D}^{2}.$

The following properties follow from the definition of the metric $D.$

######  Proposition 4.

(a) Every Euclidean rotation
$R_{\phi}:\mathbb{D}^{2}\longrightarrow\mathbb{D}^{2},\phi\in[0,2\pi],$
centered at the origin is an isometry of the metric space
$\left(\mathbb{D}^{2},D\right).$
(b) Every Euclidean reflection
$Q_{\phi}:\mathbb{D}^{2}\longrightarrow\mathbb{D}^{2}$ with respect to a line
forming an angle $\phi\in[0,2\pi]$ with the $x-$axis is an isometry of the
metric space $\left(\mathbb{D}^{2},D\right).$

###### Proof.

As
$\int_{0}^{2\pi}d_{\theta}\left(x,y\right)d\theta=\int_{\phi}^{2\pi+\phi}d_{\theta}\left(x,y\right)d\theta$
part (a) follows.

For (b), it suffices, by (a), to show that the Euclidean reflection $R_{0}$
with respect to the $x-$ axis is an isometry. Clearly, for arbitrary
$x,y\in\mathbb{D}^{2}$ and for every $\theta\in[0,2\pi]$ we have

$d_{\theta}\left(x,y\right)=d_{2\pi-\theta}\left(R_{0}(x),R_{0}(y)\right)$

which implies that $D(x,y)=D\left(R_{0}(x),R_{0}(y)\right).$ ∎

We now proceed to show that the metric space $\left(\mathbb{D}^{2},D\right)$
is geodesically complete, that is, every geodesic segment extends uniquely to
a geodesic line of infinite length.

If $\xi$ a point on the boundary, $d_{\theta}\left(O,\xi\right)$ can be
defined via projections as before and it is a positive real for all $\theta,$
except for a single value in $\left[0,\pi\right).$ Thus, the integral
$\int_{0}^{\pi}d_{\theta}\left(O,\xi\right)d\theta$ makes sense and we have

######  Lemma 5.

If $O$ is the center of the unit disk and $\xi$ a point on the boundary, the
integral $\displaystyle\int_{0}^{\pi}d_{\theta}\left(O,\xi\right)d\theta$ is
not bounded.

###### Proof.

As above, the point $O$ in the diameter $\Delta_{\theta}$ is identified with
$0\in(-1,1)$ and, if $\theta_{\xi}$ is the angle determined by $\xi,$ for all
$\theta\neq\theta_{\xi}$ we have $\displaystyle d_{\theta}\left(O,\xi\right)$
$\displaystyle=d_{I}\left(0,\|O\xi\|\cos\left(\theta-\theta_{\xi}\right)\right)=\left|f\left(\cos\left(\theta-\theta_{\xi}\right)\right)\right|$
$\displaystyle=\left|\frac{\cos\left(\theta-\theta_{\xi}\right)}{1-\left|\cos\left(\theta-\theta_{\xi}\right)\right|}\right|=\frac{1}{1-\left|\cos\left(\theta-\theta_{\xi}\right)\right|}-1.$
$\xi$$\color[rgb]{1,0,0}\scriptstyle\Delta_{\theta}$θ$\color[rgb]{1,0,0}\xi_{\theta}$$\scriptstyle
O$$\theta_{\xi}$
Then the integral $\int_{0}^{\pi}d_{\theta}\left(O,\xi\right)d\theta$ equals

$\displaystyle\hskip-11.38109pt\frac{1}{2}\int_{0}^{2\pi}$
$\displaystyle\left[\frac{1}{1-\left|\cos\left(\theta-\theta_{\xi}\right)\right|}-1\right]d\theta=\frac{1}{2}\int_{0}^{2\pi}\left[\frac{1}{1-\left|\cos\theta\right|}-1\right]d\theta=$
$\displaystyle\hskip
62.59605pt=\frac{1}{2}4\int_{0}^{\pi/2}\left[\frac{1}{1-\cos\theta}-1\right]d\theta>2\int_{0}^{\pi/2}\left[\frac{1}{\theta}-1\right]d\theta=\infty.$

∎

The above Lemma permits us to say that the $D-$length of a Euclidean ray (or,
diameter) is infinite. In a similar manner, if $\xi,\eta$ are two points on
the boundary, $d_{\theta}\left(\xi,\eta\right)$ is a positive real for all
$\theta,$ except for two values of $\theta$ in $\left[0,\pi\right).$ Thus, the
integral $\int_{0}^{\pi}d_{\theta}\left(\xi,\eta\right)d\theta$ makes sense
and we have

######  Lemma 6.

If $\xi,\eta$ are two points on the boundary, the integral
$\displaystyle\int_{0}^{\pi}d_{\theta}\left(\xi,\eta\right)d\theta$ is not
bounded.

###### Proof.

By Proposition 4, we may assume that the geodesic line determined by
$\xi,\eta$ is perpendicular to the $x-$axis with intersection point, say, $A.$
It suffices to show that
$\displaystyle\int_{0}^{\pi}d_{\theta}\left(A,\xi\right)d\theta$ is not
bounded.

Let $\theta_{\xi}$ be the angle formed by the $x-$axis and the geodesic ray
joining $O$ with $\xi.$ For all
$\theta\in\left[0,\pi/2\right]\setminus\left\\{\theta_{\xi}\right\\}$ we have

$d_{\theta}\left(O,\xi\right)=d_{\theta}\left(O,A\right)+d_{\theta}\left(A,\xi\right)$

and for all $\theta\in\left[0,\pi/2\right]$ we have

$d_{\theta}\left(O,\xi\right)<d_{\theta}\left(A,\xi\right)<d_{\theta}\left(O,A\right)+d_{\theta}\left(A,\xi\right).$

Therefore, the triangle inequality

$d_{\theta}\left(O,\xi\right)\leq
d_{\theta}\left(O,A\right)+d_{\theta}\left(A,\xi\right)$

holds for all
$\theta\in\left[0,\right]\setminus\left\\{\theta_{\xi}\right\\}.$ By Lemma 5,
$\int_{0}^{\pi}d_{\theta}\left(O,\xi\right)d\theta$ is not bounded and
$\int_{0}^{\pi}d_{\theta}\left(O,A\right)d\theta=D(O,A)$ is a positive real,
thus,

$\displaystyle\int_{0}^{\pi}d_{\theta}\left(A,\xi\right)d\theta$

cannot be bounded. ∎

######  Remark 7.

There exists a large family of metrics making the open unit disk a geodesic
metric space satisfying Propositions 3, 4 and Lemmata 5, 6. In fact, for every
strictly increasing function $g:(-1,1)\rightarrow\mathbb{R}$ which satisfies

$\lim_{t\rightarrow
1}\int_{0}^{2\pi}\left|g\left(t\cos\theta\right)\right|d\theta=\infty$

we may apply the above construction using $g$ instead of $f$ and obtain a
metric on the open unit disk with the above mentioned properties.

##  3 Further properties of $\left(\mathbb{D}^{2},D\right)$

For any $\kappa>0,$ let $\mathbb{H}_{-\kappa^{2}}^{2}$ denote the standard
hyperbolic model of constant negative curvature on the open unit disk with
distance function $d_{\kappa}.$ For a convex domain $U$ in $\mathbb{R}^{2}$
denote by $d_{\mathcal{H}}$ the Hilbert metric for which we refer the reader
to [5, Ch.5, Section 6].

######  Theorem 8.

For all $\kappa>0,$ the metric spaces $\left(\mathbb{D}^{2},D\right)$ and
$\left(\mathbb{H}_{-\kappa^{2}}^{2},d_{\kappa}\right)$ are not isometric.
Moreover, for any convex domain $U$ in $\mathbb{R}^{2}$ equipped with the
Hilbert metric $d_{\mathcal{H}}$ the metric spaces
$\left(\mathbb{D}^{2},D\right)$ and $\left(U,d_{\mathcal{H}}\right)$ are not
isometric.

###### Proof.

Assume $F_{\kappa}:\mathbb{D}^{2}\longrightarrow\mathbb{H}_{-\kappa^{2}}^{2}$
is such an isometry which, by homogeneity of $\mathbb{H}_{-\kappa^{2}}^{2},$
we may assume that $F_{\kappa}$ preserves the center $O$ of the disk.
Moreover, as the image of the $x-$axis under $F_{\kappa}$ is a line in
$\mathbb{H}_{-\kappa^{2}}^{2}$ containing the origin, by composing
$F_{\kappa}$ with a rotation in $\mathbb{H}_{-\kappa^{2}}^{2},$ we may assume
that $F_{\kappa}$ preserves the $x-$axis. We next show that $F_{\kappa}$
necessarily preserves the $y-$axis as well, by showing that the image of the
$y-$axis under $F_{\kappa}$ is a line perpendicular to the $x-$axis. To see
this, for any $x\in\left(0,1\right)$ consider the quadrilateral $XYZW$ where
$X=(x,0),$ $Y=(0,x),$ $Z=\left(-x,0\right)$ and $W=\left(0,-x\right).$
Clearly, $XYZW$ is a square with respect to the metric $D$ and, hence, so must
be its image under $F_{\kappa}.$ The geodesic segments
$\left[F_{\kappa}(X),F_{\kappa}\left(Z\right)\right]$ and
$\left[F_{\kappa}(Y),F_{\kappa}\left(W\right)\right]$ intersect at
$O=F_{\kappa}(O)$ which is the midpoint for both segments. Moreover, these
segments must form a right angle at $O$, otherwise the quadrilateral
$F_{\kappa}(X)F_{\kappa}\left(Y\right)F_{\kappa}(Z)F_{\kappa}\left(W\right)$
cannot be a square.

Define the function $h:\left[0,\infty\right]\longrightarrow\mathbb{R}$ where
$h(b),$ for $b\in\left[0,\infty\right),$ is the $d_{1}-$length of the height
of the right angle hyperbolic triangle in $\mathbb{H}_{-1}^{2}$ with side
lengths equal to $b.$ For $b=\infty,$ $h\left(\infty\right)$ is the
$d_{1}-$length of the height of the right angle ideal hyperbolic triangle in
$\mathbb{H}_{-1}^{2}$ with vertices $O,(1,0)$ and $(0,1).$ By elementary
calculations, $h\left(\infty\right)$ is the $d_{1}-$length of the segment with
endpoints $O$ and $\left(1-\sqrt{2}/2,1-\sqrt{2}/2\right).$ Thus its Euclidean
length is $\sqrt{2}-1$ and

$h\left(\infty\right)=\log\frac{1+\left(\sqrt{2}-1\right)}{1-\left(\sqrt{2}-1\right)}=\log\left(1+\sqrt{2}\right)\approx
0.8813735870$ (3)

For $b\in\left[0,\infty\right)$ let $B$ be the point on the positive $x-$axis
so that $OB$ has hyperbolic length $b$ and denote by $C$ the trace from the
origin of the height of the right angle hyperbolic triangle in
$\mathbb{H}_{-1}^{2}$ with side lengths equal to $b.$ Then the triangle
$\triangle(OCB)$ has a right angle at $C$ and $\widehat{COB}=\pi/4.$ Using the
formula

$\cos\frac{\pi}{4}=\frac{\tanh\left(h\left(b\right)\right)}{\tanh\left(b\right)}$

we find

$h\left(b\right)=\frac{1}{2}\log\frac{1+(\sqrt{2}/2)\tanh
b}{1-(\sqrt{2}/2)\tanh b}$ (4)

We next define an analogous function $h^{D}$ for the metric space
$\left(\mathbb{D}^{2},D\right)$ with a different domain

$h^{D}:\left[0,1\right]\longrightarrow\mathbb{R}$

where for $x\in\left[0,1\right),$ $h^{D}(x)$ is the $D-$length of the height
of the right angle geodesic triangle $T_{x}$ with vertices $O,(x,0)$ and
$\left(0,x\right).$ For $x=1,$ $h^{D}\left(1\right)$ is the length of the
height of the corresponding ideal geodesic triangle. As all rotations of the
unit disk are isometries of $\left(\mathbb{D}^{2},D\right)$ we have

$h^{D}\left(x\right)=D\left(O,\frac{\sqrt{2}}{2}x\right)\mathrm{\ \ for\ all\
}x\in[0,1].$ (5)

We will explicitly compute the $D-$lengths of the heights of the triangles
$T_{\sqrt{2}/2}$ and $T_{1}$ and compare them with the corresponding
$d_{k}-$lengths of the heights of the triangles
$F_{\kappa}\left(T_{\sqrt{2}/2}\right)$ and $F_{\kappa}\left(T_{1}\right).$
The comparison of the lengths of the heights of $T_{1}$ and
$F_{\kappa}\left(T_{1}\right)$ will suffice to reach a contradiction for the
case $\kappa=1.$ The triangles $T_{\sqrt{2}/2}$ and
$F_{\kappa}\left(T_{\sqrt{2}/2}\right)$ are deployed in order to reach a
contradiction for all $\kappa.$

We first compute the $D-$lengths of the heights of the triangles
$T_{\sqrt{2}/2}$ and $T_{1}.$
For the triangle $T_{\sqrt{2}/2},$ by (5) and using the easily verified fact
that the derivative with respect to $\theta$ of the function
$\frac{4\tan^{-1}\left(\frac{1}{\sqrt{3}}\tan\left(\frac{\theta}{2}\right)\right)}{\sqrt{3}}$
is $\frac{2}{2-\cos x},$ we have

$\displaystyle h^{D}\left(\frac{\sqrt{2}}{2}\right)$
$\displaystyle=D\left(O,\frac{\sqrt{2}}{2}\,\,\frac{\sqrt{2}}{2}\right)=D\left(O,\frac{1}{2}\right)=\int_{0}^{\pi}\left[\frac{1}{1-\frac{1}{2}\left|\cos\theta\right|}-1\right]d\theta$
$\displaystyle=\int_{0}^{\pi/2}\left[\frac{1}{1-\frac{1}{2}\cos\theta}-1\right]d\theta+\int_{\pi/2}^{\pi}\left[\frac{1}{1+\frac{1}{2}\cos\theta}-1\right]d\theta$
$\displaystyle=2\left[\frac{4\tan^{-1}\left(\sqrt{3}\tan\left(\frac{\theta}{2}\right)\right)}{\sqrt{3}}-\theta\right]_{\theta=0}^{\theta=\pi/2}$
$\displaystyle=2\left[\frac{8\sqrt{3}-9}{18}\pi\right]=\frac{8\sqrt{3}-9}{9}\pi\hskip
5.69046pt\approx 1.695205651$ (6)

A similar computation, using again the easily verified fact that the
derivative with respect to $\theta$ of the function
$2\sqrt{2}\tan^{-1}\left(\frac{\left(\sqrt{2}+2\right)\tan\left(\frac{\theta}{2}\right)}{\sqrt{2}}\right)$
is $\frac{2+\sqrt{2}}{-(1+\sqrt{2})\cos\theta+\sqrt{2}+2},$ shows

$h^{D}\left(1\right)=D\left(O,\frac{\sqrt{2}}{2}\right)=\frac{3\sqrt{2}-2}{2}\pi\approx
3.522731754$ (7)

The above calculation along with (3) shows that

$h^{D}\left(1\right)=\frac{3\sqrt{2}-2}{2}\pi\neq\log\left(1+\sqrt{2}\right)=h\left(\infty\right)$
(8)

and, thus, the triangles $T_{1}$ and $F_{1}\left(T_{1}\right)$ cannot be
isometric. It follows that $F_{\kappa}$ cannot be an isometry in the case
$\kappa=1.$

Before proceeding with the general case, we compute the $d_{1}-$length of the
height of the triangle $F\left(T_{\sqrt{2}/2}\right).$ This triangle, being
isometric to $T_{\sqrt{2}/2},$ has side lengths
$D\left(O,\frac{\sqrt{2}}{2}\right)$ and, using (4), its height has
$d_{1}-$length

$\displaystyle h\left(D\left(O,\frac{\sqrt{2}}{2}\right)\right)$
$\displaystyle=h\left(\frac{3\sqrt{2}-2}{2}\pi\right)$
$\displaystyle=\frac{1}{2}\log\frac{1+(\sqrt{2}/2)\tanh\left(\frac{3\sqrt{2}-2}{2}\pi\right)}{1-(\sqrt{2}/2)\tanh\left(\frac{3\sqrt{2}-2}{2}\pi\right)}\approx
0.8789154496$ (9)

We proceed now with the general case. Recall that geodesic lines in
$\mathbb{H}_{-\kappa^{2}}^{2}$ and $\mathbb{H}_{-1}^{2}$ coincide as subsets
of $\mathbb{D}^{2}$ and lengths are multiplied by $\kappa.$ Therefore,
$F_{\kappa}\left(T_{1}\right)=F_{1}\left(T_{1}\right)$ and the
$d_{\kappa}-$length of the height of the triangle
$F_{\kappa}\left(T_{1}\right)$ is equal to $\kappa\,h\left(\infty\right).$ By
(8), it follows that $F_{\kappa}$ can be an isometry only for the model
$\mathbb{H}_{-\kappa^{2}_{0}}^{2}$ where

$\kappa_{0}=\frac{h^{D}\left(1\right)}{h\left(\infty\right)}=\frac{3\sqrt{2}-2}{2\log\left(1+\sqrt{2}\right)}\pi>3.$
(10)

To rule out this last case we compare the $D-$length of the height of the
triangle $T_{\sqrt{2}/2}$ (computed in (6) above) with the
$d_{\kappa_{0}}-$length of the height of the triangle
$F_{\kappa_{0}}\left(T_{\sqrt{2}/2}\right).$ As before, the triangles
$F_{1}\left(T_{\sqrt{2}/2}\right)$ and
$F_{\kappa_{0}}\left(T_{\sqrt{2}/2}\right)$ coincide as sets and the
$d_{\kappa_{0}}-$length of its height is its $d_{1}-$length (computed in (9)
above) multiplied by $\kappa_{0}.$ Therefore, if the triangles
$T_{\sqrt{2}/2}$ and $F_{\kappa_{0}}\left(T_{\sqrt{2}/2}\right)$ were
isometric, $\kappa_{0}$ would have to satisfy

$h^{D}\left(\frac{\sqrt{2}}{2}\right)=\kappa_{0}\,\,h\left(D\left(O,\frac{\sqrt{2}}{2}\right)\right)$

which is impossible because $\kappa_{0}>3$ (see (10)) and the ratio of
$h^{D}\left(\frac{\sqrt{2}}{2}\right)$ and
$h\left(D\left(O,\frac{\sqrt{2}}{2}\right)\right)$ is, by (6) and (9), equal
to

$\frac{h^{D}\left(\frac{\sqrt{2}}{2}\right)}{h\left(D\left(O,\frac{\sqrt{2}}{2}\right)\right)}=\frac{\frac{8\sqrt{3}-9}{9}\pi}{h\left(\frac{3\sqrt{2}-2}{2}\pi\right)}\approx\frac{1.695205651}{0.8789154496}<3.$

This completes the proof for arbitrary curvature.

To see that $\left(\mathbb{D}^{2},D\right)$ is not isometric to any convex
domain $U$ equipped with its Hilbert metric $d_{\mathcal{H}}$, we will use a
result of Busemann and Kelly (see [3, §29.2]) which states the following: let
$U$ be a bounded open convex domain in $\mathbb{R}^{2}.$ Reflections with
respect to all lines in $U$ through one fixed point exist if and only if $U$
is the interior of an ellipse.

The reflection assumption holds for the metric space
$\left(\mathbb{D}^{2},D\right),$ see Proposition 4. If
$\left(\mathbb{D}^{2},D\right)$ were isometric to some
$\left(U,d_{\mathcal{H}}\right)$ then the same reflection assumption would
hold for $U$ and, hence, $U$ would have to be the interior of an ellipse
making $d_{\mathcal{H}}$ the hyperbolic metric. As shown above this cannot be
the case. ∎

We next restrict our attention to geodesic rays, that is, isometric maps
$[0,\infty)\rightarrow\mathbb{D}^{2}.$

######  Definition 9.

Two geodesic rays $r_{1},r_{2}:[0,\infty)\rightarrow\mathbb{D}^{2}$ in the
geodesic metric space $\left(\mathbb{D}^{2},D\right)$ are called asymptotic if
the distance function $t\rightarrow D(r_{1}(t),r_{2}(t))$ is bounded.

######  Remark 10.

Asymptoticity of geodesic rays may be seen as a generalization to arbitrary
metric spaces of parallelism of geodesic rays in Euclidean space. Moreover,
equivalence classes of asymptotic geodesic rays are the tool to define the
visual boundary of a geodesic metric space. It is well known, see for example
[5, Prop. 10.1.4], that two geodesic rays in a geodesic metric space are
asymptotic if and only if their images are at finite Hausdorff distance, a
notion defined below.

######  Definition 11.

For two geodesic rays $r_{1}$ and $r_{2}$ in a geodesic metric space $(X,d)$,
define their Hausdorff distance by

$d_{H}\left(r_{1},r_{2}\right)=\max\left\\{\sup_{x\in\mathop{\rm Im}\nolimits
r_{1}}d\left(x,\mathop{\rm Im}\nolimits r_{2}\right),\sup_{x\in\mathop{\rm
Im}\nolimits r_{2}}d\left(x,\mathop{\rm Im}\nolimits r_{1}\right)\right\\}$
(11)

where the distance of a point $\alpha$ from a set $B$ is
$d\left(\alpha,B\right)=\inf_{\beta\in B}d\left(\alpha,\beta\right).$

As geodesics in $\left(\mathbb{D}^{2},D\right)$ coincide with Euclidean lines,
it is natural to examine whether the natural Euclidean boundary
$\mathbb{S}^{1}$ of $\mathbb{D}^{2}$ coincides with the set of equivalence
classes of asymptotic geodesic rays in $\left(\mathbb{D}^{2},D\right).$
We say that two geodesics rays $r_{1}$ and $r_{2}$ coincide at infinity if, as
Euclidean lines, intersect the same point of
$\mathbb{S}^{1}\equiv\partial\mathbb{D}^{2}.$

######  Theorem 12.

Let $r_{1}$ and $r_{2}$ be two geodesic rays in
$\left(\mathbb{D}^{2},D\right).$ Then $r_{1}$ and $r_{2}$ coincide at infinity
if and only if they are asymptotic.

###### Proof.

The only if direction follows from Lemma 6. In view of the above Remark we
will show that $r_{1}$ and $r_{2}$ coincide at infinity then their images are
at finite Hausdorff distance.

Since rotations around the origin $O$ are isometries (see Proposition 4) we
may assume that the common point at infinity is
$(1,0)\in\partial\mathbb{D}^{2}.$ We will first examine the case where one of
the geodesic rays is the positive $x-$axis and the other one is contained in
the upper half disk forming an angle $\omega\in[0,\pi/2)$ with the $x-$axis at
$(1,0).$ Clearly, the Hausdorff distance is increasing with respect to
$\omega.$ Thus, we may consider a geodesic ray being contained entirely in the
first quadrant forming an angle $\omega\in[\pi/4,\pi/2)$ with the $x-$axis at
$(1,0).$ The general case will then follow easily.

Set $\lambda=\tan\omega$ and for any $x$ satisfying $0<x<1$ consider the point
$B=(x,0)$ on the (geodesic) $x-$axis and the point
$A=\left(x,\lambda(1-x)\right)$ on the other geodesic ray. Set $\theta_{x}$ to
be the angle formed by the segment $OA$ and the $x-$axis. Clearly the distance
$D\left(A,B\right)$ depends on $x$ and it suffices to show that

$\lim_{x\rightarrow 1}D\left(A,B\right)\mathrm{\ is\ bounded.}$ (12)

Recall that for a direction $\Delta_{\theta},$ $\theta\in\left[0,\pi\right]$
and a point $X$ in the interior of the unit disk we denote by $X_{\theta}$ its
projection on the diameter $\Delta_{\theta}.$ For the points $A$ and $B$ we
have

$\left\|OA_{\theta}\right\|=\left\|OA\right\|\cos\left(\theta-\theta_{x}\right)=x\frac{\cos\left(\theta-\theta_{x}\right)}{\cos\theta_{x}}\mathrm{\
\ and\ \ }\left\|OB_{\theta}\right\|=x\cos\theta.$

As $A$ is contained in the first quadrant, for all
$\theta\in\left[0,\frac{\pi}{2}\right]$ we obtain

$d_{\theta}\left(A,B\right)=d_{I}\left(A_{\theta},B_{\theta}\right)=\frac{1}{1-\left\|OA_{\theta}\right\|}-\frac{1}{1-\left\|OB_{\theta}\right\|}=\frac{1}{1-x\frac{\cos\left(\theta-\theta_{x}\right)}{\cos\theta_{x}}}-\frac{1}{1-x\cos\theta}.$
(13)

In a similar manner we find

$d_{\theta}\left(A,B\right)=\begin{cases}\displaystyle\frac{1}{1-x\frac{\cos\left(\theta-\theta_{x}\right)}{\cos\theta_{x}}}+\frac{1}{1+x\cos\theta},&\displaystyle\text{if
}\theta\in\left[\frac{\pi}{2},\frac{\pi}{2}+\theta_{x}\right]\\\\[22.76219pt]
\displaystyle-\frac{1}{1+x\frac{\cos\left(\theta-\theta_{x}\right)}{\cos\theta_{x}}}+\frac{1}{1+x\cos\theta},&\text{if
}\theta\in\left[\frac{\pi}{2}+\theta_{x},\pi\right]\end{cases}$

In view of (12) we will only examine the limit as $x\rightarrow 1$ of the
integral

$\int_{0}^{\theta_{0}}d_{\theta}\left(A,B\right)d\theta=\int_{0}^{\theta_{0}}\left[\frac{1}{1-x\frac{\cos\left(\theta-\theta_{x}\right)}{\cos\theta_{x}}}-\frac{1}{1-x\cos\theta}\right]d\theta$

for sufficiently small $\theta_{0}$ (to be chosen later) because all the above
expressions for $d_{\theta}\left(A,B\right)$ are continuous and bounded on
$\left[\theta_{0},\pi-\theta_{0}\right]$ and the integral
$\int_{\pi-\theta_{0}}^{\pi}d_{\theta}\left(A,B\right)d\theta$ is treated
similarly. By substituting

$\sin\theta_{x}=\frac{\lambda\left(1-x\right)}{\sqrt{x^{2}+\lambda^{2}\left(1-x\right)^{2}}}\mathrm{\
\ and\ \
}\cos\theta_{x}=\frac{x}{\sqrt{x^{2}+\lambda^{2}\left(1-x\right)^{2}}}$

in (13) we obtain

$\displaystyle d_{\theta}\left(A,B\right)$
$\displaystyle=\frac{x\sin\theta_{x}\sin\theta}{\cos\theta_{x}-x\cos\left(\theta-\theta_{x}\right)-x\cos\theta_{x}\cos\theta+x^{2}\cos\left(\theta-\theta_{x}\right)\cos\theta}$
$\displaystyle=\frac{\lambda\left(1-x\right)\sin\theta}{1-x\cos\theta-\lambda\left(1-x\right)\sin\theta-x\cos\theta+x^{2}\cos^{2}\theta+\lambda
x\left(1-x\right)\cos\theta\sin\theta}$
$\displaystyle=\frac{\lambda\left(1-x\right)\sin\theta}{\left(1-x\cos\theta\right)\left(1-x\cos\theta-\lambda\left(1-x\right)\sin\theta\right)}$
(14)

Define

$\Phi\left(\theta\right)=\frac{\sin\theta}{1-x\cos\theta-\lambda\left(1-x\right)\sin\theta}$

and a straightforward calculation shows that

$\Phi^{\prime}\left(\theta\right)=\frac{-x+\cos\theta}{\left(1-x\cos\theta-\lambda\left(1-x\right)\sin\theta\right)^{2}}.$

Therefore, there exists a unique angle
$\omega_{x}\in\left(\theta_{x},\frac{\pi}{2}\right)$ such that

$\cos\omega_{x}=x\Longleftrightarrow\Phi^{\prime}\left(\theta\right)=0.$

Using the equalities $\cos\omega_{x}=x$ and $\sin\omega_{x}=\sqrt{1-x^{2}}$ it
is easily shown that

$\displaystyle\sqrt{1-x}\Phi\left(\omega_{x}\right)$
$\displaystyle=\frac{\sqrt{1-x}\sin\omega_{x}}{1-x\cos\omega_{x}-\lambda\left(1-x\right)\sin\omega_{x}}=\frac{\sqrt{1-x}\sqrt{1-x^{2}}}{1-x^{2}-\lambda\left(1-x\right)\sqrt{1-x^{2}}}$
$\displaystyle=\frac{\sqrt{1+x}\left(1-x\right)}{\left(1-x\right)\left(1+x-\lambda\sqrt{1-x^{2}}\right)}\longrightarrow\frac{\sqrt{2}}{2}\mathrm{\
\ as\ \ }x\rightarrow 1.$ (15)

The choice of $x$ determines both $A$ and $B$ as well as $\theta_{x}$, thus,
we may choose $\theta_{0}$ such that

$\mathrm{for\ \ all\ \
}\theta\leq\theta_{0},\sqrt{1-x}\Phi\left(\omega_{x}\right)\leq 2.$ (16)

As the quantity $\Phi\left(\theta\right)$ attains its maximum at
$\theta=\omega_{x}$ we have, using (14),

$\displaystyle\int_{0}^{\theta_{0}}d_{\theta}\left(A,B\right)d\theta$
$\displaystyle=\int_{0}^{\theta_{0}}\frac{\lambda\left(1-x\right)}{\left(1-x\cos\theta\right)}\Phi\left(\theta\right)d\theta\leq\int_{0}^{\theta_{0}}\frac{\lambda\left(1-x\right)}{\left(1-x\cos\theta\right)}\Phi\left(\omega_{x}\right)d\theta$
$\displaystyle=\int_{0}^{\theta_{0}}\frac{\lambda\sqrt{1-x}}{\left(1-x\cos\theta\right)}\sqrt{1-x}\Phi\left(\omega_{x}\right)d\theta\leq
2\int_{0}^{\theta_{0}}\frac{\lambda\sqrt{1-x}}{\left(1-x\cos\theta\right)}d\theta$
(17)

where the latter inequality follows from (16). It suffices to show that
$\int_{0}^{\theta_{0}}\frac{\sqrt{1-x}}{\left(1-x\cos\theta\right)}d\theta$ is
bounded which follows from the following identity

$\int\frac{1}{1-x\cos\theta}d\theta=\frac{2}{\sqrt{1-x^{2}}}\tan^{-1}\left(\sqrt{\frac{1+x}{1-x}}\tan\frac{\theta}{2}\right)$

and the observation that the range of the inverse tangent function is a
bounded interval:

$\displaystyle\int_{0}^{\theta_{0}}\frac{\sqrt{1-x}}{\left(1-x\cos\theta\right)}d\theta$
$\displaystyle=\frac{\sqrt{1-x}}{\sqrt{1-x^{2}}}\left[\tan^{-1}\left(\sqrt{\frac{1+x}{1-x}}\tan\frac{\theta}{2}\right)\right]_{\theta=0}^{\theta=\theta_{0}}$
$\displaystyle=\frac{1}{\sqrt{1+x}}\tan^{-1}\left(\sqrt{\frac{1+x}{1-x}}\tan\frac{\theta_{0}}{2}\right).$

We now discuss the case of two arbitrary asymptotic geodesic rays $r_{1}$ and
$r_{2}.$ As mentioned at the beginning of the proof, we may assume that the
common boundary point is $(1,0)\in\partial\mathbb{D}^{2}.$ Let
$\omega_{i}\in\left(-\pi/2,\pi/2\right),$ $i=1,2$ be the angle formed by
$r_{i}$ and the $x-$axis. Denote by $r_{x}$ the geodesic ray whose image is
the positive $x-$axis in the unit disk. If both $\omega_{1},\omega_{2}$ are
positive and, say, $\omega_{1}<\omega_{2}$ then
$d_{H}\left(r_{1},r_{2}\right)<d_{H}\left(r_{x},r_{2}\right)$. If
$\omega_{1}\omega_{2}<0$ then a triangle inequality argument asserts that
$d_{H}\left(r_{1},r_{2}\right)\leq
d_{H}\left(r_{x},r_{1}\right)+d_{H}\left(r_{x},r_{1}\right).$ This completes
the proof of the theorem. ∎

Acknowledgments. The authors would like to thank the anonymous referee for
very helpful comments and suggestions which, among other things, improved the
exposition of the proof of Theorem 8.

Conflict of Interest statement: On behalf of all authors, the corresponding
author states that there is no conflict of interest.
Data Availability Statement: Data sharing not applicable to this article as no
datasets were generated or analysed during the current study.

## References

  * [1] E. Beltrami, Saggio di Interpretazione della geometria non-Euclidea. Giornaledi Matematiche vol. VI (1868). Beltrami’s Works, Vol. I, pp. 374–405. French translation: Essai d’interprétation de la géométrie noneucilienne, by J. Hoüel, Annales ENS, 1re Sér. t. 6 (1869), p. 251-288. English translation: Essay on the interpretation of non-Euclidean geometry, by J. Stillwell, in [9] pp. 7–34.
  * [2] H. Busemann, On Hilbert fourth problem, Uspechi Mat. Nauk 21 (1966), no. 1 (127), 155-164.
  * [3] H. Busemann and P.J. Kelly, Projective geometry and projective metrics, Academic Press, New York 1953.
  * [4] F. Klein, Uber die sogenannte Nicht-Euklidische Geometrie (erster Aufsatz). Math. Ann. IV (1871), 573–625. English translation by J. Stillwell in [9].
  * [5] A. Papadopoulos, Metric Spaces, Convexity and Nonpositive Curvature, IRMA Lectures in Mathematics and Theoretical Physics 6, European Mathematical Society, 2013.
  * [6] A. Papadopoulos, Hilbert’s fourth problem, Hanbook of Hilbert Geometry, Irma Lectures in Mathematics and Theoretical Physics 22, A. Papadopoulos and M. Troyanov editors, EMS 2014.
  * [7] A. V. Pogorelov, Hilbert’s fourth problem. Translated by R. A. Silverman, Edited by I. Kra, in cooperation with E. Zaustinskyi, V. H. Winston & Sons, Washington D. C., 1979. Original by Nauka.
  * [8] A. V. Pogorelov, A complete solution of Hilbert’s fourth problem, Sov. Math. 14 (1973), pp. 46–49.
  * [9] J. Stillwell, Sources of hyperbolic geometry, History of Mathematics series, Vol. 10, AMS-LMS, 1996.