# One-parameter discrete-time Calogero-Moser system

Umpon Jairuk† and Sikarin Yoo-Kong∗
† _Division of Physics, Faculty of Science and Technology,_
_Rajamangala University of Technology Thanyaburi, Rangsit-Nakornnayok Road,_
_Pathumthani, Thailand 12110._
∗ _The Institute for Fundamental Study(IF), Naresuan University(NU),_
_99 Moo 9, Tha Pho, Mueang Phitsanulok, Phitsanulok, Thailand, 65000_

###### Abstract

We present a new type of integrable one-dimensional many-body systems called a
one-parameter Calogero-Moser (CM) system. In the discrete level, the Lax pairs
with a parameter are introduced and, of course, the discrete-time equations of
motion are obtained as well as their corresponding discrete-time Lagrangian.
The integrability feature of this new system can be captured through the
discrete Lagrangian closure relation by employing a connection with the
temporal Lax matrices of the discrete-time Ruijsenaars-Schneider (RS) system,
exact solution, and the existence of the classical r-matrix. Under the
appropriate limit on the parameter, which in this case is approaching zero,
the standard CM system is retrieved both discrete-time and continuous-time.

##  1 Introduction

The Calogero-Moser (CM) system is a mathematical model that describes the
motion of one-dimensional system of particles interacting through long-range
forces [1, 2]. The CM system is, of course, an integrable system which
exhibits rich symmetries and possesses a sufficient number of conserved
quantities, according to Liouville’s integrability notion, to construct the
exact solutions. Let us give the equations motion of the CM system for the
simplest type of interaction, known as the rational case,

$\ddot{x}_{i}=\sum_{j=1}^{N}\frac{1}{(x_{i}-x_{j})^{3}}\;,\;\;\;i=1,2,3,...,N\;,$
(1.1)

where $x_{i}$ is a position of the $i^{th}$ particle.

The Ruijsenaars-Schneider (RS) system is another integrable one-dimensional
system of particles with a long-range interaction[3, 4]. In the simplest
interaction, namely the rational case, the equations of motion are given by

$\ddot{x}_{i}+\sum_{j=1}^{N}\dot{x}_{i}\dot{x}_{j}\left(\frac{1}{x_{i}-x_{j}+\lambda}+\frac{1}{x_{i}-x_{j}-\lambda}-\frac{2}{x_{i}-x_{j}}\right)\;,\;\;\;i=1,2,3,...,N;,$
(1.2)

where $\lambda$ is a parameter. Under the limit: $\lambda\to 0$, the CM system
is recovered. Then the RS system can be treated as a “one-parameter
generalisation” of the CM system.
In 1994, the time-discretised version of the CM system was introduced by
Nijhoff and Pang [5]. In the rational case, the discrete-time equations of
motion is given by

$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widetilde{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}\right)-\sum\limits_{k=1\atop
k\neq i}^{N}\frac{2}{x_{i}-x_{k}}=0,$ (1.3)

where $\widetilde{x}_{i}=x_{i}(n+1)$ is a forward shift and ${{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{i}}=x_{i}(n-1)$
is a backward shift. The integrability of the system can be captured in the
same sense with the continuous system in terms of the classical r-matrix, the
existence of the exact solution, and the existence of a set of sufficient
invariants.

Soon after, the time-discretised version of the RS system was introduced [6].
In the rational case, the discrete-time equations of motion are given by

$\prod_{j=1\atop jk\neq
i}^{N}\frac{x_{i}-x_{j}+\lambda}{x_{i}-x_{j}-\lambda}=\prod_{j=1}^{N}\frac{(x_{j}-\widetilde{x}_{j})(x_{i}-{{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{j}}+\lambda)}{(x_{j}-{{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{j}})(x_{i}-\widetilde{x}_{j}+\lambda)}\;.$
(1.4)

Again, under the limit: $\lambda\to 0$, the discrete-time CM system is
recovered. Of course, the discrete-time RS system can also be treated as the
“one-parameter generalisation” of the discrete-time CM system.

Recently, a new hallmark for integrability was promoted known as the multi-
dimensional consistency. On the level of the discrete-time equations of
motion, the multi-dimensional consistency can be inferred as the consistency
around the cube [7, 8]. On the level of the Hamiltonians, the feature can be
captured through the Hamiltonian commuting flows as a direct consequence of
the involution in Liouville’s integrability [9]. Alternatively, on the level
of Lagrangians, the multi-dimensional consistency can be expressed through the
Lagrangian closure relation as a direct result in the variation of the action
with respect to independent variables. Since the closure relation for
Lagrangian 1-form will play a major role in this paper as an integrability
criterion, then we shall spend a bit more space to derive the relation. Now
let $\boldsymbol{n}$ be a vector in the lattice and let $\boldsymbol{e}_{i}$
be a unit vector in the $i^{th}$ direction. Then an elementary shift in the
$i^{th}$ direction on the lattice is defined as
$\boldsymbol{n}\to\boldsymbol{n}+\boldsymbol{e}_{i}$. Therefore, the discrete-
time Lagrangians can be expressed in the form

$\mathcal{L}_{i}(\boldsymbol{n})=\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}))\;,$
(1.5)

where $\boldsymbol{x}=\\{x_{1},x_{2},...,x_{N}\\}$. The discrete-time action
is defined as

$S=S[\boldsymbol{x}(\boldsymbol{n}):\Gamma]=\sum_{\boldsymbol{n}\in\Gamma}\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}))\;,$
(1.6)

where $\Gamma$ is a arbitrary discrete curve, see figure 1. Next, we shall
consider another discrete curve $\Gamma^{\prime}$ sharing the same endpoints
with the discrete curve $\Gamma$ and the action is given by

$S^{\prime}=S[\boldsymbol{x}(\boldsymbol{n}):\Gamma^{\prime}]=\sum_{\boldsymbol{n}\in\Gamma^{\prime}}\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}))\;.$
(1.7)

Of course, this can be viewed as the variation of independent variables
$\boldsymbol{n}\to\boldsymbol{n}+\Delta\boldsymbol{n}$ of the action

$\displaystyle S^{\prime}$ $\displaystyle=$ $\displaystyle
S-\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}+\boldsymbol{e}_{j}))+\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}))$
(1.8)
$\displaystyle+\mathcal{L}_{j}(\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}+\boldsymbol{e}_{i}))-\mathcal{L}_{j}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}))\;.$

The least action principle requires $\delta S=S^{\prime}-S=0$ resulting in

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}+\boldsymbol{e}_{j}))-\mathcal{L}_{i}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}))$
(1.9)
$\displaystyle-\mathcal{L}_{j}(\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{i}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}+\boldsymbol{e}_{i}))+\mathcal{L}_{j}(\boldsymbol{x}(\boldsymbol{n}),\boldsymbol{x}(\boldsymbol{n}+\boldsymbol{e}_{j}))\;,$

which is the closure relation for the discrete-time Lagrangian 1-form.
Equivalently, for two-dimensional lattice, see figure 2, (1.9) can be re-
expressed in the form

$\widehat{\mathcal{L}(\boldsymbol{x},\widetilde{\boldsymbol{x}})}-\mathcal{L}(\boldsymbol{x},\widetilde{\boldsymbol{x}})-\widetilde{\mathcal{L}(\boldsymbol{x},\widehat{\boldsymbol{x}})}+\mathcal{L}(\boldsymbol{x},\widehat{\boldsymbol{x}})=0\;.$
(1.10)

$n_{i}$$n_{j}$$\Gamma$$\Gamma^{\prime}$ Figure 1: Arbitrary discrete curves on
the space of independent variables.
$n$$m$$\boldsymbol{x}$$\widetilde{\boldsymbol{x}}$$\widehat{\boldsymbol{x}}$$\widehat{\widetilde{\boldsymbol{x}}}$$\Gamma$$\Gamma^{\prime}$$\widehat{\mathcal{L}(\boldsymbol{x},\widetilde{\boldsymbol{x}})}$$\mathcal{L}(\boldsymbol{x},\widehat{\boldsymbol{x}})$$\mathcal{L}(\boldsymbol{x},\widetilde{\boldsymbol{x}})$$\widetilde{\mathcal{L}(\boldsymbol{x},\widehat{\boldsymbol{x}})}$
Figure 2: The local variation of the discrete curve on the space of two
independent variables.

In this work, we propose a new type of one-parameter CM system, apart from the
RS system and study its integrability through the existence of the exact
solution, classical r-matrix and the closure relation. Therefore, the
structure of the paper is in the following. In section 2, the two compatible
one-parameter discrete-time CM systems will be obtained from the Lax
equations. In section 3, the discrete-time Lagrangians are also established
and the closure relation is directly obtained via the connection between the
RS temporal Lax matrices and the Lagrangian. In section 4, the classical
r-matrix for the one-parameter discrete-time CM system is considered. In
section 5, the construction of the exact solution is carefully derived. In
section 6, the continuum limit will be performed on the one-parameter
discrete-time CM system resulting in the one-parameter continuous-time CM
system. The final section is a summary and possible further investigations.

##  2 One-parameter discrete-time CM system

In this section, we will construct the discrete-time CM system with a
parameter $\lambda$. First, we introduce the spatial Lax matrix:
$\boldsymbol{L}_{\lambda}$ with two temporal matrices: $\boldsymbol{M}$ and
$\boldsymbol{N}$ as follows

$\displaystyle\boldsymbol{L}_{\lambda}$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{N}\frac{1}{x_{i}-x_{j}+\lambda}E_{ij}\;,$ (2.1a)
$\displaystyle\boldsymbol{M}$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{N}\frac{1}{\widetilde{x}_{i}-x_{j}}E_{ij}\;,$
(2.1b) $\displaystyle\boldsymbol{N}$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{N}\frac{1}{\widehat{x}_{i}-x_{j}}E_{ij}\;,$ (2.1c)

