text:
\section{Introduction}

To analyze a large-scale dockless bike-sharing system (DBSS) with unusable
bikes, this paper proposes and develops a new computational method by
applying the RG-factorizations of block-structured Markov process in the
closed queueing networks. This method has the strong analytic and
computational ability to deal with the challenging fact that the unusable
bikes are not only randomly generated but also disorderly distributed in
various parking regions of the DBSS. This fact makes the geographical
structure and operations management of DBSSs more complicated so that
analysis of such a DBSS becomes more difficult and challenging from the
stochastic network theory through observing at least three key factors: (1)
The DBSS is a stochastic complex network, and its analysis is always very
difficult. (2) The bikes' failure, removing, repair, redistribution and
reuse processes need a comprehensive study of more complicated stochastic
networks, but so far there has been no existing methodology or results along
such a research line. (3) To remove the unusable bikes and redistribute the
repaired bikes efficiently, two batch policies need to be introduced and
developed from a practical point of view, but no previous work is available.
To overcome these difficulties, this paper establishes the product-form
solution of the closed queueing network by means of not only the nonlinear
matrix equation satisfied by the relative arrival rates but also the
RG-factorizations of block-structured Markov processes. Based on this, we
can provide performance analysis of the DBSS with unusable bikes. We hope
that the methodology and results given in this paper shed light on the study
of more general DBSSs with unusable bikes, and accordingly we can open a
series of promising research in this important direction.

Now, the DBSSs are an important part of smart transportation, which is based
on some key equipments with advanced information and network technologies,
such as GPS, GIS, big data, mobile payment and so forth. Note that the DBSSs
first appeared in 2014 and sprung up in 2016. Up to now, the DBSSs have been
implemented in more than 250 cities in 20 countries over the world, for
example, China, America, France, England, Singapore, India and so on. During
the development of DBSSs, a prominent issue is that there are a great deal
of damaged and unusable bikes stacked in parking regions or roadsides, which
not only negatively decrease the quality of service but also seriously
affect the urban green environment. For example, for the Citi Bike system in
New York in May 2019, there were 11,422 bikes and the maintenance shop
completed 21,353 repairs for unusable bikes, so every bike is repaired 1.87
times on average in that month. See \cite{Cit:2019} for more details.

To deal with a lot of unusable bikes of the DBSS, there are some interesting
practical problems that need to be answered, for example, how to estimate
the effect of unusable bikes on quality of service and further on system
performance? How to develop effective removing policies of unusable bikes
such that system performance can further be improved? However, so far little
work has studied the interesting problems of bike-sharing systems with
unusable bikes. Therefore, the aim of this paper is to develop a stochastic network
model to discuss the bikes' failure, removing, repair, redistribution and
reuse processes under two batch policies, and to further evaluate performance
measures of the DBSS with unusable bikes.

It is necessary and useful to compare the methodology and results given in
this paper with that in the literature of bike-sharing systems. Our analysis
have three key features: (1) A nonlinear matrix equation is set up to determine the relative arrival rates of the virtual closed queueing network, (2) the block-structured Markov process is
used to describe and discuss the Markov system of a node in the closed
queueing network, and (3) a novel maintenance network is established by the closed queueing network for
reliability analysis of the DBSS with unusable bikes. Based on this, together with the virtual closed queueing network, we can
provide a detailed analysis for four different research aspects as follows:

\textit{A nonlinear routing matrix equation:} In the theory of closed
queueing networks, it is well known that the relative arrival rates can
uniquely be determined by a system of linear equations $\mathbf{e}P=\mathbf{
}$ and $\mathbf{e}_{1}=1$, where $\mathbf{e}$ is the relative arrival rate
vector, $\mathbf{e}_{1}$ is the first element of the vector $\mathbf{e}$,
and $P$ is the routing matrix, e.g., see Bolch et al. \cite{Bol:2006} and
Serfozo \cite{Ser:1999} for more details. However, when studying the DBSS
with unusable bikes, this paper is the first to find a new fundamental
result: The relative arrival rate vector must be determined by a nonlinear
routing matrix equation $\mathbf{e}P\left( \mathbf{e}\right) =\mathbf{e}$,
where some entries of the routing matrix $P\left( \mathbf{e}\right) $ depend
on the relative arrival rate vector so that the routing matrix equation
\mathbf{e}P\left( \mathbf{e}\right) =\mathbf{e}$ is nonlinear in $\mathbf{e}
. Therefore, the nonlinearity of the routing matrix equation opens a new
research avenue in the study of closed queueing networks, and specially, it
will be useful in the analysis of the DBSSs with unusable bikes.