where $x_{i}(n,m)$ is the position of the $i^{th}$ particle, $N$ is the number
of particles in the system and $E_{ij}$ is the matrix with entries
$(E_{ij})_{kl}=\delta_{ik}\delta_{jl}$. Here, $\widehat{x}_{i}=x_{i}(m+1)$ is
a forward shift and ${{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{i}}=x_{i}(m-1)$
is a backward shift.

_Discrete flow- $n$ direction_: The compatibility between (2.1a) and (2.1b)
gives us

$\displaystyle\widetilde{\boldsymbol{L}_{\lambda}}\boldsymbol{M}$
$\displaystyle=$ $\displaystyle\boldsymbol{M}\boldsymbol{L}_{\lambda}$
$\displaystyle\sum\limits_{i,j=1}^{N}\sum\limits_{k,\ell=1}^{N}\frac{1}{(\widetilde{x}_{i}-\widetilde{x}_{j}+\lambda)}\frac{1}{(\widetilde{x}_{k}-x_{\ell})}E_{ij}E_{k\ell}$
$\displaystyle=$
$\displaystyle\sum\limits_{i,j=1}^{N}\sum\limits_{k,\ell=1}^{N}\frac{1}{(\widetilde{x}_{i}-x_{j})}\frac{1}{(x_{k}-x_{\ell}+\lambda)}E_{ij}E_{k\ell}\;$
$\displaystyle\sum\limits_{i,\ell=1}^{N}\sum\limits_{k=1}^{N}\frac{1}{(\widetilde{x}_{i}-\widetilde{x}_{k}+\lambda)(\widetilde{x}_{k}-x_{\ell})}E_{i\ell}$
$\displaystyle=$
$\displaystyle\sum\limits_{i,\ell=1}^{N}\sum\limits_{k=1}^{N}\frac{1}{(\widetilde{x}_{i}-x_{k})(x_{k}-x_{\ell}+\lambda)}E_{k\ell}\;.$
(2.2a) Taking out a common factor, we obtain
$\sum\limits_{k=1}^{N}\left(\frac{1}{\widetilde{x}_{i}-\widetilde{x}_{k}+\lambda}-\frac{1}{\widetilde{x}_{i}-x_{k}}\right)=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{k}-x_{\ell}+\lambda}-\frac{1}{\widetilde{x}_{k}-x_{\ell}}\right)\;.$
(2.2b) We see that both sides of (2.2b) are independent and, therefore, it
holds if
$\sum\limits_{k=1}^{N}\left(\frac{1}{\widetilde{x}_{i}-\widetilde{x}_{k}+\lambda}-\frac{1}{\widetilde{x}_{i}-x_{k}}\right)\equiv\widetilde{p}\;,$
(2.2c) where $p=p(n)$ is independent particle indices and a function of
discrete time variable $n$. Taking a backward shift on (2.2c), we have
$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-x_{k}+\lambda}-\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}\right)=p\;.$
(2.2d) Automatically, on the right-hand side of (2.2b), we have
$p=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{\ell}-\widetilde{x}_{k}}-\frac{1}{x_{\ell}-x_{k}-\lambda}\right)\;.$
(2.2e) Now, it is not difficult to see that, from (2.2d) and (2.2e), we obtain
$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widetilde{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}\right)-\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-x_{k}+\lambda}-\frac{1}{x_{i}-x_{k}-\lambda}\right)=0\;,$
(2.2f) which will be treated as a one-parameter discrete-time CM system and,
under the limit: $\lambda\to 0$, one obtains
$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widetilde{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}\right)-\sum\limits_{k=1\atop
k\neq i}^{N}\frac{2}{x_{i}-x_{k}}=0\;,$ (2.2g)

which is nothing but a standard discrete-time CM system in the $n$ direction.

_Discrete flow- $m$ direction_: The compatibility between (2.1a) and (2.1c)
gives us

$\displaystyle\widehat{\boldsymbol{L}_{\lambda}}\boldsymbol{M}$
$\displaystyle=$ $\displaystyle\boldsymbol{M}\boldsymbol{L}_{\lambda}$
$\displaystyle\sum\limits_{i,j=1}^{N}\sum\limits_{k,\ell=1}^{N}\frac{1}{(\widehat{x}_{i}-\widehat{x}_{j}+\lambda)}\frac{1}{(\widehat{x}_{k}-x_{\ell})}E_{ij}E_{k\ell}$
$\displaystyle=$
$\displaystyle\sum\limits_{i,j=1}^{N}\sum\limits_{k,\ell=1}^{N}\frac{1}{(\widehat{x}_{i}-x_{j})}\frac{1}{(x_{k}-x_{\ell}+\lambda)}E_{ij}E_{k\ell}\;$
$\displaystyle\sum\limits_{i,\ell=1}^{N}\sum\limits_{k=1}^{N}\frac{1}{(\widehat{x}_{i}-\widehat{x}_{k}+\lambda)(\widehat{x}_{k}-x_{\ell})}E_{i\ell}$
$\displaystyle=$
$\displaystyle\sum\limits_{i,\ell=1}^{N}\sum\limits_{k=1}^{N}\frac{1}{(\widehat{x}_{i}-x_{k})(x_{k}-x_{\ell}+\lambda)}E_{k\ell}\;.$
(2.3a)

Again, taking out a common factor, we obtain

$\sum\limits_{k=1}^{N}\left(\frac{1}{\widehat{x}_{i}-\widehat{x}_{k}+\lambda}-\frac{1}{\widehat{x}_{i}-x_{k}}\right)=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{k}-x_{\ell}+\lambda}-\frac{1}{\widehat{x}_{k}-x_{\ell}}\right)\;.$
(2.3b)

The situation is similar to the previous discrete flow. Both sides of (2.3b)
are independent and it holds if

$\sum\limits_{k=1}^{N}\left(\frac{1}{\widehat{x}_{i}-\widehat{x}_{k}+\lambda}-\frac{1}{\widehat{x}_{i}-x_{k}}\right)\equiv\widehat{q}\;,$
(2.3c)

where $q=q(m)$ is independent particle indices and a function of discrete time
variable $m$. Computing a backward shift on (2.3c), we obtain

$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-x_{k}+\lambda}-\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{k}}\right)=q\;.$
(2.3d)

From the right-hand side of (2.3b), we shall have

$q=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{\ell}-\widehat{x}_{k}}-\frac{1}{x_{\ell}-x_{k}-\lambda}\right)\;.$
(2.3e)

Therefore, (2.3d) and (2.3e) give

$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widehat{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{k}}\right)-\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-x_{k}+\lambda}-\frac{1}{x_{i}-x_{k}-\lambda}\right)=0\;,$
(2.3f)

which will be treated as a one-parameter discrete-time CM system in the
$m$-direction and, under the limit: $\lambda\to 0$, we obtain

$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widehat{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{k}}\right)-\sum\limits_{k=1\atop
k\neq i}^{N}\frac{2}{x_{i}-x_{k}}=0\;,$ (2.3g)

which is a discrete-time CM system in the $m$ direction.

_Commutativity between discrete flows_ : Two discrete-time dynamics will be
consistent if the compatibility between (2.1b) and (2.1c)

$\displaystyle\widehat{\boldsymbol{M}}\boldsymbol{N}$ $\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{N}}\boldsymbol{M}\;$ (2.4a) holds. This
gives us a set of equations
$p-q=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widetilde{x}_{k}}-\frac{1}{x_{i}-\widehat{x}_{k}}\right)\;,$
(2.4b) and $p-q=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}-\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{k}}\right)\;,$
(2.4c) which will be called corner equations. Imposing (2.4b) = (2.4c), we
obtain
$\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widetilde{x}_{k}}+\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{k}}\right)=\sum\limits_{k=1}^{N}\left(\frac{1}{x_{i}-\widehat{x}_{k}}-\frac{1}{x_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widehat{}}}}}}_{k}}\right)\;,$
(2.4d)

which is a constraint equation given the connection between one discrete flow
and another discrete flow.

##  3 Integrability: the closure relation

In this section, we will show that the one-parameter discrete-time CM systems
in the previous section are integrable in the sense that their discrete-time
Lagrangians satisfy the closure relation as a consequence of the least action
principle with respect to the independent variables [10, 11, 12, 13, 14, 15,
16, 17].

It is not difficult to see that (2.2f) and (2.3f) can be obtained from the
discrete Euler-Lagrange equations [12]

$\displaystyle\widetilde{\frac{\partial\mathcal{L}_{n}(x,\widetilde{x})}{\partial
x_{i}}}+\frac{\partial\mathcal{L}_{n}(x,\widetilde{x})}{\partial{\widetilde{x}_{i}}}=0\;,$
(3.1)
$\displaystyle\widehat{\frac{\partial\mathcal{L}_{m}(x,\widehat{x})}{\partial
x_{i}}}+\frac{\partial\mathcal{L}_{m}(x,\widehat{x})}{\partial{\widehat{x}_{i}}}=0\;,$
(3.2)

where

$\displaystyle\mathcal{L}_{n}{(x,\widetilde{x})}=-\sum\limits_{i,j=1}^{N}\ln\left|x_{i}-\widetilde{x}_{j}\right|+\sum\limits_{i,j=1}^{N}\ln\left|x_{i}-x_{j}+\lambda\right|+p(\Xi-\widetilde{\Xi})\;,$
(3.3)

and

$\displaystyle\mathcal{L}_{m}{(x,\widehat{x})}=-\sum\limits_{i,j=1}^{N}\ln\left|x_{i}-\widehat{x}_{j}\right|+\sum\limits_{i,j=1}^{N}\ln\left|x_{i}-x_{j}+\lambda\right|+q(\Xi-\widehat{\Xi})\;.$
(3.4)

Here $\Xi=\sum\limits_{i=1}^{N}x_{i}$ is a centre of mass variable.

To show that the Lagrangian closure relation for the one-parameter discrete-
time CM model holds, we shall employ a connection between the temporal Lax
matrix and Lagrangian as we did have in the case of the standard discrete-time
CM model [12]. An interesting point is that, for this present system, it turns
out that one could obtain the discrete-time Lagrangian from the relation
$\mathcal{L}(x,\widetilde{x})=\ln\left|\det\boldsymbol{M}_{RS}\right|$111See
the appendix A for the explicit computation., where $\boldsymbol{M}_{RS}$ is a
temporal matrix for the RS model given by