\textit{A new class of virtual closed queueing networks:} It is an effective
method that the DBSS is described as a virtual closed queueing network,
including virtual customers, virtual servers and virtual routines. Such a
closed queueing network described in this paper is close to that in our
previous three papers by Li et al. \cite{LiF:2016,Li:2017a,Li:2017b}. It is
worthwhile to note that the bikes' failure, removing, repair, distributing
and reuse processes, together with two batch policies for removing the
unusable bikes and redistributing the repaired bikes, can substantially
change the physical structure of the virtual closed queueing network, and
hence this work is very different from those previous three papers. Although
Li et al. \cite{Li:2017a} introduced several more general random factors
(e.g., the Markovian arrival processes of customers at stations, and the
phase-type time of ridding a bike on every road), the physical structure of its
corresponding closed queueing network is not changed so that the routing
matrix is still independent of the relative arrival rates; that is, the relative
arrival rates can be determined by the system of linear equations $\mathbf{e
P=\mathbf{e}$ and $\mathbf{e}_{1}=1$. Thus this paper is the first to find such a
nonlinear routing matrix equation in the study of closed queueing networks.

\textit{The block-structured Markov process: }This paper also contributes to
the literature of closed queueing networks through extending and
generalizing a simple queueing system (e.g., the M/M/1 queue, the M/M/C
queue, the M/G/1 queue and others) of a node to a more general
block-structured Markov process. We show that such a block structure
is established by either the bike's failure and batch removing processes in
a parking region or the bike's repair and batch redistribution processes in
the maintenance shop. Moreover, for the closed queueing network corresponding
to the DBSS with unusable bikes, once the relative arrival rate is
introduced to each node, all the nodes are isolated from each other so that
the Markov system in either each parking region or the maintenance shop can
be described as a block-structured Markov process. Based on this, this paper
gives a key result: The stationary probabilities of the DBSS with unusable
bikes are proved to have the product-form solution of the closed queueing
network, and they can be expressed under a unified framework of
block-structured Markov processes. Furthermore, the paper can provide
performance analysis of the DBSS with unusable bikes through applying the
block-structured Markov processes which are studied by using the
RG-factorizations.

\textit{A novel large-scale maintenance network: }This paper shows that the
DBSS with unusable bikes is a large-scale maintenance network, in which the
bikes may fail at every parking region, the unusable bikes are batch removed
from the parking regions to the maintenance shop, and the repaired bikes are
batch redistributed for their reuse from the maintenance shop to the parking
regions. However, so far no paper has analyzed such a maintenance network
with a lot of unusable bikes. It is worth noting that a large literature
focuses on manufacturing reliability engineering, e.g., see Buzacott and
Shanthikumar \cite{Buz:1993}, Gershwin \cite{Ger:1994} and Li and Meerkov
\cite{LiM:2009}. Despite this, limited research methods can be applied to deal with
the DBSS with unusable bikes.

\vskip0.3cm

The main contributions of this paper are summarized as follows:

(1) To the best of our knowledge, among the studies of bike-sharing systems
with unusable bikes, this paper is the first to propose and develop a new
computational method by applying the RG-factorizations of block-structured
Markov processes in the closed queueing networks.

(2) The DBSS with unusable bikes, together with two batch policies (one for
removing the unusable bikes, and the other for redistributing the repaired
bikes), is described as a new virtual closed queueing network. Also, such a
virtual closed queueing network, together with the bikes' failure, removing,
repair, redistribution and reuse processes, is a necessary and useful
extension and generalization of our previous works by Li et al. \cit
{LiF:2016,Li:2017a,Li:2017b}.

In addition, we believe that the generalized virtual closed queueing network
will also be useful in the study of other stochastic systems in practice.

(3) This paper extends the theory of closed queueing networks from two
fundamental aspects: First, this paper is the first to set up a new
nonlinear matrix equation (i.e., $\mathbf{e}P\left( \mathbf{e}\right)
\mathbf{e}$) to determine the relative arrival rates, where the nonlinearity
is caused by two classes of processes: The bikes' failure and batch
removing processes in various parking regions, and the repair and batch
redistribution processes in the maintenance shop. Second, we extend and
generalize simple queueing systems of some nodes (e.g., the \textit{parking
region} nodes and the \textit{maintenance shop} node) to more general
block-structured Markov processes, whose stationary probabilities can easily
be computed by means of the RG-factorizations. Based on this, we can further
extend and generalize the product-form solution of the closed queueing
networks to a more general case that the Markov systems of the nodes are
described as the block-structured Markov processes.

(4) This paper is the first to be able to provide performance analysis of DBSSs with
unusable bikes by means of our new computational method through applying the
RG-factorizations of block-structured Markov processes in the closed
queueing networks, and some numerical experiments are also provided to
analyze performance measures of the DBSS. Therefore, this paper can provide
insights for design and operations management of the DBSSs with unusable
bikes. We hope that our computational method can open a new avenue to
quantitative evaluation of more general DBSSs with unusable bikes.

\vskip0.3cm

The remainder of this paper is organized as follows. Section 2 provides
literature review. Section 3 provides model description for a DBSS with
unusable bikes. Section 4 abstracts the DBSS as a closed queueing network,
and provides a detailed analysis for the relative arrival rates, the service
processes, the block-structured Markov processes, and the nonlinear routing
matrix equation. Section 5 provides a more general product-form solution of
the closed queueing network and gives performance analysis of the DBSS with
unusable bikes. Section 6 uses numerical experiments to verify that our
computational method is effective for performance analysis of the DBSS with unusable bikes.
Concluding remarks are given in Section 7. Finally, four necessary
appendices are provided.

\section{Literature Review}

Our current work is related to four streams of literature: DBSSs with
unusable bikes; virtual closed queueing networks of bike-sharing systems;
maintenance networks and reliability analysis; and the RG-factorizations. We
review some related literature from the four different aspects.

\textbf{DBSSs with unusable bikes:} DBSS are the fourth generation of
bike-sharing systems due to applications of some new information and network
technologies, for example, GPS, GIS, platform service, smartphone, mobile
payment and so forth. The practical use of DBSSs started from 2014, but
little research on DBSSs has been available up to now. Reiss and Bogenberger
\cite{Rei:2016} built a demand model of DBSSs to forecast the upcoming
demand of bikes at certain time and place, and their model used some
historical booking data for a few months. Pal et al. \cite{Pal:2017}
established a novel mixed integer linear program to solve the static
repositioning problem of DBSSs, and their algorithm can deal with the
increase in scale of the static repositioning problem. Pal et al. \cit
{Pal:2018} took interactions among independent variables into consideration,
and showed that more insights about the mobility patterns and imbalance can
be obtained by considering such interactions. Liu et al. \cite{Liu:2018}
studied a bike repositioning problem with multiple depots, multiple visits,
and multiple heterogeneous vehicles for a DBSS. They proposed an enhanced
version of chemical reaction optimization (CRO), and verified that the
enhanced CRO had the potential to tackle the repositioning problem for a
longer repositioning duration and more vehicle instances. Caggiani et al.
\cite{Cag:2018} proposed a comprehensive, dynamic, and operator-based bike
repositioning method. They first investigated the spatio-temporal
correlation pattern of different zones in a city by adopting clustering
techniques on the historical usage records of bike sharing demand. Then they
applied a nonlinear autoregressive neural network demand forecasting model
to forecast the number of bikes needed. Finally, they proposed a decision
support system for dynamic repositioning to determine the optimal
repositioning flows and distribution patterns.

Different from the above works in the literature, this paper focuses on
performance analysis of DBSSs with unusable bikes by means of the closed
queueing networks and the block-structured Markov processes, and develops a
new computational method by applying RG-factorizations of block-structured
Markov process in the closed queueing networks.

On the other hand, in the study of DBSSs, unusable bikes are always a key negative factor that
affects not only sustainable development of bike-sharing systems but also
the quality of service in the DBSS. To keep sustainable development, it is
necessary to set up a useful control mechanism which enables fast removal of
lots of unusable bikes from every parking region to the maintenance shop.
Recently, some authors discussed bike-sharing systems with unusable bikes in
both station-based bike-sharing systems (SBSSs) and DBSSs.

(a) \textit{Unusable bikes in the SBSSs}: Kaspi et al. \cite{Kas:2016}
estimated the probability that each bike was usable, and formulated a
Bayesian model making use of on-line transactions data to constantly update
these probabilities. Kaspi et al. \cite{Kas:2017} measured user
dissatisfaction by using a weighted sum of the expected shortages of bikes
and lockers at one station, and the expected shortages can be expressed as a
utility function of initial inventory levels of usable bikes and unusable
bikes at the station.

(b) \textit{Unusable bikes in the DBSSs}: Wang and Szeto \cite{Wan:2019}
studied the static repositioning of both usable bikes and unusable bikes in
a DBSS. They applied a mixed integer linear program to balance the demand
and supply of bikes at each parking region. Usama et al. \cite{Usa:2019}
analyzed a repositioning policy of unusable bikes that are brought to the
maintenance shops by trucks.

\textbf{Closed queueing networks:} Closed queueing networks can be
applied to the study of bike-sharing systems (SBSSs and DBSSs). Li et al.
\cite{LiF:2016} is the first to construct a virtual closed queueing network
to describe and analyze a general SBSS with exponential or Poisson factors,
in which the bikes are regarded as ``virtual customers''; while the parking
regions and roads are viewed as ``virtual nodes''. For such a virtual closed
queueing network, Li et al. \cite{Li:2017a} further extended and generalized
it from three aspects: The road structure is an irreducible graph; the real
customer arrivals at the stations are Markovian arrival processes (MAPs);
and the ridding-bike times on the roads are i.i.d. with phase-type (PH)
distributions. Further, Li et al. \cite{Li:2017b} applied the fluid and diffusion
limits to analyze the SBSS by means of the virtual closed queueing network.
Along such a research line, this paper is an interesting extension and
generalization of the virtual closed queueing network both from the nonlinear routing matrix equation and from the block-structured Markov processes.

In addition, as a closely related application of closed queueing networks,
the rental systems have been drawn high research attention. Important
examples include Adelman \cite{Ade:2007}, George and Xia \cite{Geo:2011},
George \cite{Geo:2012}, Fanti et al. \cite{Fan:2014}, Samet et al. \cit
{Sam:2018} and so forth.

\textbf{Maintenance networks and reliability analysis: }This paper shows
that the DBSS with unusable bikes is a large-scale maintenance network, in
which bikes may fail at every parking region, and the unusable bikes are
collected from the parking regions to a maintenance shop for their repair
and reuse purpose. This paper is the first to discuss the DBSSs with
unusable bikes from the research perspective of maintenance networks and
reliability analysis.

So far there has been much research on manufacturing engineering with
maintenance, reliability and safety, such as manufacturing systems with
repairable machines by Buzacott and Yao \cite{Buz:1986}, Buzacott and Shanthikumar \cite{Buz:1993}
and Gershwin \cite{Ger:1994}, production systems with repairable stations by
Li and Meerkov \cite{LiM:2009}, and parts inventory systems with repairable
parts by Park and Lee \cite{Par:2011, Par:2014}. Also, much
recent literature indicates that queueing networks are an effective
mathematical tool to deal with reliability problems of many manufacturing
systems with maintenance, reliability and safety issues.

Compared with the literature of manufacturing systems with maintenance,
reliability and safety issues, the maintenance process of DBSSs with
unusable bikes is more complicated and challenging due to the following
factors: (1) Many bikes in different parking regions may simultaneously
fail, and the failed bikes need to be removed from each parking region to
the maintenance shop. Also, the batch removing policies of unusable bikes play a
key role in the study of DBSSs with unusable bikes; (2) The unusable
bikes are collected by trucks in batch from each parking region to the
maintenance shop, and the repaired bikes in the maintenance shop also need
to be redistributed in batch by trucks to various parking regions.
Therefore, such a DBSS maintenance network with a lot of unusable bikes and
their repair and distribution, together with two batch policies, is more interesting and
challenging in the study of stochastic maintenance networks.

\textbf{The RG-factorizations:} For block-structured Markov processes, Li
\cite{Li:2010} provided a unified effective computational framework of the
RG-factorizations, which includes stationary performance analysis, transient
solution, the first passage times, reward processes, quasi-stationary
distribution and so forth. Some practical applications of the
RG-factorizations include, for example, the first passage time and the
sojourn time by Li and Cao \cite{LiC:2004}, the semiconductor production
line by Li et al. \cite{LiZ:2009}, and computer and communication networks by
Wang et al. \cite{Wan:2007, Wan:2010}.

For the DBSS with unusable bikes, this paper provides a more general
product-form solution of the closed queueing network, and gives performance
analysis of the DBSS with unusable bikes. Therefore, this paper can
successfully analyze a large-scale DBSS with unusable bikes through applying
the RG-factorizations of block-structured Markov processes in the closed
queueing networks.

\section{Model Description}

In this section, we describe a DBSS with $N$ parking regions and one
maintenance shop. The routing among them is modeled as an irreducible path
graph, and two batch policies are introduced to remove the unusable bikes
and to redistribute the repaired bikes.

We provide a more detailed description for the DBSS with unusable bikes,
including the operation mechanism, system parameters, and mathematical notation
as follows:

\textbf{(1) Parking regions:} We assume that there are $N$ parking regions
in the DBSS. The $N$ parking regions may be different due to their
geographical locations and surrounding traffic environment, such as the subway
station, commercial center, residential area, uphill road, downhill road and
so forth. We assume that the capacity of each parking region is sufficiently
large, and the total number of bikes in the DBSS is fixed at $K$.

\textbf{(2) Roads:} The parking regions and the maintenance shop are
referred to as nodes. There are roads between any two nodes. A user arrives
at a parking region, rents a usable bike if any and rides it for a while.
Once he finishes his trip, he returns the bike to a parking region and immediately
leaves the system. There might be unusable bikes at some parking regions. We
need to use trucks to remove unusable bikes from some parking regions to the
maintenance shop, and redistribute repaired bikes from the maintenance shop
to some suitable parking regions. Thus it is necessary to classify the roads as bike
roads and truck roads as follows:

\textbf{Bike roads:} Let Road $i\rightarrow j$ denote a road for riding
bikes from parking region $i$ to parking region $j$. Note that Road
i\rightarrow j$ and Road $j\rightarrow i$ may be different due to different
practical factors, for example, an uphill road versus a downhill road.

Denote all the bike-riding roads beginning from parking region $i$ for
1\leq i\leq N$ a
\begin{equation*}
R_{B}\left( i\right) =\left\{ \text{Road }i\rightarrow j:j\neq i,1\leq j\leq
N\right\} .
\end{equation*
To express all the parking regions in the downlink of parking region $i$, we
writ
\begin{equation*}
\Theta _{i}=\left\{ \text{parking region }j:\text{Road }i\rightarrow j\in
R_{B}\left( i\right) \right\} .
\end{equation*}

\textbf{Truck roads: }We divide the truck roads into two classes: (i) The
truck roads for removing unusable bikes from a parking region to the
maintenance shop; (ii) The truck roads for redistributing repaired bikes
from the maintenance shop to the parking regions. We denote the maintenance
shop by Node 0. Then the set of truck roads of class one is denoted a
\begin{equation*}
R_{E}\left( 0\right) =\left\{ \text{Road }i\rightarrow 0:1\leq i\leq
N\right\} ,
\end{equation*
and the set of truck roads of class two as
\begin{equation*}
R_{B}\left( 0\right) =\left\{ \text{Road }0\rightarrow i:1\leq i\leq
N\right\} .
\end{equation*}

\textbf{(3) User arrival processes:} The arrivals of outside users at
parking region $i$\ follow a Poisson process with arrival rate $\lambda
_{i}>0$ for $1\leq i\leq N$.

\textbf{(4) The bike-renting and -riding processes: }An outside user arrives
at parking region $i$ in order to rent a bike. If there is no usable bike
(either the parking region $i$ is empty or all the bikes in parking region
i $ are unusable), then the user immediately leaves the system. If there are
at least one usable bike in parking region $i$, then the user rents a bike
and rides it on Road $i\rightarrow j$ with
probability $p_{i,j}$ for $j\neq i$, $1\leq j\leq N$ and $\sum\nolimits_{j\in \Theta _{i}}p_{i,j}=1$ for
each $i=1,2,\ldots ,N$. We assume that the bike-riding times on Road
i\rightarrow j$ are i.i.d. and exponential with rate $\mu _{i,j}>0$.
\begin{figure}[h]
\centering         \includegraphics[width=8.5cm]{fig-1.eps} \newline
\caption{Key factors of the bike-sharing system}
\label{figure:figure-1}
\end{figure}

\textbf{(5) The lifetime of bikes:} We assume that all the bikes are
identical, and the lifetime of a bike is exponential with failure rate
\alpha >0$.

\textbf{(6) A maintenance shop:} There is a maintenance shop whose purpose
is to repair all the unusable bikes, which are collected from every parking
region. For simplicity of description, we assume that a repaired bike is the
same as an initial new bike. Furthermore, if a seriously damaged bike is
scrapped after a repair, then a new bike is supplemented immediately. Note
that such an assumption keeps the total number of bikes in the system
constant, so that we can apply the closed queueing network in our later
study. We also assume that the repair time of each unusable bike is
exponential with rate $w>0$, and there are $r$ repairmen in the maintenance
shop.

\textbf{(7) A batch removing policy:} Once the number of unusable bikes at
any parking region reaches a positive integer $M$, a truck immediately
removes the $M$ unusable bikes in a batch and transports the unusable bikes
to the maintenance shop. We assume the transportation time by truck on Road
i\rightarrow 0$ is exponential with rate $\mu _{i,0}>0$.
\begin{figure}[h]
\centering                         \includegraphics[width=9cm]{fig-2.eps}
\newline
\caption{The physical structure of the DBSS}
\label{figure:figure-2}
\end{figure}

\textbf{(8) A batch redistribution policy: }Once the number of repaired
bikes in the maintenance shop reaches a given positive integer $Z$, a truck
immediately takes the $Z$ repaired bikes in a batch away from the
maintenance shop to some parking regions, in which $Z_{i}$ repaired bikes
are sent to the parking region $i$ with $\sum_{i=1}^{N}Z_{i}=Z$. Let $\beta
_{i}=Z_{i}/Z$. The transportation time by truck from the maintenance shop to
parking region $i$ is exponential with rate $\mu _{0,i}>0$.

For convenience of expression, we assume that $Z=\psi M$, where $\psi $ is a
given positive integer. In addition, let $\phi =\lfloor K/M\rfloor $.

\textbf{(9) The departure discipline:} The departure process of users has
two different cases: (a) A user directly leaves the DBSS if he arrives at a
parking region without any usable bike. (b) When a user completes his trip,
he returns the bike to a parking region, and leaves the DBSS immediately.

We assume that all the random variables mentioned above are independent of
each other.

The key elements of the DBSS, including parking regions, unusable bikes,
trucks for transporting bikes, and a maintenance shop, are demonstrated in
Figure 1. Figure 2 illustrates the physical structure of the DBSS with
unusable bikes.

\begin{Rem}
In the above (4), the probability $p_{i,j}$ that a bike in parking region $i$
is ridden to parking region $j\in \Theta _{i},1\leq i\leq N$, can be
statistically determined from practical data in the DBSS. Note that data
analysis is based on effective statistical methods and simulation
techniques. Similarly, we can also derive $\beta _{i}$ in the above (8).
\end{Rem}

The notation used in this paper is summarized as follows:

\begin{array}{ll}
N\ \ \ \ \ \ \  & \text{Number of parking regions;} \\
K & \text{Total number of bikes;} \\
\Theta _{i} & \text{Set of parking regions in the downlink of parking region
}i\text{;} \\
\lambda _{i}\  & \text{User arrival rate at parking region }i\text{;} \\
\mu _{i,j} & \text{Riding rate of bike-riding on Road }i\rightarrow j\text{;}
\\
p_{i,j}\ \ \ \ \ \ \  & \text{Probability of entering Road }i\rightarrow
\text{ from parking region }i\text{;} \\
\alpha & \text{Failure rate of a bike;} \\
w & \text{Repair rate of an unusable bike;} \\
r & \text{Number of repairmen in the maintenance shop;} \\
M & \text{Batch size of removing unusable bikes from any parking region;} \\
Z & \text{Batch size of redistributing repaired bikes from the maintenance
shop;} \\
\mu _{i,0}\  &
\begin{array}{l}
\text{Transportation rate by truck from parking region }i\text{ to the} \\
\text{maintenance shop;
\end{array}
\\
\mu _{0,i} &
\begin{array}{l}
\text{Transportation rate by truck from the maintenance shop} \\
\text{to parking region }i\text{;
\end{array}
\\
Z_{i} & \text{Number of repaired bikes redistributed to parking region }
\text{;} \\
\beta _{i}\ \ \ \ \ \ \ \ \ \ \ \ \  & Z_{i}/Z\text{;} \\
\phi & \lfloor K/M\rfloor \text{;} \\
\psi & Z/M\text{, a positive integer.
\end{array
$

\section{A Virtual Closed Queueing Network}

In this section, we describe the DBSS with unusable bikes as a virtual closed
queueing network, and explain its physical structure and mathematical
notation.

Since the total number of bikes (usable and unusable) in the DBSS is fixed,
we can abstract the DBSS as a closed queueing network whose nodes and
customers are defined as follows.

\textbf{(1) Virtual customers:} The virtual customers are the bikes.
Obviously, bikes are either parked in parking regions, ridden or transported
on roads, or repaired in the maintenance shop. There are two states of
bikes: Usable and unusable. Thus we use $G$ ($G$: Good) and $B$ ($B$: Bad)
to denote the usable and unusable bikes, respectively.

\textbf{(2) Virtual nodes:} Note that parking regions, roads, and the
maintenance shop have different physical attributes such as functions and
geographical structures. Thus they are considered as different classes of
nodes.

\textbf{(3) An irreducible path graph:} The set of all virtual nodes in the
DBSS is given b
\begin{equation*}
\Theta =\left\{ \text{Node }i:0\leq i\leq N\right\} \underset{i=1}{\overset{
}{\cup }}R_{B}\left( i\right) \cup R_{B}\left( 0\right) \cup R_{E}\left(
0\right) .
\end{equation*
We assume that all the roads of the DBSS are connected as an irreducible
path graph whose nodes are in the set $\Theta $. In this case, we call that
the DBSS is path irreducible.

Let $Q_{G}^{\left( i\right) }\left( t\right) $ and $Q_{B}^{\left( i\right)
}\left( t\right) $ denote the numbers of usable bikes and of unusable bikes
parked in parking region $i$ at time $t\geq 0$ for $1\leq i\leq N$,
respectively; and $Q_{G}^{\left( 0\right) }\left( t\right) $ and
Q_{B}^{\left( 0\right) }\left( t\right) $ the numbers of usable bikes and of
unusable bikes in the maintenance shop at time $t\geq 0$, respectively. Let
R_{i,j}\left( t\right) $ be the number of bikes ridden on Road $i\rightarrow
j$ at time $t$ for $1\leq i\leq N$ and $j\in \Theta _{i}$. $R_{i,0}\left(
t\right) $ and $R_{0,i}\left( t\right) $ are numbers of bikes on Road
i\rightarrow 0$ and Road $0\rightarrow i$\ at time $t$ for $1\leq i\leq N$,
respectively.

We writ
\begin{equation*}
X\left( t\right) =\left( \mathbf{L}_{0}\left( t\right) ,\mathbf{L}_{1}\left(
t\right) ,\mathbf{L}_{2}\left( t\right) ,\ldots ,\mathbf{L}_{N-1}\left(
t\right) ,\mathbf{L}_{N}\left( t\right) \right) ,
\end{equation*
\newline
wher
\begin{equation*}
\mathbf{L}_{0}\left( t\right) =\left( Q_{G}^{\left( 0\right) }\left(
t\right) ,Q_{B}^{\left( 0\right) }\left( t\right) ;R_{i,0}\left( t\right)
,R_{0,i}\left( t\right) :1\leq i\leq N\right) ,
\end{equation*
and for $1\leq i\leq N
\begin{equation*}
\mathbf{L}_{i}\left( t\right) =\left( Q_{G}^{\left( i\right) }\left(
t\right) ,Q_{B}^{\left( i\right) }\left( t\right) ,R_{i,j}\left( t\right)
:j\in \Theta _{i}\right) .\text{ }
\end{equation*
Obviously, $\left\{ X\left( t\right) :t\geq 0\right\} $ is a continuous-time
irreducible Markov process. The state space of the Markov process $\left\{
X\left( t\right) :t\geq 0\right\} $ is given b
\begin{eqnarray*}
\Omega &=&\left\{ \mathbf{n}:\sum_{i=0}^{N}\left( n_{G}^{\left( i\right)
}+n_{B}^{\left( i\right) }\right) +\sum_{i=1}^{N}\sum_{j\in \Theta
_{i}}m_{i,j}+\sum_{i=1}^{N}\left( m_{i,0}+m_{0,i}\right) =K,\right. \\
&&\left. 0\leq n_{B}^{\left( i\right) }\leq M,0\leq n_{G}^{\left( i\right)
},m_{i,j}\leq K,0\leq n_{G}^{\left( 0\right) }\leq Z,0\leq n_{B}^{\left(
0\right) }\leq \phi M,\right. \\
&&\left. n_{G}^{\left( 0\right) }+n_{B}^{\left( 0\right)
}=kM,m_{i,0}=kM,m_{0,i}=lZ_{i},\text{for }1\leq i\leq N,j\in \Theta
_{i},\right. \\
&&\left. 0\leq k\leq \phi ,0\leq l\leq \phi /\psi ,0\leq Z_{i}\leq Z\right\}
,
\end{eqnarray*
wher
\begin{equation*}
\mathbf{n}=\left( \mathbf{n}_{0},\mathbf{n}_{1},\mathbf{n}_{2},\ldots
\mathbf{n}_{N-1},\mathbf{n}_{N}\right) ,
\end{equation*
an
\begin{equation*}
\mathbf{n}_{0}=\left( n_{G}^{\left( 0\right) },n_{B}^{\left( 0\right)
};m_{i,0},m_{0,i}:1\leq i\leq N\right) ,
\end{equation*
for $1\leq i\leq N$
\begin{equation*}
\mathbf{n}_{i}=\left( n_{G}^{\left( i\right) },n_{B}^{\left( i\right)
};m_{i,j}:j\in \Theta _{i}\right) .
\end{equation*
Note that $n_{G}^{\left( i\right) }$ and $n_{B}^{\left( i\right) }$\ are the
numbers of usable bikes and of unusable bikes parked in parking region $i$
for $1\leq i\leq N$, respectively; $n_{G}^{\left( 0\right) }$ and
n_{B}^{\left( 0\right) }$ are the numbers of usable bikes and of unusable
bikes in the maintenance shop, respectively; $m_{i,j}$ is the number of
bikes ridden on Road $i\rightarrow j$ for $1\leq i\leq N$, $j\in \Theta _{i}
; and $m_{i,0},m_{0,i}$ denote the numbers of unusable bikes and of repaired
bikes transported on Road $i\rightarrow 0$ and Road $0\rightarrow i$ for
1\leq i\leq N$, respectively.

Figure 3 illustrates the physical structure of the closed queueing network.
To study the closed queueing network, in the remainder of this section we
will analyze the relative arrival rates, the Markov processes at the nodes,
the service rates, the routing matrix, and computation of the relative
arrival rates.
\begin{figure}[h]
\centering                         \includegraphics[width=8.5cm]{fig-3.eps}
\newline
\caption{The structure of the virtual closed queueing network}
\label{figure:figure-3}
\end{figure}

\subsection{The relative arrival rates}

Note that the relative arrival rates play a key role in the derivation of
the product-form solution of queueing networks, e.g., see Bolch et al. \cit
{Bol:2006} and Serfozo \cite{Ser:1999} for more details. Here, we describe
how to introduce the relative arrival rates to each node.

We denote by $e_{i}$ the relative arrival rate of Node $i$ for $0\leq i\leq
N $, $e_{i,j}$ the relative arrival rate of Node $i\rightarrow j$ for $1\leq
i\leq N$ and $j\in \Theta _{i}$, $e_{i,0}$ and $e_{0,i}$ the relative
arrival rates of Node $i\rightarrow 0$ and Node $0\rightarrow i$,
respectively. Le
\begin{equation*}
\mathbb{E=}\left\{ \mathbf{e}_{0},\mathbf{e}_{1},\mathbf{e}_{2},\ldots
\mathbf{e}_{N-1},\mathbf{e}_{N}\right\} ,
\end{equation*
wher
\begin{equation*}
\mathbf{e}_{0}=\left( e_{0};e_{i,0},e_{0,i}:1\leq i\leq N\right) ,
\end{equation*
and for $1\leq i\leq N
\begin{equation*}
\mathbf{e}_{i}=\left( e_{i};e_{i,j}:\text{ }j\in \Theta _{i}\right) .
\end{equation*}

Note that the relative arrival rates only depend on the physical structure
of the closed queueing network, while they are independent of the states of
the Markov process $\left\{ X\left( t\right) :t\geq 0\right\} $.

To determine the relative arrival rates, we first need to analyze the Markov
processes of various parking regions and of the maintenance shop in
Subsections 4.2 and 4.3, respectively. Then we need to establish the
routing matrix in Subsection 4.4, so that we can set up a nonlinear routing
matrix equation to determine the relative arrival rates.

\subsection{Markov process of any parking region}

When the relative arrival rates are introduced to various nodes in the
closed queueing network, all the nodes are isolated from each other so that
the Markov processes of the nodes are independent of each other. Therefore,
the Markov system of each node is a block-structured Markov process, which
is related to the numbers of usable bikes and of unusable bikes in either each
parking region or the maintenance shop.

To set up the block-structured Markov process, it is necessary to analyze
the state transition relations of the Markov process in any parking region.
Note that $Q_{G}^{\left( i\right) }\left( t\right) $ and $Q_{B}^{\left(
i\right) }\left( t\right) $ are the numbers of usable bikes and of unusable
bikes at parking region $i$ at time $t$, respectively. Under the batch
removing policy of unusable bikes with parameter $M$, it is easy to see that
$\mathbf{Q}=\left\{ \left( Q_{G}^{\left( i\right) }\left( t\right)
,Q_{B}^{\left( i\right) }\left( t\right) \right) :t\geq 0\right\} $ is a
two-dimensional Markov process with finite levels and phases, where
Q_{G}^{\left( i\right) }\left( t\right) $ is the phase variable while
Q_{B}^{\left( i\right) }\left( t\right) $ is the level variable. Figure 4
depicts the state transition relations of the Markov process $\mathbf{Q}$.
Thus the infinitesimal generator $\mathbf{Q}$ of the Markov process is given
by

\begin{figure}[tbp]
\centering                         \includegraphics[width=13cm]{fig-4.eps}
\newline
\caption{State transition relations of the Markov process in parking region
i$ }
\label{figure:figure-4}
\end{figure}

\begin{equation}
\mathbf{Q}=\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{M-1,M-1} & Q_{M-1,M} \\
Q_{M,0} &  &  &  & Q_{M,M
\end{array
\right) .  \label{V-1}
\end{equation
Note that all the elements of the infinitesimal generator $\mathbf{Q}$ are
given in Appendix A.

The structure of the matrix $\mathbf{Q}$ is beautiful with bidiagonal blocks
and a special block $Q_{M,0}$ in the lower-left corner. Obviously, the
matrix-geometric method (see Neuts \cite{Neu:1981} or Latouche and Ramaswami
\cite{Lat:1999}) will be inapplicable for such a Markov process $\mathbf{Q}
. Thus, to compute the stationary probabilities of the Markov process
\mathbf{Q}$, we need to use the UL-type RG-factorization (see Li \cit
{Li:2010} for more details). To this end, we writ
\begin{equation*}
Q_{M,M}^{[\leq M]}=Q_{M,M},
\end{equation*
and for $1\leq k\leq M-1$,
\begin{equation*}
Q^{[\leq k]}=\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{k-1,k-1} & Q_{k-1,k} \\
Q_{k,0}^{[\leq k]} &  &  &  & Q_{k,k
\end{array
\right) .
\end{equation*}

\begin{Lem}
We hav
\begin{equation*}
Q_{M,M}^{[\leq M]}=Q_{M,M},
\end{equation*
for $1\leq k\leq M-1,
\begin{equation}
Q_{k,0}^{[\leq k]}=\prod\nolimits_{l=k}^{M-1}Q_{l,l+1}\left(
-Q_{l+1,l+1}\right) ^{-1}Q_{M,0},  \label{eq-Qk}
\end{equation
and
\begin{equation}
Q^{\left[ \leq 0\right] }=Q_{0,0}+Q_{0,1}\left( -Q_{1,1}\right) ^{-1}\left[
\prod\limits_{l=1}^{M-1}Q_{l,l+1}\left( -Q_{l+1,l+1}\right) ^{-1}\right]
Q_{M,0}.  \label{eq-Q0}
\end{equation}
\end{Lem}

\textbf{Proof: }See Appendix B.

Le
\begin{equation}
\Psi _{n}=Q_{n,n}^{[\leq n]},\text{ }0\leq n\leq M,  \label{eq-RG}
\end{equation
\begin{equation}
R_{i,j}=Q_{i,j}^{\left[ \leq j\right] }\left( -\Psi _{j}\right) ^{-1},\text{
}0\leq i<j\leq M,  \label{eq-R}
\end{equation
an
\begin{equation}
G_{i,j}=\left( -\Psi _{i}\right) ^{-1}Q_{i,j}^{\left[ \leq i\right] },\text{
}0\leq j<i\leq M.  \label{eq-G}
\end{equation}

The following theorem gives the UL-type RG-factorization of the Markov
process $\mathbf{Q}$. It is easy to see that the RG-factorization has a
beautiful block structure.

\begin{The}
The UL-type RG-factorization of the continuous-time Markov process $\mathbf{
}$ is given b
\begin{equation*}
\mathbf{Q}=\left( I-R_{U}\right) \Psi _{D}\left( I-G_{L}\right) ,
\end{equation*
wher
\begin{equation*}
R_{U}=\left(
\begin{array}{ccccccc}
0 & R_{0,1} &  &  &  &  &  \\
& 0 & R_{1,2} &  &  &  &  \\
&  & 0 & R_{2,3} &  &  &  \\
&  &  & \ddots & \ddots &  &  \\
&  &  &  & 0 & R_{M-2,M-1} &  \\
&  &  &  &  & 0 & R_{M-1,M} \\
&  &  &  &  &  &
\end{array
\right) ,
\end{equation*
\begin{equation*}
\Psi _{D}=\text{diag}\left( \Psi _{0},\Psi _{1},\Psi _{2},\ldots ,\Psi
_{M-1},\Psi _{M}\right) ,
\end{equation*
an
\begin{equation*}
G_{L}=\left(
\begin{array}{ccccccc}
0 &  &  &  &  &  &  \\
G_{1,0} & 0 &  &  &  &  &  \\
G_{2,0} &  & 0 &  &  &  &  \\
\vdots &  &  & \ddots &  &  &  \\
G_{M-2,0} &  &  &  & 0 &  &  \\
G_{M-1,0} &  &  &  &  & 0 &  \\
G_{M,0} & 0 & 0 & \cdots & 0 & 0 &
\end{array
\right) .
\end{equation*}
\end{The}

\textbf{Proof:} See Appendix C.

Since the Markov process $\mathbf{Q}$ is irreducible, finite states, and
\mathbf{Q}\mathbbm{1}=0$ where $\mathbbm{1}$ is a column vector of ones, it
must be positive recurrent. Let $\mathbf{\pi }_{B}^{\left( i\right) }=\left(
\mathbf{\pi }_{B,0}^{\left( i\right) },\mathbf{\pi }_{B,1}^{\left( i\right)
},\mathbf{\pi }_{B,2}^{\left( i\right) },\ldots ,\mathbf{\pi
_{B,M-1}^{\left( i\right) },\mathbf{\pi }_{B,M}^{\left( i\right) }\right) $
be the stationary probability vector of the Markov process $\mathbf{Q}$,
wher
\begin{equation*}
\left\{
\begin{array}{l}
\mathbf{\pi }_{B,0}^{\left( i\right) }=\left( \pi _{G,0;B,0}^{\left(
i\right) },\pi _{G,1;B,0}^{\left( i\right) },\pi _{G,2;B,0}^{\left( i\right)
},\ldots ,\pi _{G,K-1;B,0}^{\left( i\right) },\pi _{G,K;B,0}^{\left(
i\right) }\right) , \\
\mathbf{\pi }_{B,1}^{\left( i\right) }=\left( \pi _{G,0;B,1}^{\left(
i\right) },\pi _{G,1;B,1}^{\left( i\right) },\pi _{G,2;B,1}^{\left( i\right)
},\ldots ,\pi _{G,K-2;B,1}^{\left( i\right) },\pi _{G,K-1;B,1}^{\left(
i\right) }\right) , \\
\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\vdots \\
\mathbf{\pi }_{B,M}^{\left( i\right) }=\left( \pi _{G,0;B,M}^{\left(
i\right) },\pi _{G,1;B,M}^{\left( i\right) },\pi _{G,2;B,M}^{\left( i\right)
},\ldots ,\pi _{G,K-M-1;B,M}^{\left( i\right) },\pi _{G,K-M;B,M}^{\left(
i\right) }\right)
\end{array
\right.
\end{equation*
Then using the UL-type RG-factorization, we obtai
\begin{equation}
\left\{
\begin{array}{l}
\mathbf{\pi }_{B,0}^{\left( i\right) }=\kappa \mathbf{x}_{0}, \\
\mathbf{\pi }_{B,k}^{\left( i\right) }=\mathbf{\pi }_{B,k-1}^{\left(
i\right) }R_{k-1,k},\text{ }1\leq k\leq M
\end{array
\right.  \label{eq-iM}
\end{equation
where $\mathbf{x}_{0}$ is the stationary probability vector of the censored
Markov chain $\Psi _{0}$ to level $0$, the scalar $\kappa $ is uniquely
determined by $\sum\nolimits_{k=0}^{M}\mathbf{\pi }_{B,k}^{\left( i\right)
\mathbbm{1}=1$, and $\mathbbm{1}$ is a column vector of ones. Note that the
censored Markov process $Q^{\left[ \leq 0\right] }$ has the stationary
probability vector $\mathbf{x}_{0}$, and thus we hav
\begin{equation}
\left\{
\begin{array}{l}
\mathbf{x}_{0}Q^{\left[ \leq 0\right] }=0, \\
\mathbf{x}_{0}\mathbbm{1}=1
\end{array
\right.  \label{eq-Q0x}
\end{equation}

\subsection{Markov process of the maintenance shop}

The repair behavior of the maintenance shop can be analyzed by a
two-dimensional continuous-time Markov process $\left\{ \left( Q_{B}^{\left(
0\right) }\left( t\right) ,Q_{G}^{\left( 0\right) }\left( t\right) \right)
t\geq 0\right\} $, where $Q_{B}^{\left( 0\right) }\left( t\right) $ is the
phase variable and $Q_{G}^{\left( 0\right) }\left( t\right) $ is the level
variable. When the number of repaired bikes amounts to $Z$, the $Z$ repaired
bikes are taken away by truck in batch from the maintenance shop to the
parking regions. We assume that the proportion of repaired bikes
redistributed to parking region $i$ is $\beta _{i}$, where $\beta
_{i}=Z_{i}/Z$ for $1\leq i\leq N$ and the transportation time of the $Z_{i}$
repaired bikes is exponential with rate $\mu _{0,i}$. We write $\mu
_{0}=\sum\nolimits_{i=1}^{N}\beta _{i}\mu _{0,i}$. Note that there are $r$
repairmen at the maintenance shop and the repair time of each unusable bike
is exponential with repair rate $w$. The repair rate of the maintenance shop
is given b
\begin{equation*}
\overline{w}\left( n_{B}^{\left( 0\right) }\right) =\min \left\{
n_{B}^{\left( 0\right) },r\right\} w,
\end{equation*
where $n_{B}^{\left( 0\right) }$ is the number of unusable bikes in the
maintenance shop.
\begin{figure}[h]
\centering                         \includegraphics[width=15cm]{fig-5.eps}
\newline
\caption{State transition relations of the Markov process in the maintenance
shop}
\label{figure:figure-5}
\end{figure}

Figure 5 depicts the state transition relations of the Markov process
\left\{ \left( Q_{B}^{\left( 0\right) }\left( t\right) ,Q_{G}^{\left(
0\right) }\left( t\right) \right) :\right. \newline
\left. t\geq 0\right\} $ in the maintenance shop. Obviously, this Markov
process is irreducible and positive recurrent. The infinitesimal generator
\mathbf{T}$ of the Markov process $\left\{ \left( Q_{B}^{\left( 0\right)
}\left( t\right) ,Q_{G}^{\left( 0\right) }\left( t\right) \right) :t\geq
0\right\} $ is given by
\begin{equation}
\mathbf{T}=\left(
\begin{array}{cccccccc}
T_{0,0} & T_{0,1} &  &  &  &  &  &  \\
& T_{1,1} & T_{1,2} &  &  &  &  &  \\
&  & \ddots & \ddots &  &  &  &  \\
&  &  & T_{M,M} & T_{M,M+1} &  &  &  \\
&  &  &  & T_{M+1,M+1} & T_{M+1,M+2} &  &  \\
&  &  &  &  & \ddots & \ddots &  \\
&  &  &  &  &  & T_{Z-1,Z-1} & T_{Z-1,Z} \\
T_{Z,0} &  &  &  &  &  &  & T_{Z,Z
\end{array
\right) ,  \label{V-2}
\end{equation
whose block elements are given in Appendix D.

For the two-dimensional Markov process $\mathbf{T}$, we write
T_{Z,Z}^{[\leq Z]}=T_{Z,Z}$. For $1\leq k\leq Z-1$, we can iteratively ge
\begin{equation*}
T^{[\leq k]}=\left(
\begin{array}{ccccc}
T_{0,0} & T_{0,1} &  &  &  \\
& T_{1,1} & T_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & T_{k-1,k-1} & T_{k-1,k} \\
\Xi _{k} &  &  &  & T_{k,k
\end{array
\right) ,
\end{equation*
where $\Xi _{k}=\prod\nolimits_{l=k}^{Z-1}\left[ T_{l,l+1}\left(
-T_{l+1,l+1}\right) ^{-1}\right] T_{Z,0}$.\ Note that $T^{\left[ \leq
\right] }=T^{\left[ 0\right] }$, and we obtai
\begin{equation}
T^{\left[ \leq 0\right] }=T_{0,0}+T_{0,1}\left( -T_{1,1}\right) ^{-1}\left[
\prod\limits_{l=1}^{Z-1}T_{l,l+1}\left( -T_{l+1,l+1}\right) ^{-1}\right]
T_{Z,0}.  \label{eq-T0}
\end{equation}

Le
\begin{equation*}
U_{n}=T_{n,n}^{[\leq n]},\text{ }0\leq n\leq Z,
\end{equation*
\begin{equation*}
\mathbf{R}_{i,j}=T_{i,j}^{\left[ \leq j\right] }\left( -U_{j}\right) ^{-1}
\text{ }0\leq i<j\leq Z,
\end{equation*
an
\begin{equation*}
\mathbf{G}_{i,j}=\left( -U_{i}\right) ^{-1}T_{i,j}^{\left[ \leq i\right] }
\text{ }0\leq j<i\leq Z.
\end{equation*
Then the UL-type RG-factorization of the Markov process $\mathbf{T}$ is
given b
\begin{equation*}
\mathbf{T}=\left( I-\mathbf{R}_{U}\right) U_{D}\left( I-\mathbf{G
_{L}\right) ,
\end{equation*
wher
\begin{equation*}
\mathbf{R}_{U}=\left(
\begin{array}{ccccccc}
0 & \mathbf{R}_{0,1} &  &  &  &  &  \\
& 0 & \mathbf{R}_{1,2} &  &  &  &  \\
&  & 0 & \mathbf{R}_{2,3} &  &  &  \\
&  &  & \ddots & \ddots &  &  \\
&  &  &  & 0 & \mathbf{R}_{Z-2,Z-1} &  \\
&  &  &  &  & 0 & \mathbf{R}_{Z-1,Z} \\
&  &  &  &  &  &
\end{array
\right) ,
\end{equation*
\begin{equation*}
U_{D}=\text{diag}\left( U_{0},U_{1},\ldots ,U_{Z-1},U_{Z}\right) ,
\end{equation*
an
\begin{equation*}
\mathbf{G}_{L}=\left(
\begin{array}{cccccc}
0 &  &  &  &  &  \\
\mathbf{G}_{1,0} & 0 &  &  &  &  \\
\mathbf{G}_{2,0} &  & 0 &  &  &  \\
\vdots &  &  & \ddots &  &  \\
\mathbf{G}_{Z-1,0} &  &  &  & 0 &  \\
\mathbf{G}_{Z,0} & 0 & \cdots & 0 & 0 &
\end{array
\right) .
\end{equation*}

Since the Markov process $\mathbf{T}$ is irreducible and positive recurrent,
there exists the stationary probability vector. Let $\mathbf{\pi
_{G}^{\left( 0\right) }=\left( \pi _{G,0}^{\left( 0\right) },\pi
_{G,1}^{\left( 0\right) },\pi _{G,2}^{\left( 0\right) },\ldots ,\pi
_{G,Z-1}^{\left( 0\right) },\pi _{G,Z}^{\left( 0\right) }\right) $ be the
stationary probability vector of the Markov process $\mathbf{T}$, wher
\begin{equation*}
\left\{
\begin{array}{l}
\mathbf{\pi }_{G,k}^{\left( 0\right) }=\left( \pi _{B,0;G,k}^{\left(
0\right) },\pi _{B,M;G,k}^{\left( 0\right) },\ldots ,\pi _{B,\left( \phi
-1\right) M-k;G,k}^{\left( 0\right) },\pi _{B,\phi M-k;G,k}^{\left( 0\right)
}\right) ,\text{ for }k=lM,\text{ }0\leq l\leq \psi , \\
\mathbf{\pi }_{G,k}^{\left( 0\right) }=\left( \pi _{B,\left( l+1\right)
M-k;G,k}^{\left( 0\right) },\pi _{B,\left( l+2\right) M-k;G,k}^{\left(
0\right) },\ldots ,\pi _{B,\left( \phi -1\right) M-k;G,k}^{\left( 0\right)
},\pi _{B,\phi M-k;G,k}^{\left( 0\right) }\right) ,\text{ for }lM \\
\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+1\leq k<\left( l+1\right) M,0\leq l\leq \psi -1
\end{array
\right.
\end{equation*}
Then by using the UL-type RG-factorizations, we obtai
\begin{equation}
\left\{
\begin{array}{l}
\mathbf{\pi }_{G,0}^{\left( 0\right) }=\sigma \mathbf{y}_{0}, \\
\mathbf{\pi }_{G,k}^{\left( 0\right) }=\mathbf{\pi }_{G,k-1}^{\left(
0\right) }\mathbf{R}_{k-1,k},\text{ }1\leq k\leq Z
\end{array
\right.  \label{eq-0Z}
\end{equation
where $\mathbf{y}_{0}$ is the stationary probability vector of the censored
Markov process $U_{0}$ to level $0$, and the scalar $\sigma $ is uniquely
determined by $\sum\nolimits_{k=0}^{Z}\mathbf{\pi }_{G,k}^{\left( 0\right)
\mathbbm{1}=1$. The censored Markov process $T^{\left[ \leq 0\right] }$ has
the stationary probability vector $\mathbf{y}_{0}$ such tha
\begin{equation}
\left\{
\begin{array}{l}
\mathbf{y}_{0}T^{\left[ \leq 0\right] }=0, \\
\mathbf{y}_{0}\mathbbm{1}=1
\end{array
\right.  \label{eq-T0y}
\end{equation}

\subsection{The routing matrix}

The routing probabilities of the closed queueing network may be complicated,
and they depend on the states of each node via considering the bikes'
failure, removing, repair, redistribution and reuse processes under two
batch policies. Also, the routing probabilities for the three classes of
nodes are depicted in Figure 6.
\begin{figure}[h]
\centering                         \includegraphics[width=10cm]{fig-6.eps}
\newline
\caption{Routing probabilities of nodes}
\label{figure:figure-6}
\end{figure}

For the three classes of nodes, we write the routing matrix as $P=\left(
f_{i,j}\right) $, where $f_{i,j}$\ is given b
\begin{equation*}
f_{i,j}=\left\{
\begin{array}{ll}
q_{0,i}, & \text{Node }0\text{ to Node }0\rightarrow i\text{, }1\leq i\leq
\text{,} \\
q_{i,0}, & \text{Node }i\text{ to Node }i\rightarrow 0\text{, }1\leq i\leq
\text{,} \\
q_{i,i}, & \text{Node }i\text{ to Node }i\text{, }0\leq i\leq N\text{,} \\
q_{i,j}, & \text{Node }i\text{ to Node }i\rightarrow j\text{, }j\in \Theta
_{i}\text{, }1\leq i\leq N\text{,} \\
1, & \text{Node }i\rightarrow j\text{ to Node }j\text{, }j\in \Theta _{i
\text{, }1\leq i\leq N\text{,} \\
1, & \text{Node }i\rightarrow 0\text{ to Node }0\text{, }1\leq i\leq N\text{
} \\
1, & \text{Node }0\rightarrow i\text{ to Node }i\text{, }1\leq i\leq N\text{
} \\
0, & \text{otherwise,
\end{array
\right.
\end{equation*
wher
\begin{equation*}
q_{0,0}=1-\left( \sum\limits_{l=0}^{\phi -\psi }\pi _{B,lM;G,Z}^{\left(
0\right) }\right) ,
\end{equation*
for $1\leq i\leq N
\begin{equation*}
q_{0,i}=\left( \sum\limits_{l=0}^{\phi -\psi }\pi _{B,lM;G,Z}^{\left(
0\right) }\right) \beta _{i},\text{ }q_{i,0}=\sum\limits_{k=0}^{K-M}\pi
_{G,k;B,M}^{\left( i\right) },\text{ }q_{i,i}=\sum\limits_{k=0}^{M-1}\pi
_{G,0;B,k}^{\left( i\right) },
\end{equation*
and for $1\leq i\leq N,$ $j\in \Theta _{i}$
\begin{equation*}
q_{i,j}=\left( 1-q_{i,0}-q_{i,i}\right) p_{i,j}.
\end{equation*
Note that the routing probabilities are expressed by the stationary
probabilities of the two Markov processes $\mathbf{Q}$ and $\mathbf{T}$,
both of which are determined by the relative arrival rate vector $\mathbf{e}
. Thus the routing matrix depends on $\mathbf{e}$, and it can be written as
P\left( \mathbf{e}\right) $.

\textbf{A nonlinear routing matrix equation to determine the relative
arrival rates:} Note that the relative arrival rate vector $\mathbf{e}$
satisfies a nonlinear routing matrix equation
\begin{equation}
\mathbf{e}=\mathbf{e}P\left( \mathbf{e}\right) .  \label{ep}
\end{equation
To compute the relative arrival rates, we develop Algorithm 1 to give an
iterative approximate solution to Equation (\ref{ep}).
\begin{algorithm}
\caption{An Iterative Algorithm for Computing the Relative Arrival Rates}\label{alg:1

\begin{algorithmic}
\State \textbf{Step 0: Initialization}

Give an initial vector $_{0}\mathbf{e}$.

\State \textbf{Step 1: The first iterative computation}

(a) Using $_{0}\mathbf{e}$, compute stationary probability vector $_{0
\mathbf{\pi }_{B}^{\left( i\right) }$ for $1\leq i\leq N$ and $_{0}\mathbf
\pi }_{G}^{\left( 0\right) }$ by UL-type RG-factorization.

(b) Obtain the routing matrix $P\left( _{0}\mathbf{e
\right) $ based on $_{0}\mathbf{\pi }_{B}^{\left( i\right) }$ and $_{0}\mathbf{\pi
}_{G}^{\left( 0\right) }$.

(c) Obtain $_{1}\mathbf{e}$ by the equation $_{1}\mathbf{e}= {_{0}\mathbf{e}
P\left( _{0}\mathbf{e}\right) $.

\State \textbf{Step 2: The second iterative computation}

(a) Using $_{1}\mathbf{e}$, compute stationary probability vector $_{1}\pi
_{B}^{\left( i\right) }$ for $1\leq i\leq N$ and $_{1}\mathbf{\pi
_{G}^{\left( 0\right) }$ by UL-type RG-factorization.

(b) Obtain routing matrix $P\left( _{1}\mathbf{e}\right)
$ based on $_{1}\mathbf{\pi }_{B}^{\left( i\right) }$ and $_{1}\mathbf{\pi
}_{G}^{\left( 0\right) }$.

(c) Obtain $_{2}\mathbf{e}$ by the equation $_{2}\mathbf{e}= {_{1}\mathbf{e}
P\left( _{1}\mathbf{e}\right) $.

\State\textbf{Step 3: The $(k+1)$st iterative computation for $k\geq 2$}

(a) Using $_{k}\mathbf{e}$, compute stationary probability vector $_{k
\mathbf{\pi }_{B}^{\left( i\right) }$ for $1\leq i\leq N$ and $_{k}\mathbf
\pi }_{G}^{\left( 0\right) }$ by UL-type RG-factorization.

(b) Obtain routing matrix $P\left( _{k}\mathbf{e}\right) $ based on $_{k
\mathbf{\pi }_{B}^{\left( i\right) }$ and $_{k}\mathbf{\pi }_{G}^{\left(
0\right) }$.

(c) Obtain $_{k+1}\mathbf{e}$ by the equation $_{k+1}\mathbf{e}= {_{k
\mathbf{e}}P\left( _{k}\mathbf{e}\right) $.

\State\textbf{Step 4: Convergence check}

If there exists a relative arrival rate vector $_{n}\mathbf{e}$ such tha
\begin{equation*}
\sqrt{\left( {_{n+1}\mathbf{e}_{0}-}\text{ }_{n}\mathbf{e}_{0}\right)
^{2}+\left( {{_{n+1}\mathbf{e}_{1}-}}\text{ }_{n}\mathbf{e}_{1}\right)
^{2}+\cdots +\left( {_{n+1}\mathbf{e}_{N}-}\text{ }_{n}\mathbf{e}_{N}\right)
^{2}}<\varepsilon ,
\end{equation*
(called a stop condition), for a given precision $\varepsilon =1e-10$, then
the computation stops. In this case, $_{n+1}\mathbf{e}\approx $ $_{n}\mathbf
e}$. Otherwise, return to \textbf{Step 3} until this stop condition is
satisfied.

\State \textbf{Step 5: Output}

Obtain the relative arrival rate vector $_{n}\mathbf{e}$.
\end{algorithmic}
\end{algorithm}

To understand Algorithm 1, we first provide a simple example to show the
routing matrix as well as how to iteratively compute the relative arrival
rates.

\textbf{Example one:} We analyze a DBSS with two parking regions and a
maintenance shop whose structure is depicted in Figure 7.

\begin{figure}[tbp]
\centering                         \includegraphics[width=10cm]{fig-7.eps}
\newline
\caption{A closed queueing network of a DBSS with two parking regions}
\label{figure:figure-7}
\end{figure}

We record the nodes in the following order: Node $0$, Node $1$, Node $2$,
Node $0\rightarrow 1$, Node $1\rightarrow 0$, Node $0\rightarrow 2$, Node
2\rightarrow 0$, Node $1\rightarrow 2$ and Node $2\rightarrow 1$. Then the
routing matrix of the DBSS is given b
\begin{equation*}
P\left( \mathbf{e}\right) =\left(
\begin{array}{ccccccccc}
q_{0,0} &  &  & q_{0,1} &  & q_{0,2} &  &  &  \\
& q_{1,1} &  &  & q_{1,0} &  &  & q_{1,2} &  \\
&  & q_{2,2} &  &  &  & q_{2,0} &  & q_{2,1} \\
& 1 &  &  &  &  &  &  &  \\
1 &  &  &  &  &  &  &  &  \\
&  & 1 &  &  &  &  &  &  \\
1 &  &  &  &  &  &  &  &  \\
&  & 1 &  &  &  &  &  &  \\
& 1 &  &  &  &  &  &  &
\end{array
\right) ,
\end{equation*
where
\begin{equation*}
q_{0,0}=1-\left( \sum\limits_{l=0}^{\phi -\psi }\pi _{B,lM;G,Z}^{\left(
0\right) }\right) ,\text{\ }q_{0,1}=\left( \sum\limits_{l=0}^{\phi -\psi
}\pi _{B,lM;G,Z}^{\left( 0\right) }\right) \beta _{1},\text{ }q_{0,2}=\left(
\sum\limits_{l=0}^{\phi -\psi }\pi _{B,lM;G,Z}^{\left( 0\right) }\right)
\beta _{2},
\end{equation*
\begin{equation*}
q_{1,1}=\sum\limits_{k=0}^{M-1}\pi _{G,0;B,k}^{\left( 1\right) },\text{
q_{1,0}=\sum\limits_{l=0}^{K-M}\pi _{G,k;B,M}^{\left( 1\right) },\text{\
q_{1,2}=\left( 1-q_{1,0}-q_{1,1}\right) p_{1,2},\text{ }
\end{equation*
\begin{equation*}
\text{\ }q_{2,2}=\sum\limits_{k=0}^{M-1}\pi _{G,0;B,k}^{\left( 2\right) }
\text{ }q_{2,0}=\sum\limits_{l=0}^{K-M}\pi _{G,k;B,M}^{\left( 2\right) }
\text{ }q_{2,1}=\left( 1-q_{2,0}-\text{\ }q_{2,2}\right) p_{2,1},\text{\ }
\end{equation*
and the omitted elements are all zeros.

The parameters of Example one are given in Table 1.
\begin{table}[h]
\caption{The parameters in Example one}
\begin{center}
\begin{tabular}{ccccccccccccc}
\toprule $\lambda_1$ & $\lambda_2$ & $\mu_{1,0}$ & $\mu_{2,0}$& $\mu_{1,0}$ & $\mu_{2,0}$ & $\alpha$ & $\beta_1$ & $\beta_2$ & $w$ & $r$ &M&Z\\ \hline
10 & 8& 0.3 & 0.3 & 0.2 & 0.2 & 0.01 & 0.5 & 0.5 & 1 & 2&5 & 10  \\
\bottomrule &  &  &  &  &  &  &  &  &  &  &  &
\end{tabular
\end{center}
\end{table}

By using Algorithm 1, the relative arrival rates of the 9 nodes in the
closed queueing network are given in Table 2.
\begin{table}[h]
\caption{The relative arrival rates are computed by using Algorithm 1}
\begin{center}
\begin{tabular}{ccccccccc}
\toprule $e_{0}$ & $e_{1}$ & $e_{2}$ & $e_{0,1}$ & $e_{1,0}$ & $e_{0,2}$ &
e_{2,0}$ & $e_{1,2}$ & $e_{2,1}$ \\ \hline
0.1400 & 3.8582 & 1.9789 & 0.0152 & 0.0240 & 0.0152 & 0.0065 & 1.4766 &
1.4854 \\
\bottomrule &  &  &  &  &  &  &  &
\end{tabular
\end{center}
\end{table}

\section{A Product-Form Solution and Performance Analysis}

In this section, we provide the product-form solution for the stationary
joint probabilities of queue lengths in the closed queueing network. Based
on this, we can define and analyze some useful performance measures of the
DBSS with unusable bikes.

\subsection{The product-form solution}

Note that the numbers of usable bikes and of unusable bikes in any parking
region and in the maintenance shop are described as the block-structured
Markov processes, and thus the stationary probabilities of Markov processes
of the two classes of nodes can be computed by the UL-type RG-factorization.

It is easy to see that the two classes of nodes in the closed queueing
networks are extended and generalized from simple queueing systems to
block-structured Markov processes.

For the closed queueing network, the following theorem provides the
product-form solution of the stationary joint probabilities $\pi \left(
\mathbf{n}\right) $ of queue lengths for $\mathbf{n}\in \Omega $.

\begin{The}
\label{The:PS} For the closed queueing network of the DBSS with unusable
bikes, the stationary joint probability $\pi \left( \mathbf{n}\right) $ is
given by
\begin{equation*}
\pi \left( \mathbf{n}\right) =\frac{1}{\mathbf{C}}\prod\limits_{i=0}^{N}
\left( n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) \prod\limits
_{\substack{ i=1  \\ j\in \Theta _{i}}}^{N}H\left( m_{i,j}\right)
\prod\limits_{i=1}^{N}H\left( m_{i,0}\right) H\left( m_{0,i}\right) ,
\end{equation*
where
\begin{eqnarray*}
H\left( n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) &=&\left\{
\begin{array}{ll}
\pi _{G,n_{G}^{\left( i\right) };B,n_{B}^{\left( i\right) }}^{\left(
i\right) }, & 1\leq i\leq N,0\leq n_{G}^{\left( i\right) }\leq K,0\leq
n_{B}^{\left( i\right) }\leq M, \\
\pi _{B,n_{B}^{\left( 0\right) };G,n_{G}^{\left( 0\right) }}^{\left(
0\right) }, & \text{ }i=0,0\leq n_{B}^{\left( 0\right) }\leq Z,0\leq
n_{G}^{\left( 0\right) }\leq K
\end{array
\right. \\
H\left( m_{i,j}\right) &=&\left( \frac{e_{i,j}}{\mu _{i,j}}\right) ^{m_{i,j}
\frac{1}{m_{i,j}!},\text{ }0\leq m_{i,j}\leq K, \\
H\left( m_{i,0}\right) &=&\left( \frac{e_{i,0}}{\mu _{i,0}}\right) ^{m_{i,0}
\frac{1}{m_{i,0}!},\text{ }m_{i,0}=kM\text{ for }0\leq k\leq \phi , \\
H\left( m_{0,i}\right) &=&\left( \frac{e_{0,i}}{\mu _{0,i}}\right) ^{m_{0,i}
\frac{1}{\left( m_{0,i}\right) !},\text{ }m_{0,i}=lZ_{i}\text{ for }0\leq
l\leq \frac{\phi }{\psi },
\end{eqnarray*
and $\mathbf{C}$ is a normalization constant, given b
\begin{equation*}
\mathbf{C}=\sum\limits_{\mathbf{n}\in \Omega }\prod\limits_{i=0}^{N}H\left(
n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) \prod\limits
_{\substack{ i=1  \\ j\in \Theta _{i}}}^{N}H\left( m_{i,j}\right)
\prod\limits_{i=1}^{N}H\left( m_{i,0}\right) H\left( m_{0,i}\right) .
\end{equation*}
\end{The}

\textbf{Proof:} We first classify the nodes in the closed queueing network
as three classes: Node $i$ for $0\leq i\leq N$; Node $i\rightarrow j$ for
1\leq i\leq N$ and $j\in \Theta _{i}$; and Node $i\rightarrow 0$ and Node
0\rightarrow i$. Based on this, we prove the product-form solution from the
following three different parts.

\textit{Part one}: Node $i$ for $0\leq i\leq N$. By using the stationary
probabilities of block-structured Markov processes of parking region $i$
with $1\leq i\leq N$ given in Section 4.2, we get
\begin{equation*}
H\left( n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) =\pi
_{G,n_{G}^{\left( i\right) };B,n_{B}^{\left( i\right) }}^{\left( i\right) }.
\end{equation*
Similarly, from the stationary probabilities of block-structured Markov
process of the maintenance shop given in Section 4.3, we obtain
\begin{equation*}
H\left( n_{G}^{\left( 0\right) },n_{B}^{\left( 0\right) }\right) =\pi
_{B,n_{B}^{\left( 0\right) };G,n_{G}^{\left( 0\right) }}^{\left( 0\right) }.
\end{equation*
Thus we have
\begin{equation*}
H\left( n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) =\left\{
\begin{array}{ll}
\pi _{G,n_{G}^{\left( i\right) };B,n_{B}^{\left( i\right) }}^{\left(
i\right) }, & 1\leq i\leq N,0\leq n_{G}^{\left( i\right) }\leq K,0\leq
n_{B}^{\left( i\right) }\leq M, \\
\pi _{B,n_{B}^{\left( 0\right) };G,n_{G}^{\left( 0\right) }}^{\left(
0\right) }, & \text{ }i=0,0\leq n_{B}^{\left( 0\right) }\leq Z,0\leq
n_{G}^{\left( 0\right) }\leq K
\end{array
\right.
\end{equation*}

\textit{Part two}: For Node $i\rightarrow j$ with $1\leq i\leq N$ and $j\in
\Theta _{i}$.

Since the number of servers in Node $i\rightarrow j$ may be infinite,
applying Gorden and Newell \cite{Gor:1967} and Subsection 7.35 in Bolch et
al. \cite{Bol:2006}, the function $H\left( m_{i,j}\right) $ is given by

\begin{equation*}
H\left( m_{i,j}\right) =\left( \frac{e_{i,j}}{\mu _{i,j}}\right) ^{m_{i,j}
\frac{1}{m_{i,j}!},\text{ }0\leq m_{i,j}\leq K.
\end{equation*}

\textit{Part three}: For Node $i\rightarrow 0$ (resp. Node $0\rightarrow i$)
with $1\leq i\leq N$.

The unusable bikes (resp. the repaired bikes) are transported on roads by
the batch removing (resp. batch redistribution) policy. By using Henderson
and Taylor \cite{Hen:1990}, we get
\begin{eqnarray*}
H\left( m_{i,0}\right) &=&\left( \frac{e_{i,0}}{\mu _{i,0}}\right) ^{m_{i,0}
\frac{1}{m_{i,0}!},\text{ }m_{i,0}=kM\text{ for }0\leq k\leq \phi , \\
H\left( m_{0,i}\right) &=&\left( \frac{e_{0,i}}{\mu _{0,i}}\right) ^{m_{0,i}
\frac{1}{\left( m_{0,i}\right) !},\text{ }m_{0,i}=lZ_{i}\text{ for }0\leq
l\leq \frac{\phi }{\psi }.
\end{eqnarray*}

Finally, for the closed queueing network, we can express the product-form
solution in Theorem 2, and the normalization constant is given by
\begin{equation*}
\mathbf{C}=\sum\limits_{\mathbf{n}\in \Omega }\prod\limits_{i=0}^{N}H\left(
n_{G}^{\left( i\right) },n_{B}^{\left( i\right) }\right) \prod\limits
_{\substack{ i=1  \\ j\in \Theta _{i}}}^{N}H\left( m_{i,j}\right)
\prod\limits_{i=1}^{N}H\left( m_{i,0}\right) H\left( m_{0,i}\right) .
\end{equation*
This completes the proof. \textbf{{\rule{0.08in}{0.08in}}}

The following theorem establishes the useful relations between the marginal
stationary probabilities and the joint stationary probability.

\begin{The}
\label{The:DN} For the closed queueing network, the marginal probabilities
of the system can be determined by the joint stationary probabilities as
follows:

(1) For Node $i,1\leq i\leq N,0\leq l\leq M,0\leq k\leq K-l,$
\begin{equation*}
\pi \left( \mathbf{n:}n_{G}^{\left( i\right) }=k,n_{B}^{\left( i\right)
}=l\right) =\sum\limits_{\substack{ \mathbf{n}\in \Omega  \\ n_{G}^{\left(
i\right) }=k,n_{B}^{\left( i\right) }=l}}\pi \left( \mathbf{n}\right) =\pi
_{G,k;B,l}^{\left( i\right) }\frac{\widetilde{\mathbf{C}}\left(
n_{G}^{\left( i\right) }=k,n_{B}^{\left( i\right) }=l\right) }{\mathbf{C}},
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( n_{G}^{\left( i\right) }=k,n_{B}^{\left(
i\right) }=l\right) =\sum\limits_{\substack{ \mathbf{n}\in \Omega  \\
n_{G}^{\left( i\right) }=k,n_{B}^{\left( i\right) }=l}}\prod\limits
_{\substack{ j=0  \\ j\neq i}}^{N}H\left( n_{G}^{\left( j\right)
},n_{B}^{\left( j\right) }\right) \prod\limits_{\substack{ j=1  \\ h\in
\Theta _{j}}}^{N}H\left( m_{j,h}\right) \prod\limits_{j=1}^{N}H\left(
m_{j,0}\right) H\left( m_{0,j}\right) .
\end{equation*}

(2) For Node $i,i=0,0\leq k\leq Z,k+l=hM,0\leq h\leq \phi ,$
\begin{equation*}
\pi \left( \mathbf{n:}n_{G}^{\left( 0\right) }=k,n_{B}^{\left( 0\right)
}=l\right) =\sum\limits_{\substack{ \mathbf{n}\in \Omega  \\ n_{G}^{\left(
0\right) }=k,n_{B}^{\left( 0\right) }=l}}\pi \left( \mathbf{n}\right) =\pi
_{B,l;G,k}^{\left( 0\right) }\frac{\widetilde{\mathbf{C}}\left(
n_{G}^{\left( 0\right) }=k,n_{B}^{\left( 0\right) }=l\right) }{\mathbf{C}},
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( n_{G}^{\left( 0\right) }=k,n_{B}^{\left(
0\right) }=l\right) =\sum\limits_{\substack{ \mathbf{n}\in \Omega  \\
n_{G}^{\left( 0\right) }=k,n_{B}^{\left( 0\right) }=l}}\pro
\limits_{j=1}^{N}H\left( n_{G}^{\left( j\right) },n_{B}^{\left( j\right)
}\right) \prod\limits_{\substack{ j=1  \\ h\in \Theta _{j}}}^{N}H\left(
m_{j,h}\right) \prod\limits_{j=1}^{N}H\left( m_{j,0}\right) H\left(
m_{0,j}\right) .
\end{equation*}

(3) For Node $i\rightarrow j,1\leq i\leq N,j\in \Theta _{i},0\leq h\leq K,$
\begin{equation*}
\pi \left( \mathbf{n:}m_{i,j}=h\right) =\sum\limits_{\substack{ \mathbf{n
\in \Omega  \\ m_{i,j}=h}}\pi \left( \mathbf{n}\right) =\frac{1}{h!}\left(
\frac{e_{i,j}}{\mu _{i,j}}\right) ^{h}\frac{\widetilde{\mathbf{C}}\left(
m_{i,j}=h\right) }{\mathbf{C}},
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( m_{i,j}=h\right) =\sum\limits_{\substack{
\mathbf{n}\in \Omega  \\ m_{i,j}=h}}\prod\limits_{k=0}^{N}H\left(
n_{G}^{\left( k\right) },n_{B}^{\left( k\right) }\right) \prod\limits
_{\substack{ k=1  \\ l\in \Theta _{k}  \\ \left( k,l\right) \neq \left(
i,j\right) }}^{N}H\left( m_{k,l}\right) \prod\limits_{k=1}^{N}H\left(
m_{k,0}\right) H\left( m_{0,k}\right) .
\end{equation*}

(4) For Node $i\rightarrow 0$, $1\leq i\leq N,0\leq h\leq \phi ,$
\begin{equation*}
\pi \left( \mathbf{n:}m_{i,0}=hM\right) =\sum\limits_{\substack{ \mathbf{n
\in \Omega  \\ m_{i,0}=hM}}\pi \left( \mathbf{n}\right) =\frac{1}{\left(
hM\right) !}\left( \frac{e_{i,0}}{\mu _{i,0}}\right) ^{hM}\frac{\widetilde
\mathbf{C}}\left( m_{i,0}=hM\right) }{\mathbf{C}},
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( m_{i,0}=hM\right) =\sum\limits_{\substack{
\mathbf{n}\in \Omega  \\ m_{i,0}=hM}}\prod\limits_{k=0}^{N}H\left(
n_{G}^{\left( k\right) },n_{B}^{\left( k\right) }\right) \prod\limits
_{\substack{ k=1  \\ l\in \Theta _{k}}}^{N}H\left( m_{k,l}\right)
\prod\limits _{\substack{ k=1  \\ k\neq i}}^{N}H\left( m_{k,0}\right)
\prod\limits_{k=1}^{N}H\left( m_{0,k}\right) .
\end{equation*}

(5) For Node $0\rightarrow i,1\leq i\leq N,0\leq l\leq \phi /\psi ,$
\begin{equation*}
\pi \left( \mathbf{n:}m_{0,i}=lZ_{i}\right) =\sum\limits_{\substack{ \mathbf
n}\in \Omega  \\ m_{0,i}=lZ_{i}}}\pi \left( \mathbf{n}\right) =\frac{1}
\left( lZ_{i}\right) !}\left( \frac{e_{0,i}}{\mu _{0,i}}\right) ^{lZ_{i}
\frac{\widetilde{\mathbf{C}}\left( m_{0,i}=lZ_{i}\right) }{\mathbf{C}},
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( m_{0,i}=lZ_{i}\right) =\sum\limits_{\substack{
\mathbf{n}\in \Omega  \\ m_{0,i}=lZ_{i}}}\prod\limits_{k=0}^{N}H\left(
n_{G}^{\left( k\right) },n_{B}^{\left( k\right) }\right) \prod\limits
_{\substack{ k=1  \\ l\in \Theta _{k}}}^{N}H\left( m_{k,l}\right)
\prod\limits_{k=1}^{N}H\left( m_{k,0}\right) \prod\limits_{\substack{ k=1
\\ k\neq i}}^{N}H\left( m_{0,k}\right) .
\end{equation*}
\end{The}

\textbf{Proof}: This proof is obvious by Section 7 in Bolch et al. \cit
{Bol:2006}, and hence it is omitted here. \textbf{{\rule{0.08in}{0.08in}}}

\subsection{Performance Analysis}

In this subsection, we provide some useful performance measures of the DBSS
with unusable bikes by using the stationary joint probability vector $\pi
\left( \mathbf{n}\right) $ for $\mathbf{n}\in \Omega $.

\textbf{(a) The stationary proportion of unusable bikes in the DBSS}

In the DBSS with unusable bikes, we care for the stationary proportion of
unusable bikes. Note that the unusable bikes are distributed in the $N$
parking regions, the maintenance shop, and the roads among various parking
regions and the maintenance shop. The stationary average number of unusable
bikes $E\left[ \Im \right] $ can be computed by means of the marginal
probabilities as follows
\begin{equation*}
E\left[ \Im \right] =\sum\limits_{i=1}^{N}\sum\limits_{l=0}^{M}l\mathbf{\pi
\left( \mathbf{n:}n_{B}^{\left( i\right) }=l\right) +\sum\limits_{l=0}^{K}
\mathbf{\pi }\left( \mathbf{n:}n_{B}^{\left( 0\right) }=l\right)
+\sum\limits_{i=1}^{N}\sum\limits_{l=0}^{\phi }lM\mathbf{\pi }\left( \mathbf
n:}m_{i,0}=lM\right) ,
\end{equation*
where the three terms are the stationary average numbers of unusable bikes
in the $N$ parking regions, the maintenance shop, and Road $i\rightarrow 0$
for $1\leq i\leq N$, respectively. Based on Theorem 3, we obtai
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{B}^{\left( i\right) }=l\right) =\frac{1}
\mathbf{C}}\sum\limits_{k=0}^{K-l}\left[ \pi _{G,k;B,l}^{\left( i\right)
\widetilde{\mathbf{C}}\left( n_{G}^{\left( i\right) }=k,n_{B}^{\left(
i\right) }=l\right) \right] ,
\end{equation*
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{B}^{\left( 0\right) }=l\right) =\frac{1}
\mathbf{C}}\sum\limits_{h=0}^{\phi }\sum\limits_{k=0}^{Z}\sum\limits_{k+l=hM
\left[ \pi _{B,l;G,k}^{\left( 0\right) }\widetilde{\mathbf{C}}\left(
n_{G}^{\left( 0\right) }=k,n_{B}^{\left( 0\right) }=l\right) \right] ,
\end{equation*
an
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}m_{i,0}=lM\right) =\frac{1}{\left( lM\right)
}\left( \frac{e_{i,0}}{\mu _{i,0}}\right) ^{lM}\frac{\widetilde{\mathbf{C}
\left( m_{i,0}=lM\right) }{\mathbf{C}}.
\end{equation*
Thus the stationary proportion of unusable bikes in the DBSS with respect to the total number $K$ of bikes in the DBSS is given b
\begin{equation*}
\eta =\frac{E\left[ \Im \right] }{K}.
\end{equation*}

\textbf{(b)} \textbf{The stationary proportion of usable bikes in either various parking regions or the roads}

We concern the usable bikes in either various parking regions or the roads. Note that the roads with usable bikes are two different cases: The roads with riding bikes, and the roads from the maintenance shop to parking regions. We denote the stationary average number $E\left[ \varpi \right] $ of usable
bikes in either the parking regions or the roads can be
calculated by means of the marginal probabilities as follows:

\begin{equation*}
E\left[ \varpi \right] =\sum\limits_{i=1}^{N}\sum\limits_{k=0}^{K}k\mathbf
\pi }\left( \mathbf{n:}n_{G}^{\left( i\right) }=k\right)
+\sum\limits_{i=1}^{N}\sum\limits_{j\in \Theta _{j}}\sum\limits_{k=0}^{K}
\mathbf{\pi }\left( \mathbf{n:}m_{i,j}=k\right)
+\sum\limits_{i=1}^{N}\sum\limits_{k=0}^{\lfloor K/Z\rfloor }k\beta _{i}
\mathbf{\pi }\left( \mathbf{n:}m_{0,i}=k\beta _{i}Z\right) ,
\end{equation*
where the three terms are the stationary average number of usable bikes in
the $N$\ parking regions, the Road $i\rightarrow j$, and Road $0\rightarrow i$ for $1\leq i\leq N$ and
j\in \Theta _{i}$, respectively. Based on
Theorem 3, we obtai
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{G}^{\left( i\right) }=k\right) =\frac{1}
\mathbf{C}}\sum\limits_{l=0}^{K-k}\left[ \pi _{G,k;B,l}^{\left( i\right)
\widetilde{\mathbf{C}}\left( n_{G}^{\left( i\right) }=k,n_{B}^{\left(
i\right) }=l\right) \right] ,
\end{equation*
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}m_{i,j}=k\right) =\frac{1}{k!}\left( \frac
e_{i,j}}{\mu _{i,j}}\right) ^{k}\frac{\widetilde{\mathbf{C}}\left(
m_{i,j}=k\right) }{\mathbf{C}},
\end{equation*
an
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}m_{0,i}=k\beta _{i}Z\right) =\frac{1}{\left(
k\beta _{i}Z\right) !}\left( \frac{e_{0,i}}{\mu _{0,i}}\right) ^{k\beta
_{i}Z}\frac{\widetilde{\mathbf{C}}\left( m_{0,i}=k\beta _{i}Z\right) }
\mathbf{C}}.
\end{equation*
Thus the stationary proportion of usable bikes in either various parking regions or the roads (with respect to the total number $K$ of bikes in the DBSS) is given b
\begin{equation*}
\xi =\frac{E\left[ \varpi \right] }{K}.
\end{equation*}
Note that the stationary proportion $\xi$ measures the availability of bikes that can directly serve the arriving users.

\textbf{(c) The busy probability of the maintenance shop}

We write $F_{0}$ as the stationary probability that there is no unusable
bike in the maintenance shop. Hence we hav
\begin{equation*}
F_{0}=\mathbf{\pi }\left( \mathbf{n:}n_{B}^{\left( 0\right) }=0\right)
\frac{1}{\mathbf{C}}\sum\limits_{k=0}^{\psi }\left[ \pi _{B,0;G,kM}^{\left(
0\right) }\widetilde{\mathbf{C}}\left( n_{G}^{\left( 0\right)
}=kM,n_{B}^{\left( 0\right) }=0\right) \right] ,
\end{equation*
wher
\begin{equation*}
\widetilde{\mathbf{C}}\left( n_{G}^{\left( 0\right) }=kM,n_{B}^{\left(
0\right) }=0\right) =\sum\limits_{\substack{ \mathbf{n}\in \Omega  \\
n_{G}^{\left( 0\right) }=kM,n_{B}^{\left( 0\right) }=0}}\pro
\limits_{j=1}^{N}H\left( n_{G}^{\left( j\right) },n_{B}^{\left( j\right)
}\right) \prod\limits_{\substack{ j=1  \\ h\in \Theta _{j}}}^{N}H\left(
m_{j,h}\right) \prod\limits_{j=1}^{N}H\left( m_{j,0}\right) H\left(
m_{0,j}\right) .
\end{equation*
Let
\begin{equation*}
F_{A}=1-F_{0}.
\end{equation*
Then the probability $F_{A}$ is the stationary busy probability of the
maintenance shop, that is, the probability that the maintenance shop is busy for
repairing unusable bikes.

\textbf{(d) Two useful stationary proportions}

Let $\gamma _{1}$ be the proportion of the repaired bikes with respect to
all the bikes in the maintenance shop. The
\begin{equation*}
\gamma _{1}=\frac{\sum\limits_{k=0}^{Z}k\mathbf{\pi }\left( \mathbf{n:
n_{G}^{\left( 0\right) }=k\right) }{\sum\limits_{h=0}^{\phi
}\sum\limits_{k+l=hM}(k+l)\mathbf{\pi }\left( \mathbf{n:}n_{G}^{\left(
0\right) }=k,n_{B}^{\left( 0\right) }=l\right) },
\end{equation*
wher
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{G}^{\left( 0\right) }=k\right) =\frac{1}
\mathbf{C}}\sum\limits_{h=0}^{\phi }\sum\limits_{l=hM-k}\left[ \pi
_{B,l;G,k}^{\left( 0\right) }\widetilde{\mathbf{C}}\left( n_{G}^{\left(
0\right) }=k,n_{B}^{\left( 0\right) }=l\right) \right] ,
\end{equation*
an
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{G}^{\left( 0\right) }=k,n_{B}^{\left(
0\right) }=l\right) =\frac{1}{\mathbf{C}}\pi _{B,l;G,k}^{\left( 0\right)
\widetilde{\mathbf{C}}\left( n_{G}^{\left( 0\right) }=k,n_{B}^{\left(
0\right) }=l\right) .
\end{equation*}

Let $\gamma _{2}$ denote the proportion of the unusable bikes in the
maintenance shop with respect to all the unusable bikes in the DBSS
\begin{equation*}
\gamma _{2}=\frac{\sum\limits_{l=0}^{\phi M}l\mathbf{\pi }\left( \mathbf{n:
n_{B}^{\left( 0\right) }=l\right) }{E\left[ \Im \right] },
\end{equation*
wher
\begin{equation*}
\mathbf{\pi }\left( \mathbf{n:}n_{B}^{\left( 0\right) }=l\right) =\frac{1}
\mathbf{C}}\sum\limits_{h=0}^{\phi }\sum\limits_{k=hM-l}\left[ \pi
_{B,l;G,k}^{\left( 0\right) }\widetilde{\mathbf{C}}\left( n_{G}^{\left(
0\right) }=k,n_{B}^{\left( 0\right) }=l\right) \right] .
\end{equation*}

We hope that the two proportions $\gamma _{1}$ and $\gamma _{2}$ are used to
measure the repair ability of the maintenance shop and the removing ability of
unusable bikes in the DBSS, respectively.

\section{Numerical Experiments}

In this section, we provide numerical experiments to discuss how the
performance measures of the DBSS with unusable bikes depend on some key
system parameters by means of our computational method developed in the
previous sections. Furthermore, our intuitive understanding for the numerical experiments can be provided by means of the coupling method of Markov processes. Therefore, such numerical analysis will shed light on
design and operations management of DBSSs with unusable bikes.

We consider a practical DBSS with five parking regions and a maintenance
shop. For example, the physical structure of the DBSS is taken from a real
bike-sharing system with a subway station and four ordinary parking regions which are
closely distributed around the subway station. Based on this, let Node 0 be
the maintenance shop; Node 1 a parking region of the subway station; and
Nodes 2 to 5 the four ordinary parking regions. The closed queueing network structure of the DBSS is
depicted in Figure 8.

Note that the virtual nodes of the DBSS contain three different parts: the
parking regions, the roads and the maintenance shop, thus the virtual closed
queueing network given by us becomes more complicated than the real network structure of
the DBSS with unusable bikes. As such, we provide a simple example of DBSS with unusable bikes for a numerical
illustration purpose. We assume that there are 15 bikes in the DBSS and 2
repairmen in the maintenance shop. It is easy to see from Figure 8 that the
closed queueing network has 24 nodes and 36 routines. Note that the following numerical experiments are completed by means of the
computing power of a laptop.

In what follows we consider two numerical experiments with respect
to two groups of parameters: The failure rate $\alpha $ of usable bikes and the repair rate
w$ of unusable bikes in Subsection 6.1; and the batch sizes $M$ and $Z$ given for the two batch
policies in Subsection 6.2.

\begin{figure}[h]
\centering   \includegraphics[width=11cm]{fig-8.eps} \newline
\caption{The network structure of a DBSS}
\label{figure:figure-8}
\end{figure}

\subsection{Impact of $\protect\alpha $ and $w$ on system
performance}

To conduct the numerical experiments, we set $M=Z=5$ and list the common parameters of the DBSS in Table 3. Now, we apply the coupling method to observe the
effect of system parameters: failure rate $\alpha $ and repair rate $w$, on
the performance measures of the DBSS from three different aspects as follows:

\begin{table}[h]
\caption{The parameters of the DBSS with unusable bikes}
\begin{center}
\begin{tabular}{cccccccccccccc}
\toprule $\lambda _{1}$ & $\lambda _{2}$ & $\lambda _{3}$ & $\lambda _{4}$ &
$\lambda _{5}$ & $\mu _{0,1}$ & $\mu _{0,2}$ & $\mu _{0,3}$ & $\mu _{0,4}$ &
$\mu _{0,5}$ & $\mu _{1,0}$ & $\mu _{1,2}$ & $\mu _{1,3}$ & $\mu _{1,4}$ \\
\midrule 0.30 & 0.35 & 15 & 10 & 10 & 8 & 8 & 0.40 & 0.35 & 0.35 & 0.30 &
0.30 & 0.20 & 0.30 \\
\bottomrule \bottomrule $\mu _{1,5}$ & $\mu _{2,0}$ & $\mu _{2,1}$ & $\mu
_{3,0}$ & $\mu _{3,1}$ & $\mu _{4,0}$ & $\mu _{4,1}$ & $\mu _{5,0}$ & $\mu
_{5,1}$ & $\beta_1$ & $\beta_2$ & $\beta_3$ & $\beta_4$ & $\beta_5$ \\
\midrule 0.35 & 0.20 & 0.30 & 0.20 & 0.30 & 0.20 & 0.35 & 0.20 & 0.35 & 0.2
& 0.2 & 0.2 & 0.2 & 0.2 \\
\bottomrule &  &  &  &  &  &  &  &  &  &  &  &  &
\end{tabular
\end{center}
\end{table}

\textit{(a) The stationary proportion $\eta $ of unusable bikes in the DBSS}

We observe how the stationary proportion $\eta $ of unusable bikes in the
DBSS depends on the failure rate $\alpha $ and the repair rate $w$,
respectively. It is seen from Figure 9 that $\eta $ increases as $\alpha $
increases, but it decreases as $w$ increases. Such results are constant with
our intuition. If the bikes fail fast with a greater $\alpha $, then there
are more unusable bikes in the system, and $\eta $ is larger; if
the unusable bikes are repaired fast with a bigger $w$, then more unusable
bikes can be fast repaired, and $\eta $ becomes smaller.

It is seen from Figure 9(II) that when $w$ is bigger than 1.6, $\eta $
almost keeps at a constant. This implies that when the repair ability of the
maintenance shop exceeds a threshold, our further improving repair ability (e.g., introducing better repair equipments, or adding some repairmen) will not be able to reduce
unusable bikes in the DBSS. This insight is valuable in practical
operations management of DBSSs with unusable bikes.
\begin{figure}[h]
\centering   \includegraphics[width=7.3cm]{etaalpha.eps}
\includegraphics[width=7.3cm]{etaw.eps}
\caption{$\protect\eta $ vs. $\protect\alpha $ and $w$}
\label{figure:figure-9}
\end{figure}

\textit{(b) The stationary proportion $\xi $ of usable bikes in either the parking regions or the roads}

\begin{figure}[h]
\centering   \includegraphics[width=7.3cm]{xialpha.eps}
\includegraphics[width=7.3cm]{xiw.eps}
\caption{$\protect\xi $ vs. $\protect\alpha $ and $w$}
\label{figure:figure-10}
\end{figure}

We observe how the stationary proportion $\xi $ of usable bikes in either the parking regions or the roads depends on the failure rate $\alpha $ and the
repair rate $w$, respectively. It is seen from Figure 10(I)
that $\xi $ decreases as $\alpha $ increases. It is intuitive that when
\alpha $ increases, the total number of usable bikes in either the $N$ parking regions or the roads decrease.

It is seen from Figure 10(II) that $\xi $ increases as $w$ increases; and
when $w$ is bigger than 1.6, $\xi $ almost keeps at a constant. It is
intuitive that if $w$ increases, the total number of usable bikes in either the $N$ parking regions or the roads
increases. Also, when the $w$ exceeds a threshold, our further improving repair
ability (e.g., introducing better repair equipments, or adding some repairmen) will not be able to increase the
number of usable bikes in either the $N$ parking regions and the roads.

\textit{(c) The busy probability of the maintenance shop}

We observe how the busy probability $F_{A}$ of the maintenance shop depends
on the failure rate $\alpha $ and the repair rate $w$, respectively. It
is seen from Figure 11 that $F_{A}$ increases as $\alpha $ increases, and it decreases as $w$ increases.
It can intuitively be understood that if the failure rate of bikes
increases, then more unusable bikes can be fast removed to the maintenance
shop, so that the busy probability of the maintenance shop increases;
whereas if the repair rate $w$ increases, then more unusable bikes can be
fast repaired, so that the busy probability of the maintenance shop
decreases.

It is easy to see from Figure 11(II) that when $w$ is bigger than 1.6,
F_{A} $ changes very slowly. Such an observation is valuable for design of
repair ability in the maintenance shop.
\begin{figure}[h]
\centering \includegraphics[width=7cm]{FAalpha.eps}
\includegraphics[width=7cm]{FAw.eps}
\caption{$F_{A}$ vs. $\protect\alpha $ and $w$}
\label{figure:figure-11}
\end{figure}

\subsection{Impact of $M$ and $Z$ on system performance}

Note that the two batch policies: removing the unusable bikes and
redistributing the repaired bikes, play a key role in the performance of
DBSSs with unusable bikes. Thus we provide some detailed analysis on how the
two batch sizes $M$ and $Z$ influence the performance measures of the DBSS with unusable bikes.

In the numerical experiments, we set $\alpha =0.01$, $w=1$ and provide
some common parameters in Table 4, where $\mu _{0,i}=0$ implies
that the repaired bikes will not be sent to parking region $i$, and hence
\beta _{i}=0$ for $i=2$ to $5$.

Note that the positive integer $Z$ is a multiple of the positive integer $M$
(clearly, $Z\geq M$), and thus once $Z$ is fixed, the possible values of $M$
can also be determined from the simple relation. For example, if $Z=6$, then
$M\in \left\{ 1,2,3,6\right\} $. On the other hand, once $M$ is fixed, $Z=kM
, $k=1,2,3,\ldots $. For example, if $M=2$, then $Z=2,4,6,8,10$. Hereafter,
we write the two batch sizes as a pair $\left( M,Z\right) $. We discuss some
useful pairs as follows: $\left( 1,6\right) ,\left( 2,6\right) $ and $\left(
3,6\right) $ in Figures 13(I), 14(I) and 15(I); and $\left( 2,6\right)
,\left( 2,8\right) $ and $\left( 2,10\right) $ in Figure 13(II), 14(II) and
15(II). Based on this, we can observe how the performance measures of the
DBSS depend on the pair $\left( M,Z\right) $ in the following three
different aspects.

\begin{table}[h]
\caption{The parameters of the DBSS with unusable bikes}
\begin{center}
\begin{tabular}{cccccccccccccc}
\toprule $\lambda _{1}$ & $\lambda _{2}$ & $\lambda _{3}$ & $\lambda _{4}$ &
$\lambda _{5}$ & $\mu _{0,1}$ & $\mu _{0,2}$ & $\mu _{0,3}$ & $\mu _{0,4}$ &
$\mu _{0,5}$ & $\mu _{1,0}$ & $\mu _{1,2}$ & $\mu _{1,3}$ & $\mu _{1,4}$ \\
\midrule 15 & 10 & 10 & 8 & 8 & 0.40 & 0 & 0 & 0 & 0 & 0.20 & 0.30 & 0.30 &
0.35 \\
\bottomrule \bottomrule $\mu _{1,5}$ & $\mu _{2,0}$ & $\mu _{2,1}$ & $\mu
_{3,0}$ & $\mu _{3,1}$ & $\mu _{4,0}$ & $\mu _{4,1}$ & $\mu _{5,0}$ & $\mu
_{5,1}$ & $\beta_{1}$ & $\beta_{2}$ & $\beta_{3}$ & $\beta_{4}$ & $\beta_{5}$
\\
\midrule 0.35 & 0.20 & 0.30 & 0.20 & 0.30 & 0.20 & 0.35 & 0.20 & 0.35 & 1 & 0
& 0 & 0 & 0 \\
\bottomrule &  &  &  &  &  &  &  &  &  &  &  &  &
\end{tabular
\end{center}
\end{table}

\textit{(a) The stationary proportion $\eta $ of unusable bikes in the DBSS}

Figure 12(I) characterizes the effect of $M$ on the
stationary proportion $\eta $. It is easy to see that $\eta $ increases as
M$ increases. It can be intuitively understood that when $M$ becomes bigger, the rate of removing the unusable bikes from any parking region to the maintenance shop will be smaller. Thus it is seen from the coupling method that if $M$ increases,
more unusable bikes are stored in some parking regions, so that the
stationary proportion $\eta $ of unusable bikes in the DBSS increase.

Figure 12(II) shows how the different values of $Z$ affect $\eta $. It is seen that $\eta $ decreases as $Z$ increases. This
result is constant with our intuition as follows: If $Z$ increases, more repaired bikes
are stored in the maintenance shop, so that the number of unusable bikes in the DBSS decreases.

Figure 12 shows the different effect of $M$ and $Z$ on $\eta$. Based on this, we can regulate the
stationary proportion of unusable bikes in the DBSS through changing the two
batch sizes $M$ and $Z$.
\begin{figure}[h]
\centering   \includegraphics[width=7cm]{etaM.eps}
\includegraphics[width=7cm]{etaZ.eps}
\caption{$\protect\eta $ vs. $M$ and $Z$}
\label{figure:figure-12}
\end{figure}

\textit{(b) The stationary proportion $\xi $ of usable bikes in either the parking regions or the roads}

Figure 13(I) characterizes the effect of $M$ on the
stationary proportion $\xi $. It is seen that $\xi $ decreases as $M$
increases. If $M$ increases, then it is seen from the coupling method that more unusable bikes are stored in some
parking regions so that they can not be repaired in time, thus there are fewer usable bikes in either the $N$ parking regions or the roads.

Figure 13(II) shows how the different values of $Z$ affect the stationary
proportion $\xi $. It is seen that $\xi $ increases as $Z$ increases. If $Z$ increases, then it is seen from the coupling method that more repaired bikes are stored in the maintenance
shop, so that there are fewer bikes in either the $N$ parking regions or the roads.

Finally, it may be valuable from Figure 13(I) and (II) shows that the two
batch policies can improve the system performance of the DBSS through
adjusting the two batch sizes $M$ and $Z$.
\begin{figure}[h]
\centering   \includegraphics[width=7cm]{xiM.eps}
\includegraphics[width=7cm]{xiZ.eps}
\caption{$\protect\xi $ vs. $M$ and $Z$}
\label{figure:figure-13}
\end{figure}

\textit{(c) The busy probability of the maintenance shop}

It is seen from Figure 14(I)
that if $Z=6$ is fixed, $F_{A}$ decreases for $M=1,2,3$. It can be
understood by means of the coupling method that if $M$ increases, then more unusable bikes are slowly collected
to the maintenance shop, so that $F_{A}$ decreases. It is seen from Figure
14(II) that if $M=2$ is fixed, then $F_{A}$ increases for $Z=6,8,10$.

From Figure 14, it is seen that $M$ and $Z$ has different influence on the busy probability $F_{A}$ of the maintenance shop. Based on this, we can find a better $F_{A}$ by adjusting to a suitable pair ($M,Z$).

\begin{figure}[h]
\centering   \includegraphics[width=7cm]{FAM.eps}
\includegraphics[width=7cm]{FAZ.eps}
\caption{$F_{A}$ vs. $M$ and $Z$}
\label{figure:figure-14}
\end{figure}

From the above numerical experiments, we can heuristically understand the
effects of system parameters on performance measures of the DBSS with
unusable bikes. Moreover, such numerical results can also motivate us to
intuitively discuss design and operations management of the DBSS with
unusable bikes, for example, finding some better system parameters. In addition,
such numerical analysis can also help us to observe the role played by
the two batch policies: Removing the unusable bikes and redistributing the
repaired bikes, and to find some suitable values of the two batch sizes.
Therefore, our computational method developed in this paper can be a applied to improve the quality of service, system design and operations
management of the DBSSs with unusable bikes.

\section{Concluding Remarks}

This paper studies a DBSS with unusable bikes and develops a new
computational method through applying the RG-factorizations of
block-structured Markov processes in the closed queueing networks. First, we
describe the DBSS with unusable bikes as a virtual closed queueing network
that characterizes the bikes' failure, removing, redistribution and reuse
processes under two batch policies. Then we show that the relative arrival
rates satisfy a nonlinear routing matrix equation. Once the relative arrival
rates are introduced to various nodes, the parking region nodes and the
maintenance shop node are isolated from each other so that their Markov
systems are independent of each other and they are described as the
block-structured Markov processes whose stationary probability vectors can
be computed by the RG-factorizations. Finally, we provide a more general
product-form solution of the closed queueing network through extending and
generalizing from simple queueing systems of the nodes to general
block-structured Markov processes. Based on this, we provide and analyze performance
measures of the DBSS with unusable bikes, and use numerical examples to
verify the validity of our method and to show how the performance measures depend on some key system parameters.

To the best of our knowledge, this paper is the first to apply queueing
networks in the study of DBSSs with unusable bikes. Therefore, we hope that
our methodology and results given in this paper are applicable to analysis
of more general DBSSs with unusable bikes. Along the research line, there
are still a number of interesting directions for future research:

\begin{itemize}
\item In the DBSS with unusable bikes, it is interesting to develop periodic
policies both for removing the unusable bikes and redistributing the
repaired bikes;

\item developing new algorithms both for solving the nonlinear matrix
equation satisfied by the relative arrival rates, and for performance
analysis of the DBSS with unusable bikes;

\item studying periodic arrival processes of users in the DBSS with unusable
bikes in order to observe the influence of the morning or evening peak on
performance of the DBSS; and

\item discussing optimal control issues of the DBSS with unusable bikes.
\end{itemize}

\section*{Acknowledgements}

Quan-Lin Li is supported by the National Natural Science Foundation of China
under grant Nos. 71671158 and 71932002, and by the Natural Science
Foundation of Hebei province under grant No. G2017203277. Xiaole Wu is supported by the National Natural Science Foundation of China under grant Nos. 71622001, 71571047, 71531005, and 91746302.

\section*{Appendix A. The elements of matrix $\mathbf{Q}$}

Now, we provide expression for the block element $Q_{m,n}$ of matrix
\mathbf{Q}$ in Equation (\ref{V-1}), where $m$ and $n$ are two level
variables with $0\leq m,n\leq M$.

We consider parking region $i$. For $0\leq n\leq M-1$
\begin{equation*}
Q_{n,n}=\left(
\begin{array}{ccccc}
-e_{i} & e_{i} &  &  &  \\
\lambda _{i} & -\left( \lambda _{i}+e_{i}+\alpha \right) & e_{i} &  &  \\
& \ddots & \ddots &  &  \\
&  & \lambda _{i} & -\left( \lambda _{i}+e_{i}+\left( K-n-1\right) \alpha
\right) & e_{i} \\
&  &  & \lambda _{i} & -\left( \lambda _{i}+\left( K-n\right) \alpha \right
\end{array
\right) ,
\end{equation*
whose size is $\left( K-n+1\right) \times \left( K-n+1\right) $, and
\begin{equation*}
Q_{n,n+1}=\left(
\begin{array}{ccccc}
0 &  &  &  &  \\
\alpha &  &  &  &  \\
& 2\alpha &  &  &  \\
&  & \ddots &  &  \\
&  &  & \left( K-n-1\right) \alpha &  \\
&  &  &  & \left( K-n\right) \alph
\end{array
\right) _{\left( K-n+1\right) \times \left( K-n\right) };
\end{equation*
for $n=M$

\begin{equation*}
Q_{M,M}=\left(
\begin{array}{ccccc}
-\left( e_{i}+\mu _{i,0}\right) & e_{i} &  &  &  \\
\lambda _{i} & -\left( \lambda _{i}+e_{i}+\mu _{i,0}\right) & e_{i} &  &  \\
& \ddots & \ddots &  &  \\
&  & \lambda _{i} & -\left( \lambda _{i}+e_{i}+\mu _{i,0}\right) & e_{i} \\
&  &  & \lambda _{i} & -\left( \lambda _{i}+\mu _{i,0}\right
\end{array
\right) ,
\end{equation*
whose size is $\left( K-M+1\right) \times \left( K-M+1\right) $, and

\begin{equation*}
Q_{M,0}=\left(
\begin{array}{cccccccc}
\mu _{i,0} &  &  &  &  &  &  &  \\
& \mu _{i,0} &  &  &  &  &  &  \\
&  & \ddots &  &  &  &  &  \\
&  &  & \mu _{i,0} &  &  &  &  \\
&  &  &  & \mu _{i,0} & 0 & \cdots &
\end{array
\right) _{\left( K-M+1\right) \times \left( K+1\right) }.
\end{equation*}

\section*{Appendix B. Proof of Lemma 1}

\textbf{Proof: }It is obvious that $Q_{M,M}^{[\leq M]}=Q_{M,M}$. Now, we use
the inductive method to prove that for the censoring matrix $Q^{[\leq k]}$,
\ref{eq-Qk}) is true for $0\leq k\leq M-1$. Our proof contains three steps
as follows:

Step one. For $k=M-1$, referring to Section 2 of Li \cite{Li:2010}, we get

\begin{eqnarray*}
Q^{[\leq M-1]} &=&\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{M-2,M-2} & Q_{M-2,M-1} \\
&  &  &  & Q_{M-1,M-1
\end{array
\right) \\
&&+\left(
\begin{array}{c}
0 \\
0 \\
\vdots \\
0 \\
Q_{M-1,M
\end{array
\right) \left( -Q_{M,M}\right) ^{-1}\left(
\begin{array}{ccccc}
Q_{M,0} & 0 & \cdots & 0 &
\end{array
\right) \\
&=&\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{M-2,M-2} & Q_{M-2,M-1} \\
Q_{M,0}^{[\leq M-1]} &  &  &  & Q_{M-1,M-1
\end{array
\right) ,
\end{eqnarray*
in which $Q_{M,0}^{[\leq M-1]}=Q_{M-1,M}\left( -Q_{M,M}\right) ^{-1}Q_{M,0}
. Thus when $k=M-1$, (\ref{eq-Qk}) is true.

Step two. We assume that when $k=m$ for $2\leq m\leq M-2$, (\ref{eq-Qk}) is
true, i.e.,

\begin{equation*}
Q^{[\leq k]}=Q^{[\leq m]}=\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{m-1,m-1} & Q_{m-1,m} \\
Q_{m,0}^{[\leq m]} &  &  &  & Q_{m,m
\end{array
\right) ,
\end{equation*
in which $Q_{m,0}^{[\leq m]}=\prod\nolimits_{l=m}^{M-1}\left[
Q_{l,l+1}\left( -Q_{l+1,l+1}\right) ^{-1}\right] Q_{M,0}$. Now, we consider
the case with $k=m-1$ for $2\leq m\leq M-2$
\begin{eqnarray*}
Q^{[\leq m-1]} &=&\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{m-2,m-2} & Q_{m-2,m-1} \\
&  &  &  & Q_{m-1,m-1
\end{array
\right) \\
&&+\left(
\begin{array}{c}
0 \\
0 \\
\vdots \\
0 \\
Q_{m-1,m
\end{array
\right) \left( -Q_{m,m}\right) ^{-1}\left(
\begin{array}{ccccc}
Q_{m,0}^{[\leq m]} & 0 & \cdots & 0 &
\end{array
\right) \\
&=&\left(
\begin{array}{ccccc}
Q_{0,0} & Q_{0,1} &  &  &  \\
& Q_{1,1} & Q_{1,2} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & Q_{m-2,m-2} & Q_{m-2,m-1} \\
Q_{m-1,0}^{[\leq m-1]} &  &  &  & Q_{m-1,m-1
\end{array
\right) ,
\end{eqnarray*
in which $Q_{m,0}^{[\leq m-1]}=Q_{m-1,m}\left( -Q_{m,m}\right)
^{-1}Q_{m,0}^{[\leq m]}$. Since $Q_{m,0}^{[\leq
m]}=\prod\nolimits_{l=m}^{M-1}\left[ Q_{l,l+1}\left( -Q_{l+1,l+1}\right)
^{-1}\right] Q_{M,0}$, it is easy to check that $Q_{m,0}^{[\leq
m-1]}=\prod\nolimits_{l=m-1}^{M-1}\left[ Q_{l,l+1}\left( -Q_{l+1,l+1}\right)
^{-1}\right] Q_{M,0}$. Thus when $k=m-1$, (\ref{eq-Qk}) is also true.
Therefore, (\ref{eq-Qk}) is true for for $1\leq k\leq M-1$.

Step three. It is easy to see from Step two tha
\begin{equation*}
Q^{[\leq 1]}=\left(
\begin{array}{cc}
Q_{0,0} & Q_{0,1} \\
Q_{1,0}^{[\leq 1]} & Q_{1,1
\end{array
\right) ,
\end{equation*
in which $Q_{1,0}^{[\leq 1]}=\prod\nolimits_{l=1}^{M-1}\left[
Q_{l,l+1}\left( -Q_{l+1,l+1}\right) ^{-1}\right] Q_{M,0}$. Then
\begin{eqnarray*}
Q^{[\leq 0]} &=&Q_{0,0}+Q_{0,1}(-Q_{1,1})^{-1}Q_{1,0}^{[\leq 1]} \\
&=&Q_{0,0}+Q_{0,1}(-Q_{1,1})^{-1}\prod\nolimits_{l=1}^{M-1}\left[
Q_{l,l+1}\left( -Q_{l+1,l+1}\right) ^{-1}\right] Q_{M,0}.
\end{eqnarray*}

Based on the above three steps, (\ref{eq-Q0}) is true for $0\leq k\leq M$.
This completes the proof. \textbf{{\rule{0.08in}{0.08in}}}

\section*{Appendix C. Proof of Theorem 1}

\textbf{Proof:} The proof can easily be completed by means of Chapter 2 of
Li \cite{Li:2010}. Here, it is necessary to fix the two special matrices
R_{U\text{ }}$and $G_{L}$ as follows:

\textit{The R-measures}: From (\ref{eq-RG}) and (\ref{eq-R}), we obtain that
$R_{i,j}=Q_{i,j}^{\left[ \leq j\right] }\left( -Q_{j,j}^{[\leq j]}\right)
^{-1}$ for $0\leq i<j\leq M$. It is seen from (\ref{eq-Qk}) that $Q_{i,i+1}^
\left[ \leq i+1\right] }>0$, while the element $Q_{i,j}^{\left[ \leq i+
\right] }$ is zero for $i+2\leq j\leq M$. Thus for $0\leq i<j\leq M$
\begin{equation*}
R_{i,j}=\left\{
\begin{array}{ll}
R_{i,i+1}, & j=i+1, \\
0, & i+2\leq j\leq M
\end{array
\right.
\end{equation*
Thus we hav
\begin{equation*}
R_{U}=\left(
\begin{array}{ccccccc}
0 & R_{0,1} &  &  &  &  &  \\
& 0 & R_{1,2} &  &  &  &  \\
&  & 0 & R_{2,3} &  &  &  \\
&  &  & \ddots & \ddots &  &  \\
&  &  &  & 0 & R_{M-2,M-1} &  \\
&  &  &  &  & 0 & R_{M-1,M} \\
&  &  &  &  &  &
\end{array
\right) .
\end{equation*}

\textit{The G-measures}: From (\ref{eq-RG}) and (\ref{eq-G}), we obtain that
$G_{i,j}=\left( -Q_{i,i}^{[\leq i]}\right) ^{-1}Q_{i,j}^{\left[ \leq i\right]
},$ $0\leq j<i\leq M$. From (\ref{eq-Qk}), we can see $Q_{i,0}^{\left[ \leq
\right] }>0$, while the elements $Q_{i,j}^{\left[ \leq i\right] }$ is zero
for $1\leq j\leq i-1$. Thus for $0\leq j<i\leq M$
\begin{equation*}
G_{i,j}=\left\{
\begin{array}{ll}
G_{i,0}, & j=0, \\
0, & 1\leq j\leq i-1
\end{array
\right.
\end{equation*
Thus we get
\begin{equation*}
G_{L}=\left(
\begin{array}{ccccccc}
0 &  &  &  &  &  &  \\
G_{1,0} & 0 &  &  &  &  &  \\
G_{2,0} &  & 0 &  &  &  &  \\
\vdots &  &  & \ddots &  &  &  \\
G_{M-2,0} &  &  &  & 0 &  &  \\
G_{M-1,0} &  &  &  &  & 0 &  \\
G_{M,0} & 0 & 0 & \cdots & 0 & 0 &
\end{array
\right) .
\end{equation*}

Based on the above analysis, this completes the proof. \textbf{{\rul
{0.08in}{0.08in}}}

\section*{Appendix D. The block elements of Matrix $\mathbf{T}$}

Note that $\overline{w}\left( n_{B}^{\left( 0\right) }\right) =\min \left\{
n_{B}^{\left( 0\right) },r\right\} w$. For $0\leq k\leq \psi -1$ and
Equation (\ref{V-2}), we writ
\begin{equation*}
T_{kM,kM}=\left(
\begin{array}{ccccc}
-e_{0} & e_{0} &  &  &  \\
& -\left( \overline{w}\left( M\right) +e_{0}\right) & e_{0} &  &  \\
&  & \ddots & \ddots &  \\
&  &  & -\left( \overline{w}\left( \left( \phi -k-1\right) M\right)
+e_{0}\right) & e_{0} \\
&  &  &  & -\overline{w}\left( \left( \phi -k\right) M\right
\end{array
\right) ,
\end{equation*
\begin{equation*}
T_{kM,kM+1}=\left(
\begin{array}{ccccc}
0 &  &  &  &  \\
\overline{w}\left( M\right) &  &  &  &  \\
& \overline{w}\left( 2M\right) &  &  &  \\
&  & \ddots &  &  \\
&  &  & \overline{w}\left( \left( \phi -k-1\right) M\right) &  \\
&  &  &  & \overline{w}\left( \left( \phi -k\right) M\right
\end{array
\right) ;
\end{equation*
For $0\leq k\leq \psi -1$, $1\leq j\leq M-1$

\begin{equation*}
T_{kM+j,kM+j}=\left(
\begin{array}{cccc}
-\left( \overline{w}\left( M-j\right) +e_{0}\right) & e_{0} &  &  \\
& \ddots & \ddots &  \\
&  & -\left( \overline{w}\left( \left( \phi -k-1\right) M-j\right)
+e_{0}\right) & e_{0} \\
&  &  & -\overline{w}\left( \left( \phi -k\right) M-j\right
\end{array
\right) ,
\end{equation*
\begin{equation*}
T_{kM+j,kM+j+1}=\left(
\begin{array}{ccccc}
\overline{w}\left( M-j\right) &  &  &  &  \\
& \ddots &  &  &  \\
&  & \overline{w}\left( \left( \phi -k-1\right) M-j\right) &  &  \\
&  &  & \overline{w}\left( \left( \phi -k\right) M-j\right) &
\end{array
\right) ;
\end{equation*}

\begin{equation*}
T_{Z,Z}=\left(
\begin{array}{ccccc}
-\mu _{0} &  &  &  &  \\
& -\mu _{0} &  &  &  \\
&  & \ddots &  &  \\
&  &  & -\mu _{0} &  \\
&  &  &  & -\mu _{0
\end{array
\right) ,
\end{equation*
\begin{equation*}
T_{Z,0}=\left(
\begin{array}{cccccccc}
\mu _{0} &  &  &  &  &  &  &  \\
& \mu _{0} &  &  &  &  &  &  \\
&  & \ddots &  &  &  &  &  \\
&  &  & \mu _{0} &  &  &  &  \\
&  &  &  & \mu _{0} & 0 & \cdots &
\end{array
\right) .
\end{equation*}

redpajama_set_name: RedPajamaArXiv