$\displaystyle\boldsymbol{M}_{RS}=\sum\limits_{i,j=1}^{N}\frac{\widetilde{h}_{i}h_{j}}{\widetilde{x}_{i}-x_{j}+\lambda}E_{ij}\;,$
(3.5)

where $h_{i}=h_{i}(n,m)$ are auxiliary variables which can be determined [6].
Suppose there is another temporal matrix given by

$\displaystyle\boldsymbol{N}_{RS}=\sum\limits_{i,j=1}^{N}\frac{\widehat{h}_{i}h_{j}}{\widehat{x}_{i}-x_{j}+\lambda}E_{ij}\;,$
(3.6)

and both $\boldsymbol{M}_{RS}$ and $\boldsymbol{N}_{RS}$ satisfy

$\displaystyle\widehat{\boldsymbol{M}}_{RS}\boldsymbol{N}_{RS}=\widetilde{\boldsymbol{N}}_{RS}\boldsymbol{M}_{RS}.\;$
(3.7)

Taking $\det$ and $\ln$, one obtains

$\displaystyle\ln\left|\det\widehat{\boldsymbol{M}}_{RS}\right|+\ln\left|\det\boldsymbol{N}_{RS}\right|=\ln\left|\det\widetilde{\boldsymbol{N}}_{RS}\right|+\ln\left|\det\boldsymbol{M}_{RS}\right|\;,$
(3.8)

resulting in the closure relation (1.10).

##  4 Integrability: the classical r-matrix

In this section, we shall construct the classical r-matrix for the one-
parameter discrete-time CM system. We first rewrite the spatial Lax matrix as

$\displaystyle\boldsymbol{L}_{\lambda}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{N}\frac{1}{\lambda}E_{ii}-\sum\limits_{i,j=1\atop
j\neq i}^{N}\frac{1}{x_{i}-x_{j}+\lambda}E_{ij}\;.$ (4.1)

Next, we shall call the spatial Lax matrix of the standard CM system [18]
given by

$\displaystyle\boldsymbol{L}=\sum\limits_{i=1}^{N}P_{i}E_{ii}-\sum\limits_{i,j=1\atop
k\neq i}^{N}\frac{1}{x_{i}-x_{j}}E_{ij}\;,$ (4.2)

where $P_{i}$ is the momentum variable for $i^{th}$ particle. With this
structure, one finds that the classical r-matrix can be computed through the
relation

$\displaystyle\\{\boldsymbol{L}\overset{\otimes}{,}\boldsymbol{L}\\}$
$\displaystyle=$
$\displaystyle[r_{12},\boldsymbol{L}\otimes\mathds{1}]-[r_{12},\mathds{1}\otimes\boldsymbol{L}]\;,$
(4.3)

where $r_{12}$ is the classical r-matrix for the CM system. Comparing (4.1)
with (4.2), one immediately finds the classical r-matrix $r_{12}^{\lambda}$
for the one-parameter discrete-time CM system upon replacing
$P_{i}\to\frac{1}{\lambda}$ and
$\frac{1}{x_{i}-x_{j}}\to\frac{1}{x_{i}-x_{j}+\lambda}$

$\displaystyle\\{\boldsymbol{L}_{\lambda}\overset{\otimes}{,}\boldsymbol{L}_{\lambda}\\}$
$\displaystyle=$
$\displaystyle\left[r_{12}^{\lambda},\boldsymbol{L}_{\lambda}\otimes\mathds{1}\right]-\left[r_{12}^{\lambda},\mathds{1}\otimes\boldsymbol{L}_{\lambda}\right]\;.$
(4.4)

We shall note here that under the limit $\lambda\to 0$, the classical r-matrix
$r_{12}^{\lambda}$ will not yield the standard classical r-matrix. This
problem arises from the fact that the spatial Lax matrix (4.1) is a fake one
since it does not provide the integrals of motion through the relation
$I_{n}=\frac{1}{n!}Tr\boldsymbol{L_{\lambda}}^{n}$.

##  5 Integrability: the exact solution

In this section, we will construct the exact solution $\\{x_{i}(n)\\}$ with
initial values $\\{x_{i}(0)\\}$ and $\\{x_{i}(1)=\wtilde{x}_{i}(0)\\}$. We
shall first rewrite the Lax matrices in the forms
$\displaystyle\boldsymbol{X}\boldsymbol{L}-\boldsymbol{L}\boldsymbol{X}+\lambda\boldsymbol{L}=\boldsymbol{E},\;$
(5.1a)
$\displaystyle\widetilde{\boldsymbol{X}}\boldsymbol{M}-\boldsymbol{M}\boldsymbol{X}=\boldsymbol{E}\;,$
(5.1b) where $\boldsymbol{X}=\sum\limits_{i=1}^{N}x_{i}E_{ii}$ and
$\boldsymbol{E}=\sum\limits_{i=1}^{N}E_{ij}$. Moreover, we have
$\displaystyle(\widetilde{\boldsymbol{L}}-\boldsymbol{M})\boldsymbol{E}=0\;,$
(5.1c) and $\displaystyle\boldsymbol{E}(\boldsymbol{L}-\boldsymbol{M})=0\;,$
(5.1d) with also give equations of motion. Let’s write
$\boldsymbol{M}=\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}$ and
$\boldsymbol{L}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}$, where
$\boldsymbol{U}$ is an invertible matrix. (5.1b) leads to
$\displaystyle\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}-\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{X}$
$\displaystyle=$ $\displaystyle\boldsymbol{E}\;$
$\displaystyle\boldsymbol{U}^{-1}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}-\widetilde{\boldsymbol{U}}^{-1}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{X}$
$\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\;$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{U}-\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}$
$\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{U}}-\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}$
$\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\widetilde{\boldsymbol{Y}}-\boldsymbol{Y}$ $\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U},\;$
(5.1e) where $\boldsymbol{Y}=\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}$.
We also find that (5.1a) gives
$\displaystyle\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}-\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{X}+\lambda\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}$
$\displaystyle=$ $\displaystyle\boldsymbol{E}\;$
$\displaystyle\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{U}-\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}+\lambda\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{U}$
$\displaystyle=$ $\displaystyle\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Lambda}-\boldsymbol{U}^{-1}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}+\boldsymbol{U}^{-1}\lambda\boldsymbol{U}\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Lambda}-\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}+\lambda\boldsymbol{U}^{-1}\boldsymbol{U}\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}\boldsymbol{\Lambda}-\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{X}\boldsymbol{U}+\lambda\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\;$
$\displaystyle\boldsymbol{Y}\boldsymbol{\Lambda}-\boldsymbol{\Lambda}\boldsymbol{Y}+\lambda\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\;,$ (5.1f) and
(5.1c) gives
$\displaystyle\left(\widetilde{\boldsymbol{U}}\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}-\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\right)\boldsymbol{E}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\widetilde{\boldsymbol{U}}\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}-\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{E}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\widetilde{\boldsymbol{U}}\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}-\widetilde{\boldsymbol{U}}^{-1}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{E}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}-\boldsymbol{U}^{-1}\boldsymbol{E}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}$
$\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}\;.$
(5.1g) Substituting (5.1g) into (5.1f), we obtain
$\displaystyle\boldsymbol{Y}\boldsymbol{\Lambda}-\boldsymbol{\Lambda}\boldsymbol{Y}+\lambda\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}\;.$
(5.1h) To eliminate the invertible matrix $\boldsymbol{U}$ and
$\boldsymbol{E}$ on the right hand side of (5.1h), we use (5.1d) which can be
expressed in the form
$\displaystyle\boldsymbol{E}\left(\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}-\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\right)$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}-\boldsymbol{E}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{-1}\boldsymbol{U}-\boldsymbol{E}\widetilde{\boldsymbol{U}}\boldsymbol{U}^{-1}\boldsymbol{U}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}-\boldsymbol{E}\widetilde{\boldsymbol{U}}$
$\displaystyle=$ $\displaystyle 0\;$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\widetilde{\boldsymbol{U}}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}\;.$
(5.1i) Since
$\boldsymbol{U}^{-1}\boldsymbol{E}\widetilde{\boldsymbol{U}}=\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}$,
we then obtain
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}$
$\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}\;.$
(5.1j)

Substituting (5.1i) into (5.1e), one finds

$\displaystyle\widetilde{\boldsymbol{Y}}-\boldsymbol{Y}$ $\displaystyle=$
$\displaystyle\boldsymbol{U}^{-1}\boldsymbol{E}\boldsymbol{U}\boldsymbol{\Lambda}\;.\;$
(5.2)

Rearranging (5.1h), we obtain

$\displaystyle\boldsymbol{\Lambda}^{-1}\boldsymbol{Y}\boldsymbol{\Lambda}-\boldsymbol{\Lambda}^{-1}\boldsymbol{\Lambda}\boldsymbol{Y}+\boldsymbol{\Lambda}^{-1}\lambda\boldsymbol{\Lambda}$
$\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}^{-1}\boldsymbol{\Lambda}\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U},\;$
$\displaystyle\boldsymbol{\Lambda}^{-1}\boldsymbol{Y}\boldsymbol{\Lambda}-\boldsymbol{Y}+\lambda$
$\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{U}}^{-1}\boldsymbol{E}\boldsymbol{U}\;.$
(5.3)

Substituting (5.1e) into (5.3), we get

$\displaystyle\boldsymbol{\Lambda}^{-1}\boldsymbol{Y}\boldsymbol{\Lambda}-\boldsymbol{Y}+\lambda$
$\displaystyle=$ $\displaystyle\widetilde{\boldsymbol{Y}}\boldsymbol{Y},\;$
$\displaystyle\widetilde{\boldsymbol{Y}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}^{-1}\boldsymbol{Y}\boldsymbol{\Lambda}+\lambda\;.$
(5.4)

Hence if we proceed $n$ steps, we find that

$\displaystyle\widetilde{\widetilde{\boldsymbol{Y}}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}^{-1}\widetilde{\boldsymbol{Y}}\boldsymbol{\Lambda}+\lambda,\;$
$\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}^{-1}\left[\boldsymbol{\Lambda}^{-1}\boldsymbol{Y}\boldsymbol{\Lambda}+\lambda\right]\boldsymbol{\Lambda}+\lambda,\;$
$\displaystyle=$
$\displaystyle\left(\boldsymbol{\Lambda}^{-1}\right)^{2}\boldsymbol{Y}\boldsymbol{\Lambda}^{2}+2\lambda,\;$
$\displaystyle.\;$ $\displaystyle.\;$ $\displaystyle.\;$
$\displaystyle\boldsymbol{Y}(n)$ $\displaystyle=$
$\displaystyle\left(\boldsymbol{\Lambda}\right)^{-n}\boldsymbol{Y}\boldsymbol{\Lambda}^{n}+n\lambda\;$
(5.5)

and, of course, for the $m$ steps

$\displaystyle\boldsymbol{Y}(m)$ $\displaystyle=$
$\displaystyle\left(\boldsymbol{\Lambda}\right)^{-m}\boldsymbol{Y}\boldsymbol{\Lambda}^{m}+m\lambda.\;$
(5.6)

Then, at any $(n,m)$ steps, we have

$\displaystyle\boldsymbol{Y}(n,m)$ $\displaystyle=$
$\displaystyle\left(p+\boldsymbol{\Lambda}\right)^{-n}\left(q+\boldsymbol{\Lambda}\right)^{-m}\boldsymbol{Y}(0,0)\left(q+\boldsymbol{\Lambda}\right)^{m}\left(p+\boldsymbol{\Lambda}\right)^{n}+(n+m)\lambda.\;$
(5.7)

It is not difficult to find that, under the limit $\lambda\mapsto 0$, one
obtains

$\displaystyle\boldsymbol{Y}(n,m)$ $\displaystyle=$
$\displaystyle\left(p+\boldsymbol{\Lambda}\right)^{-n}\left(q+\boldsymbol{\Lambda}\right)^{-m}\boldsymbol{Y}(0,0)\left(q+\boldsymbol{\Lambda}\right)^{m}\left(p+\boldsymbol{\Lambda}\right)^{n},\;$
(5.8)

which is nothing but a standard solution of the discrete-time CM system [5].

##  6 The continuum limit

In this section, we consider the continuum limit of the one-parameter
discrete-time CM system which had been investigated in the previous sections.
Since there are two discrete-time variables $(n,m)$, we may perform a naive
continuum limit [5] on each of these variables resulting in the one-parameter
continuous-time CM system. To proceed such continuuum limit, we define
$x_{i}=Z_{i}+n\Delta$, where $\Delta$ is a small parameter. Consequently, we
also have $\widetilde{x}_{i}=\widetilde{Z}_{i}+(n+1)\Delta$ and ${\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{x}\limits_{\displaystyle\widetilde{}}}}}}_{i}={\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{i}+(n-1)\Delta$.
Then (2.2f) becomes

$\displaystyle\sum\limits_{k=1}^{N}\left(\frac{1}{Z_{i}-\widetilde{Z}_{k}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{k}+\Delta}\right)-\sum\limits_{i,k=1\atop
k\neq
i}^{N}\left(\frac{1}{Z_{i}-Z_{k}+\lambda}+\frac{1}{Z_{i}-Z_{k}-\lambda}\right)=0\;.$
(6.1)

or

$\displaystyle\left(\frac{1}{Z_{i}-\widetilde{Z}_{i}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{i}+\Delta}\right)-$
$\displaystyle\sum\limits_{i,k=1\atop k\neq
i}^{N}\left(\frac{1}{Z_{i}-\widetilde{Z}_{k}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{k}+\Delta}-\frac{1}{Z_{i}-Z_{k}+\lambda}-\frac{1}{Z_{i}-Z_{k}-\lambda}\right)=0\;.\;$
(6.2)

Taking the expansion, we get

$\displaystyle\widetilde{Z}_{i}$ $\displaystyle=$ $\displaystyle
Z_{i}+\varepsilon\frac{dZ_{i}}{dt}+\frac{\varepsilon^{2}}{2}\frac{d^{2}Z_{i}}{dt^{2}}+...\;,$
(6.3) $\displaystyle{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{i}$
$\displaystyle=$ $\displaystyle
Z_{i}-\varepsilon\frac{dZ_{i}}{dt}+\frac{\varepsilon^{2}}{2}\frac{d^{2}Z_{i}}{dt^{2}}+...\;,$
(6.4)

where $\varepsilon$ is the time-step parameter. Then, the first two terms in
(6.2) can be expressed in the form

$\displaystyle\left(\frac{1}{Z_{i}-\widetilde{Z}_{i}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{i}+\Delta}\right)$
$\displaystyle=$
$\displaystyle\frac{\varepsilon^{2}}{\Delta^{2}}\frac{d^{2}Z_{i}}{dt^{2}}+...\;.$
(6.5)

We also find that

$\displaystyle\sum\limits_{k=1\atop k\neq
i}^{N}\left(\frac{1}{Z_{i}-\widetilde{Z}_{k}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{k}+\Delta}\right)$
$\displaystyle=\sum\limits_{k=1\atop k\neq
i}^{N}\left(\frac{2}{Z_{i}-Z_{k}}+\frac{1}{(Z_{i}-Z_{k})^{3}}\left(\varepsilon^{2}\frac{dZ_{k}}{dt}+2\varepsilon\Delta\frac{dZ_{k}}{dt}+\Delta^{2}\right)+....\right)\;.$
(6.6)

If $\varepsilon\approx\Delta^{2}$, one finds that

$\displaystyle\sum\limits_{k=1\atop k\neq
i}^{N}\left(\frac{1}{Z_{i}-\widetilde{Z}_{k}-\Delta}+\frac{1}{Z_{i}-{\vrule
depth=0.0pt,width=0.0pt{\smash{{\mathop{Z}\limits_{\displaystyle\widetilde{}}}}}}_{k}+\Delta}\right)$
$\displaystyle\approx$ $\displaystyle\sum\limits_{k=1\atop k\neq
i}^{N}\left(\frac{2}{Z_{i}-Z_{k}}+\frac{2\Delta^{2}}{(Z_{i}-Z_{k})^{3}}\right)\;.$
(6.7)

Finally, the continuous version of the one-parameter CM system is given by

$\displaystyle\frac{d^{2}Z_{i}}{dt^{2}}+\sum\limits_{i,k=1\atop k\neq
i}^{N}\left(g^{\prime}\left[\frac{2}{Z_{i}-Z_{k}}-\frac{1}{Z_{i}-Z_{k}+\lambda}-\frac{1}{Z_{i}-Z_{k}-\lambda}\right]+\frac{2g}{(Z_{i}-Z_{k})^{3}}\right)=0\;,$
(6.8)

where $g\equiv\frac{\Delta^{4}}{\varepsilon^{2}}$ and
$g^{\prime}\equiv\frac{\Delta^{2}}{\varepsilon^{2}}$. Therefore, under the
limit: $\lambda\to 0$, we have

$\displaystyle\frac{d^{2}Z_{i}}{dt^{2}}+2g\sum\limits_{i,k=1\atop k\neq
i}^{N}\frac{1}{(Z_{i}-Z_{k})^{3}}$ $\displaystyle=$ $\displaystyle 0\;,$ (6.9)

which is actually a standard continuous CM system. With (6.8), the Lagrangian
is given by

$\displaystyle\mathscr{L}_{\lambda}=\sum\limits_{i=1}^{N}\frac{\partial
Z_{i}}{\partial t}-\frac{1}{2}\sum\limits_{i,k=1\atop k\neq
i}^{N}\frac{g}{(Z_{i}-Z_{k})^{2}}-g^{\prime}\sum\limits_{i,k=1\atop k\neq
i}^{N}\left(\ln\left|Z_{i}-Z_{k}+\lambda\right|+\ln\left|(Z_{i}-Z_{k})\right|\right)$
(6.10)

with the Euler-Lagrange equation

$\displaystyle\frac{\partial\mathscr{L}_{\lambda}}{\partial
Z_{i}}-\frac{\partial}{\partial
t}\left(\frac{\partial\mathscr{L}_{\lambda}}{\partial(\frac{\partial
Z_{i}}{\partial t})}\right)=0\;.$ (6.11)

Of course, under the limit $\lambda\to 0$,

$\lim_{\lambda\mapsto
0}\mathscr{L}_{\lambda}=\mathscr{L}=\sum\limits_{i=1}^{N}\frac{\partial
Z_{i}}{\partial t}+\sum\limits_{i,k=1\atop k\neq
i}^{N}\frac{g}{(Z_{i}-Z_{k})^{2}}\;,$ (6.12)

the standard Lagrangian for the CM system is recovered222We note that the CM
system in this equation comes with the opposite sign with the standard one..
In addition, the Hamiltonian of the one-parameter continuous-time of CM system
can be written in the form

$\displaystyle\mathscr{H}_{\lambda}=\sum\limits_{i=1}^{N}P_{i}^{2}+\frac{1}{2}\sum\limits_{i,k=1\atop
k\neq i}^{N}\frac{g}{(Z_{i}-Z_{k})^{2}}+g^{\prime}\sum\limits_{i,k=1\atop
k\neq
i}^{N}\left(\ln\left|Z_{i}-Z_{k}+\lambda\right|+\ln\left|Z_{i}-Z_{k}\right|\right),\;$
(6.13)

where $P_{i}=\frac{\partial Z_{i}}{\partial t}$ is the momentum variable for
the $i^{th}$ particle.

##  7 Summary

In the present work, we propose a new type of integrable one-dimensional many-
body system called a one-parameter or a deformed discrete-time CM system.
Under the limit: $\lambda\to 0$, a standard CM system is recovered in both
discrete and continuous cases. In figure 3, we provide a diagram of the
connection for all CM-type systems. We would rank our model on the same level
as the RS system since both systems contain a parameter.

One-parameter CM systemRS system$\lambda$ $\to$ $0$$\lambda$ $\to$
$0$$\lambda$ $\to$ $\infty$CM systemGoldfish System Figure 3:  The connection
among one-parameter CM, RS, CM and Goldfish systems.

We also would like to note that the continuous system obtained in section 6 is
just the first one in CM hierarchy [12]. A question can be addressed here is
that “how is the other one deformed in the hierarchy?” Moreover, one also can
try to study the integrability condition as well as the quantum property of
the system. Further investigation is needed and we shall answer these points
elsewhere.

##  Appendix A The connection between Lagrangian and $\boldsymbol{M}_{RS}$ of
the RS model

In this appendix, we will derive the connection between one-parameter
discrete-time Lagrangian and $\boldsymbol{M}_{RS}$

$\displaystyle\boldsymbol{M}_{RS}=\sum\limits_{i,j=1}^{N}\frac{\widetilde{h}_{i}h_{j}}{\widetilde{x}_{i}-x_{j}+\lambda}E_{ij}\;.$
(A.1)

For simplicity, we shall first start with the case of $2\times 2$ matrix given
by

$\boldsymbol{M}_{RS}=\begin{bmatrix}\frac{\widetilde{h}_{1}h_{1}}{\widetilde{x}_{1}-x_{1}+\lambda}&\frac{\widetilde{h}_{1}h_{2}}{\widetilde{x}_{1}-x_{2}+\lambda}\\\
\frac{\widetilde{h}_{2}h_{1}}{\widetilde{x}_{2}-x_{1}+\lambda}&\frac{\widetilde{h}_{2}h_{2}}{\widetilde{x}_{2}-x_{2}+\lambda}\end{bmatrix}\;.$

Then, we compute the determinant

$\displaystyle\det\boldsymbol{M}_{RS}$ $\displaystyle=$
$\displaystyle\frac{\widetilde{h}_{1}h_{1}\widetilde{h}_{2}h_{2}}{(\widetilde{x}_{1}-x_{1}+\lambda)(\widetilde{x}_{2}-x_{2}+\lambda)}-\frac{\widetilde{h}_{2}h_{1}\widetilde{h}_{1}h_{2}}{(\widetilde{x}_{2}-x_{1}+\lambda)(\widetilde{x}_{1}-x_{2}+\lambda)},\;$
$\displaystyle=$ $\displaystyle
h_{1}\widetilde{h}_{1}h_{2}\widetilde{h}_{2}\left[\frac{1}{(\widetilde{x}_{1}-x_{1}+\lambda)(\widetilde{x}_{2}-x_{2}+\lambda)}-\frac{1}{(\widetilde{x}_{2}-x_{1}+\lambda)(\widetilde{x}_{1}-x_{2}+\lambda)}\right],\;$
$\displaystyle=$ $\displaystyle
h_{1}\widetilde{h}_{1}h_{2}\widetilde{h}_{2}\left[\frac{(\widetilde{x}_{2}-x_{1}+\lambda)(\widetilde{x}_{1}-x_{2}+\lambda)-(\widetilde{x}_{1}-x_{1}+\lambda)(\widetilde{x}_{2}-x_{2}+\lambda)}{\prod\limits_{i,j=1}^{2}(\widetilde{x}_{i}-x_{j}+\lambda)}\right].\;$

(LABEL:AOCM) can be further simplified as follows

$\displaystyle\det\boldsymbol{M}_{RS}$ $\displaystyle=$ $\displaystyle
h_{1}\widetilde{h}_{1}h_{2}\widetilde{h}_{2}\left[\frac{(\widetilde{x}_{2}-\widetilde{x}_{1})(x_{1}-x_{2})}{\prod\limits_{i,j=1}^{2}(\widetilde{x}_{i}-x_{j}+\lambda)}\right]\;.$
(A.3)

Recalling the relations [13]

$\displaystyle h^{2}_{i}$ $\displaystyle=$
$\displaystyle-\frac{\prod\limits_{j=1}^{N}(x_{i}-x_{j}+\lambda)(x_{i}-\widetilde{x}_{j}-\lambda)}{\prod\limits_{\mathop{i,j=1}\limits_{j\neq
i}}^{N}(x_{i}-x_{j})\prod\limits_{j=1}^{N}(x_{i}-\widetilde{x}_{j})},\;$ (A.4)
$\displaystyle{\widetilde{h}_{i}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{\prod\limits_{j=1}^{N}(\widetilde{x}_{i}-x_{j}+\lambda)(\widetilde{x}_{i}-\widetilde{x}_{j}-\lambda)}{\prod\limits_{\mathop{i,j=1}\limits_{j\neq
i}}^{N}(\widetilde{x}_{i}-\widetilde{x}_{j})\prod\limits_{j=1}^{N}(\widetilde{x}_{i}-x_{j})}\;,$
(A.5)

then, for $i,j=1,2$, we have

$\displaystyle h^{2}_{1}$ $\displaystyle=$
$\displaystyle-\frac{(x_{1}-x_{1}+\lambda)(x_{1}-x_{2}+\lambda)(x_{1}-\widetilde{x}_{1}-\lambda)(x_{1}-\widetilde{x}_{2}-\lambda)}{(x_{1}-x_{2})(x_{1}-\widetilde{x}_{1})(x_{1}-\widetilde{x}_{2})},\;$
(A.6) $\displaystyle{\widetilde{h}_{1}}^{2}$ $\displaystyle=$
$\displaystyle\frac{(\widetilde{x}_{1}-x_{1}+\lambda)(\widehat{x}_{1}-x_{2}+\lambda)(\widetilde{x}_{1}-\widetilde{x}_{1}-\lambda)(\widetilde{x}_{1}-\widetilde{x}_{2}-\lambda)}{(\widetilde{x}_{1}-\widetilde{x}_{2})(\widetilde{x}_{1}-x_{1})(\widetilde{x}_{1}-x_{2})},\;$
(A.7) $\displaystyle h^{2}_{1}$ $\displaystyle=$
$\displaystyle-\frac{(x_{2}-x_{1}+\lambda)(x_{2}-x_{2}+\lambda)(x_{2}-\widetilde{x}_{1}-\lambda)(x_{2}-\widetilde{x}_{2}-\lambda)}{(x_{2}-x_{1})(x_{2}-\widetilde{x}_{1})(x_{2}-\widetilde{x}_{2})},\;$
(A.8) $\displaystyle{\widetilde{h}_{1}}^{2}$ $\displaystyle=$
$\displaystyle\frac{(\widetilde{x}_{2}-x_{1}+\lambda)(\widehat{x}_{2}-x_{2}+\lambda)(\widetilde{x}_{2}-\widetilde{x}_{1}-\lambda)(\widetilde{x}_{2}-\widetilde{x}_{2}-\lambda)}{(\widetilde{x}_{2}-\widetilde{x}_{1})(\widetilde{x}_{2}-x_{1})(\widetilde{x}_{2}-x_{2})}\;.$
(A.9)

Taking $\ln$, we get

$\displaystyle\ln|h_{1}|$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\ln|\lambda|+\ln|x_{1}-x_{2}+\lambda|+\ln|x_{1}-\widetilde{x}_{1}-\lambda|\right.$
(A.10)
$\displaystyle\left.+\ln|x_{1}-\widetilde{x}_{2}-\lambda|-\ln|x_{1}-x_{2}|-\ln|x_{1}-\widetilde{x}_{1}|-\ln|x_{1}-\widetilde{x}_{2}|\right],\;\;\;$
$\displaystyle\ln|\widetilde{h}_{1}|$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\ln|\lambda|+\ln|\widetilde{x}_{1}-x_{1}+\lambda|+\ln|\widetilde{x}_{1}-x_{2}+\lambda|\right.$
(A.11)
$\displaystyle\left.-\ln|\widetilde{x}_{1}-\widetilde{x}_{2}-\lambda|-\ln|\widetilde{x}_{1}-\widetilde{x}_{2}|-\ln|\widetilde{x}_{1}-x_{1}|-\ln|\widetilde{x}_{1}-x_{2}|\right],\;\;\;$
$\displaystyle\ln|h_{2}|$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\ln|\lambda|+\ln|x_{2}-x_{1}+\lambda|+\ln|x_{2}-\widetilde{x}_{1}-\lambda|\right.$
(A.12)
$\displaystyle\left.+\ln|x_{2}-\widetilde{x}_{2}-\lambda|-\ln|x_{2}-x_{1}|-\ln|x_{2}-\widetilde{x}_{1}|-\ln|x_{2}-\widetilde{x}_{2}|\right],\;\;\;$
$\displaystyle\ln|\widetilde{h}_{2}|$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\ln|\lambda|+\ln|\widetilde{x}_{2}-x_{1}+\lambda|+\ln|\widetilde{x}_{2}-x_{2}+\lambda|\right.$
(A.13)
$\displaystyle\left.+\ln|\widetilde{x}_{2}-\widetilde{x}_{1}-\lambda|-\ln|\widetilde{x}_{2}-\widetilde{x}_{1}|-\ln|\widetilde{x}_{2}-x_{1}|-\ln|\widetilde{x}_{2}-x_{2}|\right]\;.\;\;\;$

Hence,

$\displaystyle\det\boldsymbol{M}_{RS}$ $\displaystyle=$
$\displaystyle\ln|h_{1}|+\ln|\widetilde{h}_{1}|+\ln|h_{2}|+\ln|\widetilde{h}_{2}|+\ln|\widetilde{x}_{2}-\widetilde{x}_{1}|\;$
(A.14)
$\displaystyle+\ln|x_{1}-x_{2}|-\sum\limits_{i,j=1}^{2}\ln|\widetilde{x}_{i}-x_{j}+\lambda|\;$
$\displaystyle=2\ln|\lambda|-\sum\limits_{i,j=1}^{2}\ln|\widetilde{x}_{i}-x_{j}+\lambda|\;$
$\displaystyle=\sum\limits_{i,j=1}^{2}\ln|x_{i}-x_{j}+\lambda|-\sum\limits_{i,j=1}^{2}\ln|x_{i}-\widetilde{x}_{j}|.\;$

Obviously, for the $N$ particles or $N\times N$ matrix, we have

$\det\boldsymbol{M}_{RS}=\sum\limits_{i,j=1}^{N}\ln|x_{i}-x_{j}+\lambda|-\sum\limits_{i,j=1}^{N}\ln|x_{i}-\widetilde{x}_{j}|\;,$
(A.15)

which is indeed the discrete-time Lagrangian for the one-parameter CM system.

Acknowledgements

> Umpon Jairuk would like to thank the Rajamangala University of Technology
> Thanyaburi(RMUTT) for financial support under the Personnel Development Fund
> in 2023.

## References

  * [1] Calogero F, 1975, _Exactly Solvable One-Dimensional Many-Body Problems_ , Lett. nuovo cimento, 13(11), pp. 411-416.
  * [2] Moser J, 1975, _Three Integrable Hamiltonian Systems Connected with Isospectral Deformations_ , Adv. Math, 16, pp. 197-220.
  * [3] Ruijsenaars S M N, 1987, _Complete Integrability of Relativistic Calogero-Moser Systems and Elliptic Function Identities_ , Commun. Math. Phys., 110, pp. 191-213.
  * [4] Schneider H, 1987, _Integrable Relativistic N-Particle Systems in an External Potential_ , Phys. D: Nonlinear Phenom., pp. 203-209.
  * [5] Nijhoff F W and Pang G, 1994, _A Time-Discretized Version of the Calogero-Moser Model_ , Phys. Lett. A, 191, pp. 101-107.
  * [6] Nijhoff F W, Ragnisco O and Kuznetsov V 1996, _Integrable Time-Discretisation of the Ruijsenaars-Schneider Model_ , Commun. Math. Phys., 176, pp. 681-700.
  * [7] Nijhoff F W, and Walker A J, 2001, _The Discrete and Continuous Painleve’ VI Hierarchy and the Garnier Systems_ , Glasgow Math. J., A43, pp. 109-123.
  * [8] Nijhoff F W, 2002, _Lax Pair for the Adler(Lattice Krichever–Novikov) System_ , Phys. Lett., A297, pp. 49-58.
  * [9] Babelon O, Bernard D and Talon M, 2003, _Introduction to Classical Integrable Systems_ , Cambridge University ISBN:9780521822671.
  * [10] Lobb S B, Nijhoff F W and Quispe G R W, 2009, _Lagrangian Multiform Structure for the Lattice KP System_ , J. Phys. A: Math. Theor., 42, 472002.
  * [11] Lobb S B and Nijhoff F W, 2010, _Lagrangian Multiform Structure for the Lattice Gel’Fand–Dikii Hierarchy_ , J. Phys. A: Math. Theor., 43, 072003.
  * [12] Yoo-Kong S, Lobb S and Nijhoff F, 2011, _Discrete-time Calogero–Moser System and Lagrangian 1-form Structure_ , J. Phys. A: Math. Theor., 44, 365203.
  * [13] Yoo-Kong S and Nijhoff F, 2013, _Discrete-time Ruijsenaars-Schneider system and Lagrangian 1-form structure_ , arXiv:1112.4576v2 nlin.SI .
  * [14] Jairuk U, Yoo-Kong S and Tanasittikosol M, 2016, _The Lagrangian Structure of Calogero’s Goldfish Model_ Theoret. and Math. Phys., 183(2), pp.665-683.
  * [15] Jairuk U, Yoo-Kong S and Tanasittikosol M, 2017, _On the Lagrangian 1-form Structure of the Hyperbolic Calogero-Moser System_ , Rep. Math. Phys., 79(3), pp.299-330.
  * [16] Piensuk W and Yoo-Kong S, 2021, _Geodesic Compatibility: Goldfish Systems_ , Rep. Math. Phys., 87(1), pp.45-58.
  * [17] Boll R, Petrera M and Suris Y B, 2015, _Multi-time Lagrangian 1-forms for families of Bäcklund Transformations. Relativistic Toda-type Systems_ , J. Phys. A: Math. Theor., 48 085203\.
  * [18] Avan J and Talon M, 1993, _Classical R-matrix Structure for the Calogero Model_ , Phys. Lett., B303 33-37.


# Mobility-Aware Smart Charging of Electric
Bus Fleets

Ahmadreza Moradipari21, Nathaniel Tucker21, Tuo Zhang2, Gustavo Cezar3 and
Mahnoosh Alizadeh2  2Department of Electrical and Computer Engineering,
University of California, Santa Barbara, California, 93106, USA  3SLAC
National Accelerator Laboratory, GISMo Group, California, 94025, USA  1Authors
have equal contribution

###### Abstract

We study the joint route assignment and charge scheduling problem of a transit
system dispatcher operating a fleet of electric buses in order to maximize
solar energy integration and reduce energy costs. Specifically, we consider a
complex bus transit system with preexisting routes, limited charging
infrastructure, limited number of electric buses, and time-varying electricity
rates. We present a mixed integer linear program (MILP) that yields the
minimal cost daily operation strategy for the fleet (i.e., route assignments
and charging schedules using daily solar forecasts). We present numerical
results from a real-world case study with Stanford University’s Marguerite
Shuttle (a large-scale electric bus fleet) to demonstrate the validity of our
solution and highlight the significant cost savings compared to the status
quo.

##  I Introduction

Due to the potential reduction in operational costs [1], elimination of
tailpipe emissions [2], and encouragement from government agencies [3],
transit systems have started to purchase electric buses over the traditional
diesel or compressed natural gas (CNG) buses. At surface level, replacing
traditional buses with electric buses might seem like a simple task; however,
there are many obstacles preventing a transit system from simply assigning
electric buses to existing routes that were previously served by diesel buses.

The two most fundamental obstacles are the restricted travel distance and
lengthy recharge time of electric buses. Even with recent advances in electric
transportation and battery technology, modern electric buses are commonly
restricted to operate within 20%-95% state of charge (SOC) to prevent
stressing the batteries and reducing lifespan [4]. Combining this SOC
limitation with the high cost of large battery packs, most electric buses are
currently inferior to diesel/CNG buses in operational range. Second, the
recharging process of an electric bus takes significantly more time than the
refueling process of a diesel/CNG bus[4]. Additionally, due to the lengthy
recharge time and limited charging infrastructure, the transit system
dispatcher must be mindful of how the fleet’s recharging infrastructure is
managed in order to provide adequate energy to serve routes.

Despite the aforementioned challenges, the promise of eliminating large
amounts of greenhouse gas emissions from transit buses has enticed early
adopters to operate fleets of electric buses since the early 21st century [1];
however, it is likely that these electric bus fleets are operating
suboptimally in their recharging strategies and route assignments [5].
Accordingly, there has been increasing interest in the optimal operation and
infrastructure planning of electric bus fleets.

The first category of work that studies optimized charging for electric bus
fleets considers the assignment of buses to routes as given, i.e., the times
at which each bus is parked and is available to recharge is predetermined.
Specifically, the authors of [6] present an optimization model for installing
charging infrastructure and sizing batteries for a cost-effective electric bus
fleet. Similarly, the authors of [5] consider infrastructure planning as well
as fleet composition and the recharging process, with the goal of minimizing
total cost of ownership (TOC) of the fleet. Moving away from infrastructure
planning, the authors of [7] present a method to minimize battery aging costs
of an electric bus fleet recharging at nighttime. The authors of [8] present
the cost savings from controlling the charging thresholds for a fleet of
electric buses serving one route continuously in Tallahassee, Florida.
Similarly, the authors of [9] present a MILP framework for scheduling bus
charging and show the potential cost savings from an electric bus fleet in
Davis, California. Furthermore, [10] presents a charging strategy for electric
buses with fast charging infrastructure.

Considering both route assignment and charge scheduling (i.e., the mobility-
aware setting) the authors of [11] present a k-greedy solution method to
maximize travel distance of each electric bus within the fleet. A work similar
to ours, [12], presents a linear formulation for route assignment and charge
scheduling; however, the aim is to minimize the number of electric buses
needed to replace an existing diesel fleet. Hence, the variability of
electricity costs are not considered.

Similar to the aforementioned papers, the work presented in this manuscript
considers both the route assignment and charge scheduling problem of an
electric bus fleet. However, the presented approach is able to improve upon
previous mobility-aware work by accounting for time-varying electricity
prices, utilizing on-site solar energy generation, and providing a minimal
cost schedule for the fleet’s daily operation.

Organization: Section II describes the problem of a fleet dispatcher operating
a fleet of electric buses and proposes a mixed integer linear program (MILP)
formulation that solves for the minimal cost route assignments and recharging
schedule. Section III presents the results of the MILP for the real-world
example of Stanford’s Marguerite Shuttle Transit System.

Figure 1: Primary service area for Stanford University’s Marguerite Shuttle.
Trip origins at Caltrain Palo Alto Transit Center (star). Full system map
available at: https://transportation.stanford.edu/marguerite

##  II Problem Description

We consider a fleet dispatcher attempting to optimize an electric bus transit
system. Specifically, the fleet dispatcher aims to assign electric buses to
serve the daily trips and schedule the recharging of the buses to minimize
electricity cost (e.g., recharging during the inexpensive electricity rates of
nighttime or when solar generation is abundant while still fulfilling all
required bus routes). In the following, we consider the case where the
physical infrastructure (e.g., buses, chargers, parking spots, etc.) and time-
tables (e.g., routes, stops, start/end times, etc.) are already established
within the transit system, but not yet optimized for the aforementioned
objective (as is the case for the Stanford University Marguerite Shuttle,
discussed in Section III). Given the transit system’s fixed time-table and
electric bus infrastructure, the fleet dispatcher seeks to answer questions
such as the following:

  1. 1.

Which electric bus should be assigned to each route at each time?

  2. 2.

When should each electric bus be recharged?

  3. 3.

Does the system need to utilize spare diesel buses to supplement the electric
buses?

  4. 4.

Would more infrastructure benefit the daily operation of the electric bus
fleet?

  5. 5.

What size of on-site solar generation system is needed to fully supply the
fleet with renewable energy?

Let us consider the Stanford Marguerite Shuttle Transit System (Figure 1)
which consists of 38 electric buses, 23 diesel buses, 23 electric bus
chargers, and total of 20 daily routes. Currently, the assignment of buses to
routes and their recharging strategy follows rules adopted by operators that
work well in practice by ensuring sufficient charge is available for service.
However, as we demonstrate in our numerical case study, the current assignment
results in significant losses for the transit system in terms of daily
operational costs and can be improved upon through a joint charge and route
assignment policy. As such, in order to optimize the decision making problem
of the fleet dispatcher, we formulate a mixed-integer-linear-program (MILP) to
solve for both the optimal recharging schedules and route assignments for an
electric bus transit system.

###  II-A MILP Formulation

In the electric bus transit system, we consider one central transit center
(i.e., bus depot) from which all the buses start and finish their routes as
well as recharge. The buses are required to serve numerous routes throughout
the service area, and each route must be served multiple times each day (i.e.,
the electric bus fleet is required to fulfill multiple trips for each route).
We denote $\mathcal{S}$ as the set of scheduled trips across all routes that
need to be fulfilled. For each trip $i\in\mathcal{S}$, let $a_{i}$ and $b_{i}$
denote the start and end time of trip $i$. More specifically, these are the
times that a bus leaves the depot and later returns if serving trip $i$. If
trip $i$ is a one-way route that does not loop back to the depot, we account
for the extra duration for the bus to return to the depot in $b_{i}$
accordingly (i.e., the trip end time $b_{i}$ accounts for “deadhead” travel).
Similarly, if a route does not start at the depot, we account for the deadhead
travel time to the starting location in $a_{i}$.

In order to capture the state of charge of each bus at any time $t$, we
discretize the day into $T$ time steps (e.g., five minute intervals) and
$\mathcal{T}$ is the set of time steps for an entire day. Furthermore, let
$d_{i}$ be the energy consumption per time step for a bus serving trip $i$
(while we assume that varying traffic conditions across different routes can
affect energy consumption rates, we assume that the buses are identical in
their energy consumption when they serve the same route). Let $\mathcal{K}$ be
the set of electric buses and $\mathcal{N}$ be the set of electric bus
chargers installed at the central depot. For each charger $n\in\mathcal{N}$,
$u_{n}$ is the charging rate. Additionally, let ${\bf
p}=[p(t)]_{t\in\mathcal{T}}$ be the vector of electricity prices for an entire
day. We denote as $E^{k}_{min}$ and $E^{k}_{max}$ the minimum and maximum
energy levels for bus $k$, respectively. The fleet dispatcher usually sets
$E^{k}_{min}>0,\forall k\in\mathcal{K}$ for safety precautions. Let $g(t)$ be
the available on-site solar generation at time $t$, which we assume is known
at the time of dispatch. Moreover, we assume that the electricity used from
the on-site solar generation is free for the operator. Last, we denote the
initial energy level of bus $k$ as $e_{0}^{k}$.

Next, we describe the decision variables used in the MILP formulation. We set
the binary variable $X_{i}^{k}(t)$ to $1$ if bus $k$ is serving trip $i$ at
time $t$ and $0$ otherwise. We set the binary variable $Z_{k}(t)$ to $1$ if
bus $k$ is charging at time $t$ and $0$ otherwise. We set the binary variable
$Y_{n}^{k}(t)$ to $1$ if bus $k$ is occupying charger $n$ at time $t$ and $0$
otherwise. We use the variable $E^{k}(t)$ to track the energy level of bus $k$
at time $t$. Lastly, let $V(t)$ be the total amount of electricity that the
dispatcher purchases from the grid at time $t$, and $S(t)$ be the amount of
electricity that buses obtain from the available on-site solar generation at
time $t$. With the necessary notation and decision variables, the joint
charging and routing MILP for the electric bus fleet can be formulated as
follows:

$\displaystyle\text{Minimize }\hskip 10.0pt\sum_{t\in\mathcal{T}}\hskip
10.0ptp(t)V(t)$ (1a) Subject to: $\displaystyle
Z^{k}(t)+\sum_{i\in\mathcal{S}}X^{k}_{i}(t)\leq 1,$ $\displaystyle\forall
k\in\mathcal{K},t\in\mathcal{T}$ (1b)
$\displaystyle\sum_{k\in\mathcal{K}}X_{i}^{k}(t)=1,$ $\displaystyle\forall
i\in\mathcal{S},t\in[a_{i},b_{i}]$ (1c) $\displaystyle
X_{i}^{k}(t+1)=X_{i}^{k}(t),$ $\displaystyle\forall
i\in\mathcal{S},k\in\mathcal{K},t\in[a_{i},b_{i}-1]$ (1d)
$\displaystyle\sum_{k\in\mathcal{K}}Y_{n}^{k}(t)\leq 1,$ $\displaystyle\forall
n\in\mathcal{N},t\in\mathcal{T}$ (1e)
$\displaystyle\sum_{n\in\mathcal{N}}Y_{n}^{k}(t)=Z^{k}(t),$
$\displaystyle\forall k\in\mathcal{K},t\in\mathcal{T}$ (1f) $\displaystyle
E^{k}(t)$
$\displaystyle=E^{k}(t-1)+\sum_{n\in\mathcal{N}}u_{n}Y_{n}^{k}(t)-\sum_{i\in\mathcal{S}}d_{i}X_{i}^{k}(t),$
(1g) $\displaystyle\hskip 89.0pt\forall k\in\mathcal{K},t\in\mathcal{T}$
$\displaystyle\sum_{n\in\mathcal{N}}\sum_{k\in\mathcal{K}}Y_{n}^{k}(t)u_{n}=V(t)+S(t),$
$\displaystyle\forall t\in\mathcal{T}$ (1h) $\displaystyle E^{k}_{min}\leq
E^{k}(t)\leq E^{k}_{max},$ $\displaystyle\forall
k\in\mathcal{K},t\in\mathcal{T}\hskip 58.0pt$ (1i) $\displaystyle
X_{i}^{k}(t)\in\\{0,1\\},$ $\displaystyle\forall
i\in\mathcal{S},k\in\mathcal{K},t\in\mathcal{T}$ (1j) $\displaystyle
Y_{n}^{k}(t)\in\\{0,1\\},$ $\displaystyle\forall
n\in\mathcal{N},k\in\mathcal{K},t\in\mathcal{T}$ (1k) $\displaystyle
Z^{k}(t)\in\\{0,1\\},$ $\displaystyle\forall k\in\mathcal{K},t\in\mathcal{T}$
(1l) $\displaystyle 0\leq S(t)\leq g(t),$ $\displaystyle\forall
t\in\mathcal{T}$ (1m) $\displaystyle E^{k}(0)=e_{0}^{k},$
$\displaystyle\forall k\in\mathcal{K}$ (1n) $\displaystyle
E^{k}(T)=e_{0}^{k},$ $\displaystyle\forall k\in\mathcal{K}.$ (1o)

The objective in equation (1a) aims to minimize the daily electricity cost of
recharging the bus fleet. Constraint (1b) ensures that a bus is either
charging, serving a trip, or parked in the depot (without charging).
Constraint (1c) ensures that all the required daily trips will be served by a
bus. Constraint (1d) ensures that one unique bus will serve each trip (i.e., a
trip cannot be interrupted to switch buses). Constraint (1e) ensures that a
bus can only occupy one charger per time slot. Constraint (1f) guarantees that
if a bus is occupying a charger, then it is charging. Constraint (1g)
calculates the energy level of each bus in each time epoch. Specifically, the
energy level at time $t$ is equal to the energy level at time $t-1$ plus the
charged energy if the bus was charging or minus the spent energy if the bus
was serving a trip. Constraint (1h) ensures that buses obtain electricity from
either the grid or on-site solar. Constraint (1i) ensures that the buses
operate above a desired minimum energy threshold. Constraints (1j)-(1l) are
binary constraints on the decision variables. Constraint (1m) ensures that the
solar energy used by the bus fleet is less than or equal to available solar
generation at time $t$. Lastly, constraint (1n) sets the initial energy of
each bus and constraint (1o) ensures that the energy level of the fleet
returns to the initial value so the same route assignments and charge schedule
can be used for the next day.

###  II-B Behind-the-Meter Solar Integration

To exploit free on-site solar energy and to avoid injecting excess power back
into the distribution grid, the fleet dispatcher prioritizes recharging the
buses during periods when solar generation is available. Only if there is not
enough solar energy, then the fleet dispatcher should purchase electricity
from the grid. As stated in Section II-A, to accommodate behind-the-meter
solar integration, the dispatcher’s MILP formulation makes use of a daily
solar forecast, $g(t)|_{t=1,\dots,T}$. This can be estimated from forecast
models, including those that use weather forecasts, and previous years’ solar
irradiance data. We note that if the solar generation is over-estimated, then
the fleet will have to purchase more expensive grid energy potentially during
peak times such as midday. As such, a conservative estimate is preferred as
cheaper electricity can be procured in the late night period. Future work
could investigate moving-horizon solution methods to account for stochastic
solar generation and update the route and charge assignments in real-time as
solar energy data becomes available.

##  III Case Study

As stated in the introduction, the motivation for the proposed MILP for
electric bus fleets is the real-world Stanford Marguerite Shuttle Transit
System (Figure 1). The Marguerite Shuttle System is free, open to the public,
and operates seven days a week all year traversing the Stanford campus and
surrounding areas. More specific information can be found at
https://transportation.stanford.edu/marguerite.

###  III-A Stanford Marguerite Shuttle System Information

Currently, the Marguerite fleet consists of 23 diesel buses and 38 electric
buses from BYD split into 10 K7 models with battery capacity of 197kWh, 10 K9
models and 18 K9M models, both with 324kWh battery capacity. Additionally, the
central depot is equipped with 23 double port electric bus chargers where each
port can deliver up to 40kW. Each bus can be charged from one or two ports for
a total power of 80kW. For the electricity rates, we consider PG&E’s E-20
electricity rate structure for off-peak, partial-peak, and peak hours. The
electricity rates are given in Table I.

TABLE I: PG&E E-20 Rate Structure Time Interval | Label | Price
---|---|---
12:00am-8:30am | Off-Peak | $0.08422/kWh
8:30am-12:00pm | Partial-Peak | $0.11356/kWh
12:00pm-6:00pm | Peak | $0.16127/kWh
6:00pm-9:30pm | Partial-Peak | $0.11356/kWh
9:30pm-12:00am | Off-Peak | $0.08422/kWh

Furthermore, the Marguerite Shuttle system serves up to 20 unique routes on
any given day. Across all 20 routes, 15 of them are mainly fulfilled by
electric buses, meaning that the electric bus fleet is required to make 352
trips per day, during weekdays. The specific routes and mileages are listed in
Table II. For the purposes of this numerical example, the solar forecast used
was an average daily solar generation calculated from October 2019 with a
maximum generation of 1 MW. The solar forecast is displayed in Figure 2.

TABLE II: Stanford Marguerite Shuttle Route Information Route Name | Daily Trips | Trip Miles
---|---|---
C Line | 33 | 7.00
C Limited | 11 | 4.60
MC Line (AM/PM) | 46 | 3.00
MC Line (Mid Day) | 11 | 5.10
P Line (AM/PM) | 56 | 2.50
P Line (Mid Day) | 11 | 4.00
Research Park (AM/PM) | 24 | 10.40
X Express (AM) | 12 | 1.20
X Line | 44 | 4.60
X Limited (AM) | 10 | 2.00
X Limited (PM) | 10 | 1.50
Y Express (PM) | 20 | 1.20
Y Line | 44 | 4.60
Y Limited (AM) | 10 | 2.40
Y Limited (PM) | 10 | 2.00
Totals | 352 trips/day | 1431.50 miles/day
Figure 2: Average daily solar generation for a 1 MW on-site installation. Data
averaged from CAISO renewable database in October 2019.

###  III-B Simulation Results

The proposed MILP was implemented in Matlab making use of CVX and Mosek. All
numerical experiments were run on a laptop with 16 GB of RAM and 3.5 GHz Intel
i7 processor. This section reports on the charging schedule, route
assignments, and cost savings when comparing the proposed MILP solution with
on-site solar generation, without on-site solar generation, and the status quo
(i.e., the status quo is the actual operations of the Stanford Marguerite
Fleet from 7-October-2019) which does not yet exploit free on-site solar
generation.

Figure 3 presents the energy levels of each bus in the fleet during the day
when the dispatch is generated through our proposed MILP. Time on the x-axis
begins at 5:00am, as this is the start of the earliest route that must be
fulfilled. The left plot shows the energy levels of the buses when the MILP is
not utilizing on-site solar generation. The right plot shows the battery
levels of the buses when the MILP accounts for on-site solar generation. It
will become more clear when examining Figure 4 that the buses charge more
during midday in the right plot than the left, to make use of the free on-site
solar.

Figure 3: Left: Battery levels for each electric bus when considering a fleet
without available on-site solar. Right: Battery levels for each electric bus
when optimizing with available on-site solar generation.

Figure 4 presents the total charging power of the fleet across the entire day.
The red curve presents the total charging power for the MILP solution that
does not exploit on-site solar generation. Conversely, the blue plot shows the
fleet’s total charging power from the MILP solution that does account for on-
site solar generation. It is clear from this plot that the solution that
accounts for on-site solar (blue) is able to charge in the middle of the day
when solar is abundant; however, the solution that does not exploit solar
(red) does not charge during the midday as the electricity prices are highest
at this time. Instead, the fleet has a spike in charging power in the evening
when electricity rates are decreased. This large transient in the evening
could be detrimental to grid stability, increase in harmonics, accelerate
aging of grid assets (i.e. transformers) and could potentially lead to demand
charges for the fleet dispatcher due to high power consumption. As such, the
solution making use of on-site solar generation with a forecasting method is
preferable.

Figure 4: Total charging power of the fleet throughout the day. Blue: Solution
accounting for on-site solar generation. Red: Solution does not include on-
site solar generation.

Last, Figure 5 presents the daily electricity costs for the three different
test cases. Case A: Status Quo. We had access to the data from the operations
of the Stanford Marguerite fleet on 7-October-2019 and calculated the cost of
charging the fleet under the E-20 rate structure. As such, under normal
operation, the daily operational cost was $715.10 USD. Case B corresponds to
the solution of the proposed MILP with the same routes, buses, and chargers as
Case A; however, the mobility-aware solution reassigned buses to new trips and
rescheduled the charging of each bus. In Case B, the MILP solution did not
account for on-site solar and the daily cost was $267.90 USD. Last, Case C was
identical to Case B; however, the MILP accounted for the on-site solar
generation and had access to the daily solar forecast. As such, the daily cost
was reduced to $61.89 USD. From these results, it is evident that the fleet
dispatcher benefits from the MILP formulation for routing and charging ($55\%$
decrease in cost in Case B).

Figure 5: Price Comparison for 3 difference regimes: Case 1: Status Quo,
electric bus charging data obtained from real-implementation (Stanford
Marguerite Shuttle) on 7-Oct-2019. Case 2: Mobility-Aware MILP solution for
same routes and buses as Case A, without on-site solar generation. Case 3:
Mobility-Aware MILP solution for same routes and buses as Case A, with on-site
solar generation.

##  IV Conclusion

In this paper, we investigated the joint route assignment and charge
scheduling problem of a transit system dispatcher operating a fleet of
electric buses in order to maximize solar energy integration and reduce energy
costs. We considered a complex bus transit system with preexisting routes,
limited charging infrastructure, limited number of electric buses, and time-
varying electricity rates. We presented a mixed integer linear program (MILP)
that yields route assignments and charging schedules using daily solar
forecasts. We presented numerical results from a real-world case study with
Stanford University’s Marguerite Shuttle to demonstrate the cost-saving
benefits of our solution and highlight the significant cost savings compared
to the status quo.

Future work includes investigating a moving-horizon solution approach to
account for stochastic solar generation. Additionally, we would like to add
traditional diesel routes to the optimization to further minimize emissions
and to expand the clean operation of the electric bus fleet. Further future
work can include performing field test experiments with real buses during
operational hours, determining the optimal solar capacity to fully charge the
electric bus fleet, and quantify the value and size of onsite solar and
battery combination for resiliency.

## Acknowledgment

The authors would like to thank the Stanford Transportation team for the
support, discussions, and information about operations. This work was funded
by the California Energy Commission under grant EPC-17-020. SLAC National
Accelerator Laboratory is operated for the US Department of Energy by Stanford
University under Contract DE-AC02-76SF00515.

## References

  * [1] J. Horrox and M. Casale, “Electric Buses: Lessons from Cities Pioneering Clean Transportation,” October 2019, U.S. PIRG Education Fund, Environment America Research & Policy Center, Frontier Group.
  * [2] “Greenhouse gases, Regulated Emissions, and Energy use in Transportation (GREET) Model,” Argonne National Laboratory. https://greet.es.anl.gov/.
  * [3] “California transitioning to all-electric public bus fleet by 2040,” December 2018, California Air Resources Board. https://ww2.arb.ca.gov/.
  * [4] “Battery Bus Range - It’s All in the Math,” Mass Transit. https://www.masstransitmag.com/bus/article/12131451/battery-bus-range-its-all-in-the-math.
  * [5] M. Rogge, E. van der Hurk, A. Larsen, and D. U. Sauer, “Electric bus fleet size and mix problem with optimization of charging infrastructure,” _Applied Energy_ , vol. 211, pp. 282 – 295, 2018. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0306261917316355
  * [6] A. Kunith, R. Mendelevitch, and D. Goehlich, “Electrification of a city bus network—an optimization model for cost-effective placing of charging infrastructure and battery sizing of fast-charging electric bus systems,” _International Journal of Sustainable Transportation_ , vol. 11, no. 10, pp. 707–720, 2017. [Online]. Available: https://doi.org/10.1080/15568318.2017.1310962
  * [7] A. Houbbadi, R. Trigui, S. Pelissier, E. Redondo-Iglesias, and T. Bouton, “Optimal scheduling to manage an electric bus fleet overnight charging,” _Energies_ , vol. 12, no. 14, p. 2727, 2019.
  * [8] N. Qin, A. Gusrialdi, R. P. Brooker, T. Ali _et al._ , “Numerical analysis of electric bus fast charging strategies for demand charge reduction,” _Transportation Research Part A: Policy and Practice_ , vol. 94, pp. 386–396, 2016.
  * [9] Y. Wang, Y. Huang, J. Xu, and N. Barclay, “Optimal recharging scheduling for urban electric buses: A case study in davis,” _Transportation Research Part E: Logistics and Transportation Review_ , vol. 100, pp. 115–132, 2017.
  * [10] H. Chen, Z. Hu, Z. Xu, J. Li, H. Zhang, X. Xia, K. Ning, and M. Peng, “Coordinated charging strategies for electric bus fast charging stations,” in _2016 IEEE PES Asia-Pacific Power and Energy Engineering Conference (APPEEC)_. IEEE, 2016, pp. 1174–1179.
  * [11] T. Paul and H. Yamada, “Operation and charging scheduling of electric buses in a city bus route network,” in _17th International IEEE Conference on Intelligent Transportation Systems (ITSC)_ , Oct 2014, pp. 2780–2786.
  * [12] M. Janovec and M. Koháni, “Exact approach to the electric bus fleet scheduling,” _Transportation Research Procedia_ , vol. 40, pp. 1380–1387, 2019.