text:
\section{Introduction}

The Belle Collaboration, together with BaBar, has made essential
contributions to establish the theory of Kobayashi and Maskawa, who
explain all known CP violation phenomena within the framework of the
Standard Model (SM) by a single irreducible phase appearing in the
quark mixing matrix. For their outstanding achievement Kobayashi and
Maskawa were awarded the Nobel Prize of 2008. Although the SM has
been extremely
successful in describing virtually all data, most importantly the CP
violation phenomena of the $K$- and the $B$-systems (for a
comprehensive overview see, e.g.~\cite{buras}), there are a
number of arguments why the SM cannot be the regarded as a complete
theory. In fact, there is clear evidence for physics beyond the SM, as
suggested by the non-vanishing mass of the neutrinos, the
matter-antimatter asymmetry observed in the universe, and the apparent
necessity for dark matter. Most likely, the effects mentioned have to
do with CP violation of a yet undiscovered source. The ``New Physics'' (NP)
generating these sources is expected to appear at large, so far
unreached (multi-TeV) energy scales. While the discovery and
exploration of New
Physics is the central motivation for the LHC program, flavor physics is
expected to play a key role in unraveling possible NP at this scale
and to solve the puzzle of CP violation.

Colliding $e^+$ and $e^-$ beams with different (``asymmetric'') energies
to produce the $\Upsilon(4S)$ resonance, with
just enough energy for the creation of a pair of $B$-mesons (or background
with other quark flavors), is an alternative approach to the
high-energy frontier experiments at the LHC. With high luminosity, and
consequently large statistics, as achieved in the next generation of
flavor factories (``Super Flavor Factory'',
SFF)~\cite{superkekb,superb}, very large energy
scales can be reached, when quantum loop corrections to the SM are
considered (see fig.~\ref{fig:penguin}). Depending on the flavor
changing couplings of the NP particle spectrum, the sensitivity to large
mass scales in a SFF may be from many hundred GeV
up to tens of TeV. In this respect a SFF is truely complementary to
the LHC. A recent review of the physics potential of a future high
luminosity B factory can be found, e.g., in~\cite{newphys}. One should note
here that the discovery potential of a future SFF is indeed
extraordinary and might reach even beyond the LHC.

There are several distinct approaches to look for NP at
SFF's: Many of them concentrate on the precise measurement of the
angles and sides of the unitarity triangle for $B$-meson couplings,
derived from the Cabibbo-Kobayashi-Maskawa (CKM) matrix which connects
the mass eigenstates to the flavor eigenstates of the
down-type quarks. Within the SM the CKM matrix is unitary and their
elements describe the coupling strengths of the flavor changing
currents at the quark level, such as $b \rightarrow c$ or $b \rightarrow u$. The
unitarity of the CKM
matrix gives rise to a total of six
so-called unitarity triangles, one of them involving the $b$-quark
couplings (``$B$-triangle''). All triangles have the same area, but the
$B$-triangle has all three sides of the same order, corresponding to
large angles and thus giving rise to large CP violating effects. Within
the SM the
$B$-triangle is highly over-constrained (5 observables for only 2
independent quantities), so a precise measurement of all the three
angles and the two sides is a crucial check of the validity of the SM:
If the $B$-triangle ``does not close'', New Physics must be the
reason. One should mention that the all present measurements of the
CKM unitarity are in agreement with the SM, although a few ``tensions''
have become noticeable (for details see, e.g. ~\cite{CKMfitter}).

Another way of gaining sensitivity at the SFF to NP is the study of
rare decays of $B$ mesons and $\tau$ leptons. Some examples should
illustrate this point, such as $B \rightarrow X_d \nu \bar{\nu}$ or
$\tau \rightarrow \mu \gamma$. These decays involve (several) neutral
particles in the final state and can therefore only be measured at
SFF's. Such decays are highly suppressed (in the $B$ case) or even
completely forbidden (the $\tau$ case). With a SFF, branching fractions
down to several $10^{-9}$ can be probed. More details on rare
decays and their potential to search for New Physics can be found in
the proposals for the SFFs~\cite{superkekb,superb}.

\begin{figure}
  \includegraphics[height=.17\textheight]{Loop_Diagram}
  \caption{\label{fig:penguin}Example of a SM process at the quantum
  loop level (``penguin diagram'', left) with ``New Physics'' contributing
  (right) to $B$-meson decay amplitudes.}
\end{figure}

One of the flagship measurements at the SFF is the precise
determination of the time-dependent CP violating asymmetries. These
asymmetries are measured by observing the decay rate $\Gamma$, as
function of time, of a $B^0$ meson decaying into a specific CP eigenstate
$f_{CP}$, as compared to the same final state coming from the
${\bar{B}}^0$. The CP violating time-dependent asymmetry is defined as:
$$
{\cal A}(f_{CP},\Delta t) =\frac{\Gamma({\bar{B}}^0 \rightarrow f_{CP};\Delta
t)-\Gamma(B^0 \rightarrow f_{CP};\Delta t)}{\Gamma({\bar{B}}^0 \rightarrow f_{CP};\Delta
t)+\Gamma(B^0 \rightarrow f_{CP};\Delta t)},
$$
where $\Gamma$ is the decay rate of the $B^0({\bar{B}}^0)$ into the CP
eigenstate $f_{CP}$ (``CP side'') within some time interval $\Delta
t$, to be explained below. In order to determine which of the
two flavors
($B^0$ or ${\bar{B}}^0$) has decayed into the common final state $f_{CP}$,
the other $B$-decay (``tag side'') is analyzed for a specific flavor,
using, e.g., semi-leptonic decays ($B^0 \rightarrow X l^- \bar{\nu}, {\bar{B}}^0 \rightarrow
X l^+ \nu$). The
charge of the lepton uniquely identifies (``tags'') the flavor of the
$B$ meson, where a positive (negative) lepton signals a
${\bar{B}}^0$($B^0$). Since the two $B$ mesons are produced in an entangled
state by virtue of the quantum numbers of the $\Upsilon (4S)$, the
tag side uniquely\footnote{strictly speaking, the quantum entanglement
is broken when the first of the $B$-mesons decays. From then on the
other $B$-meson is freely oscillating between $B^0$ and ${\bar{B}}^0$, but with
a known time dependence.} determines the flavor of the $B$-meson that
decayed into the CP
eigenstate. The time difference $\Delta t$ is given by the difference
in decay times between the tag side and the CP side. Note that $\Delta
t$ can be positive or negative.

Depending on the CP eigenstate chosen, any of the three angles $\phi_1
(\beta), \phi_2(\alpha)$ or $\phi_3 (\gamma)$ of the $B$-triangle can
be measured. Chosing,
e.g., the final state $J/\psi K^0$, the angle $\phi_1$ (or $\beta$) is
determined. A recent measurement~\cite{jpsiks} of the time-dependent CP
asymmetry for the $J/\psi K^0$ channel is shown in
fig.~\ref{fig:jpsiksacp}. Here, the asymmetry is given for both odd
($J/\psi K_S$) and even ($J/\psi K_L$) CP final states. Since the $\Delta t$
distributions for the CP eigenstates are different, CP violation
is established. Recent summaries on the measurements of the angles
$\phi_1, \phi_2, \phi_3$ ($\alpha, \beta, \gamma$) have been presented at
this conference~\cite{angle_summ}.

\begin{figure}
  \includegraphics[height=.3\textheight]{jpsiksacp}
  \caption{Recent measurements of the time-dependent CP asymmetry
  using the final states $J/\psi K_S$ and $J/\psi K_L$ (see~\cite{jpsiks}).}

\label{fig:jpsiksacp}
\end{figure}

\section{Machine Upgrade: SuperKEKB}

While no significant deviations from the SM predictions have
been observed so far, there are some tantalizing hints
for possible New Physics in $B$ decays (see,
e.g.,~\cite{hints}). Clarification can
only come with a new generation of SFF's, which should aim at
integrated luminosities in excess of 50 /ab (the present world record
KEKB accelerator is about to accumulate 1 /ab). Such large integrated
luminosities require an equally large increase of the instantaneous
luminosity ${\cal L}$ which, in its simplified form, is given by
$$
{\cal L}=\frac{N_1N_2f}{4\pi\sigma_x\sigma_y}.
$$
Here, $N_i$ are the numbers of particles in each of the two colliding
bunches, $f$ is the bunch collision frequency, and $\sigma_{x,y}$ are
the transverse dimensions of the colliding bunches.

At KEK, an extremely strong
accelerator research program is focussing on an asymmetric $e^+e^-$
collider with instantaneous luminosities in excess of 8 $\times
10^{35}$ cm$^{-2}$ s$^{-1}$ (which is about 40 times the present
world record luminosity of 2.11 $\times 10^{34}$ cm$^{-2}$ s$^{-1}$,
reached in May 2009 with KEKB). The new machine, called ``SuperKEKB'',
is an upgrade of the present KEKB machine and should start producing
luminosity by the year 2013. According to the current plan the
KEKB accelerator should stop running by the end of 2009, so that the
construction work for SuperKEKB can start in 2010. The expected
luminosity development of the SuperKEKB machine is shown in
fig.~\ref{fig:lumidev}.
\begin{figure}
  \includegraphics[height=.35\textheight]{Luminosity}
  \caption{Expected development of the luminosity for the SuperKEKB
  machine.}
\label{fig:lumidev}
\end{figure}

Two options for the SuperKEKB machine have been discussed: The
``high-current'' (HC) option and the ``nano-beam'' (NB)
option (see Table~\ref{tab:machine}). Initially, the HC option was
favored. It was characterized by
a mild decrease of the beta function in the low energy ring (LER), and
a dramatic increase of the beam current in the LER and the high energy
ring (HER). In addition, an ingenious scheme was developed to counteract
the luminosity loss due to the finite angle under which the two
colliding bunches cross each other, i.e. the crab crossing scheme. In
this scheme special crab cavities before and after the interaction
region rotate the bunches, so that they collide head-on, instead of at
an angle. The crab crossing scheme was verified with the running KEKB
accelerator, increasing the instantaneous luminosity by about 20
percent to a new world record. However, it became apparent that a new
effect would make the HC option difficult: For the crab crossing to
work efficiently the bunch length must be smaller than the $\beta$
function at the interaction point (IP), i.e. $\sigma_z
< \beta^*_y$. This condition, imposed to avoid the so-called hour-glass
effect, creates no problems at low beam currents. However, at
the large beam currents of the HC option (see
Table~\ref{tab:machine}) coherent synchrotron radiation becomes an
issue which has the effect of lengthening the bunch, thus running into
the hour-glass effect. The unavoidable bunch lengthening leads to a
decrease of the obtainable maximum luminosity to slightly over 50
$\times 10^{34}$. Failing the goal of ${\cal L} \ge 80 \times
10^{34}$, this option is disfavored now.

Following the ideas of the final focus system envisaged for a
future linear collider, transferred to a circular machine as
laid down in~\cite{superb}, the NB option has become the
baseline~\cite{machine_par}. Here, the key parameters are a factor of
two increase in the
beam currents with respect to the present KEKB, but strongly reduced
$\beta$
functions at IP (factor of about 20 smaller than in the HC option). In
addition, low emittance beams are necessary to achieve the desired
small (less than 100 nano meters) transverse bunch sizes. Such a beam
size has been achieved at the ATF damping ring facility at KEK. While
a low emittance electron beam can be prepared by virtue of a carefully
designed injection system, the transverse phase-space of the positrons
needs to be cooled in a special new damping ring.

\begin{table}
\begin{tabular}{l|c|c|c|c|}
\hline
  & \tablehead{1}{c}{b}{KEKB  \\ Design}
  & \tablehead{1}{c}{b}{KEKB achieved \\ (): with crab}
  & \tablehead{1}{c}{b}{SuperKEKB \\ High-Current Option}
  & \tablehead{1}{c}{b}{SuperKEKB \\ Nano-Beam Option} \\
\hline
$\beta^*_y$ (mm) (LER/HER) & 10 / 10 & $\begin{array}{c} 6.5 / 5.9 \\
(5.9 / 5.9) \end{array}$ & 3 / 6  & 0.26 / 0.26 \\

\hline
$\epsilon_x$ (nm) & 18 / 18 & 18 / 24 & 24 / 18  & 2.8 / 2.0 \\
\hline
$\sigma_y$ ($\mu$m) & 1.9 & 1.9 (0.94) & 0.85 / 0.73  &  0.073 / 0.097 \\
\hline
$\xi_y$ & 0.052 & $\begin{array}{c}0.108 / 0.057 \\ (0.129 /
0.090) \end{array}$ & 0.3 / 0.51  & 0.079 / 0.079 \\

\hline
$\sigma_z$ ($\mu$m) & 4 & ~7 & 5 (LER) / 3 (HER)  & 5 \\
\hline
$I_{\rm beam}$(A) & 2.6/1.1& $\begin{array}{c}1.66 / 1.34 \\ (1.64 /
1.19) \end{array}$ & 9.1 / 4.1 & 3.84 / 2.21 \\

\hline
$N_{\rm bunches}$ & 5000 & 1388 (1585) & 5000 & 2252 \\
\hline
$\begin{array}{l}{\rm Luminosity} \\ (10^{34} {\rm cm}^{-2} {\rm
s}^{-1}) \end{array} $& 1 & 1.76 (2.11) & 53 & 80 \\

\hline
\end{tabular}
\caption{Comparison of parameters for the two high luminosity options
of the SuperKEKB machine, compared to the presently running KEKB
factory. The high current option for SuperKEKB has been discarded
recently (see text) and the nano-beam option is now the baseline.}

\label{tab:machine}
\end{table}

Table~\ref{tab:machine} gives an overview of the machine parameters
for the presently running KEKB accelerator, and for both the HC and NB
options (status of July 2009). Note that the presently achieved
luminosity of KEKB is a factor of two larger than
the original design value. Part of
it is due to the new crab crossing scheme, which is being successfully
tested with the running machine. It should be noted, however, that the
luminosity increase predicted in simulation for the crab crossing
scheme has not yet
been reached in the real machine. The discrepancy is attributed to
so-called ``machine errors'' which may result from a finite precision in
aligning the focussing quadrupoles along the ring. Small (unknown)
deviations from the ideal position can lead to a coupling of the
horizontal and vertical betatron oscillations and therefore to an
increase in the vertical beam size at the IP, reducing the
luminosity. During the last winter shutdown some errors of this kind
have been diagnosed and corrected for by installing a set of so-called
``skew sextupoles'' around the ring. As a consequence, the luminosity
could be immediately increased by more than 15 percent, leading to
new world records for the instantaneous luminosity (see
table~\ref{tab:machine}).

Preliminary parameters of both collision schemes for the SuperKEKB
machine are given in the table~\cite{machine_par}. Intense
optimization studies for
the new machine are going on, concentrating now mainly on the NB
option. One should note, however, that the parameters given in the
table have not yet been consolidated. In addition, the so-called crab
waist scheme, as proposed in~\cite{superb}, has not yet been included
in the present design. While the crab waist scheme is expected to
improve the instantaneous luminosity only by a few percent, it is very
effective in reducing the coupling between vertical and horizontal
betatron oscillations and is therefore quite useful to stabilize the
machine operation.

Some of the virtues of the NB option are quite evident: There is no
need to have short bunches (the condition $\sigma_z < \beta^*_y$ vital
for the HC option to avoid the hour glass effect does not apply here),
and the beam-beam parameter $\xi_y$ can be relaxed substantially (it
should be noted that the large $\xi_y$ parameter required for the HC
option has never been reached so far). Also
the synchrotron radiation is less of an issue due to the smaller beam
currents. This is certainly good news for the detectors close to the
beampipe (e.g. the silicon detectors). On the other hand, a new source
of background may arise from the Touschek effect~\cite{touschek}, a
kind of intra-beam Coulomb scattering, which couples betatron and
synchrotron oscillations. Due to this effect particles may be sent to
off-momentum orbits, causing particle loss and a severe decrease of the
beam lifetime. To counteract the shorter lifetimes (order of a few
minutes) an elaborate scheme is envisaged for SuperKEKB which keeps
the bunch currents almost constant by permanent re-injection at a
rate of about 50 Hertz. With the NB option another important physics
parameter may be chosen more favorably, i.e. the radius of the
beampipe, which may be as small as roughly 1 cm, whereas in the HC
option 1.5 cm had to be chosen. It is expected that the final design
for the SuperKEKB machine in its nano-beam version will be completed
in the fall of 2009.

\section{Detector Upgrade: BELLE-II}

Parallel to the machine studies, a Letter of Intent~\cite{belleloi}
for the Belle upgrade
was issued in the year 2004, with a recent update from
2008~\cite{sbelle}.  The main (re)design goals are to cope with the
much higher physics rates and the much larger backgrounds to be
expected, as well as improving the overall physics performance.

A comparison of Belle and its upgraded version (``Belle-II'') is shown
in a sideview in fig.~\ref{belleii}. We will first give a short
overview of the different detector components and then pick out two
specific systems which will be very important for the precise
measurement of the CP quantities. As a general side condition, the
performance of the new Belle-II detector should be as good or better
than Belle. This is a non-trivial requirement in view of the
anticipated very high background at the SuperKEKB machine, which is
estimated to be roughly an order of magnitude larger as compared to
the present KEKB machine (for details, see~\cite{sbelle}).
\begin{figure}
  \includegraphics[height=0.47\textheight]{Belle_BelleII}
  \caption{Upgrade scenario from the Belle detector (lower half) to the
  Belle-II detector (upper half).}
\label{belleii}
\end{figure}

The tracking system in the ``old'' Belle detector consisted of 4 layer of
Si strip vertex detectors (SVD), followed by a Central Drift Chamber
(CDC). Due to the largely increased background, strip detectors are no
longer an option for the innermost Si layers. Instead a two-layer
pixel detector (PXD) for the innermost Si layers is planned (see
below). The SVD will be replaced entirely as well as the CDC: due to
the harsh backgrounds the inner radius of the CDC has to be moved out
and the two outer layers of the new SVD will cover the gap. The
momentum resolution of charged particles will be improved by extending
the CDC to a larger radius. Since the magnet and the barrel part of
the electromagnetic calorimeter (ECL) will not be changed, the
particle identification system has to be replaced by a thinner,
low material budget detector (see below). In the endcap parts of the ECL
the present CsI(Tl) crystals will be replaced by pure CsI that
provide faster signals, and the forward part of the KLM ($K_L$ and
muon detector) in the iron flux return yoke will be instrumented by
scintillator strips with SiPM readout, replacing the present RPCs.

To extend the physics reach for Belle-II, the $K/\pi$ separation
ability will be improved by a new particle identification (PID)
system. Two types of detectors are being proposed, one is a
time-of-propagation (TOP) counter for the barrel region (Barrel PID),
and the other a proximity-focussing Cherenkov ring imaging counter
with aerogel radiators (ARICH) for the endcap region (Endcap PID).
\begin{figure}
  \includegraphics[width=\textwidth]{PID}
  \caption{Comparison of the PID system performance of Belle (left) and
  Belle-II (right) for the decay $B^0 \rightarrow \rho^0 \gamma$. The abscissa is
  $\Delta E= E_{\rm beam}-E_{B^0}$ (in GeV) in the
  $\Upsilon(4S)$ center-of-mass system (see text).}
\label{fig:pid}
\end{figure}
The present time-of-flight and aerogel Cherenkov counters in the
barrel region of Belle are replaced by a TOP counter made from
quartz radiator bars (thickness 2 cm), in which the time of
propagation of Cherenkov photons is measured, which are internally
reflected and focussed onto micro-channel plate
(MCP) PMTs at the end surfaces of the quartz bars. The MCP-PMTs have
an excellent time resolution of about 50 ps, so that the difference in
arrival times of the Cherenkov photons, radiated by  pions or
kaons, can be determined. This time difference results from the
different Cherenkov emission angles (for equal particle momenta) and
consequently different path lengths for the photons which undergo multiple
internal reflections on their way to the ends of the quartz bars.

The ARICH counters are located in front of the endcap ECL, where the
space is quite limited. For this reason a proximity focussing scheme
is envisaged for the Cherenkov photons, with an expansion thickness of
only 20 cm. In the present design three layers of silica aerogel, each
10 mm thick, are foreseen, with refractive indices varying
between 1.045 and 1.055, so that the photons emitted from the three
regions produce overlapping images on the photon detector surface. Since
the photon detectors have to work in a strong magnetic field of 1.5
Tesla, candidates are a hybrid avalanche photon detector (HAPD) or a
MCP-PMT.

An example of the expected performance of the new PID system relative
to the present Belle ACC system is shown in
fig.~\ref{fig:pid}. Here, the decay $B^0 \rightarrow \rho^0 \gamma$ is studied at
the $\Upsilon(4S)$ resonance and the distribution of the quantity
$\Delta E=E_{\rm beam}-E_{B^0}$ is chosen for reference, where $E_{\rm beam}$ is
the beam energy in the $\Upsilon(4S)$ center-of-mass system and $E_{B^0}$ is
the energy of the $B$ meson, reconstructed from the two pions forming
the $\rho$ and the photon.  An overwhelming
background from $B \rightarrow K^* \gamma$ is expected due to the CKM
couplings (about a factor 40). While for the case of Belle (left side)
the background from misidentified kaons is very large, the improved PID in
Belle-II (right) largely reduces the background and shows a much clearer
signal.
\begin{figure}
  \includegraphics[width=0.27\textwidth]{depfet_schematics}
  \caption{Principle of operation of the DEPFET pixel sensor.}
  \label{fig:depfet}
\end{figure}

A second example for the improved instrumentation at Belle-II is the
new two-layer Si pixel detector, located closest to the beampipe, for
precise vertexing. As mentioned above, the expected strong increase of
background relative to the present KEKB machine, proportional to
$1/r^2$ where $r$ is the radial distance from the
interaction region, creates large occupancies within the strip
detectors, making an efficient reconstruction of tracks and vertices
impossible. The solution is to use a pixel detector which
intrinsically provides three-dimensional space points. However,
pixel segmentation is not the only requirement to be fulfilled at
Belle-II. Due
to the small (transverse) momenta of the $B$ decay products, the
momentum and vertex resolutions are dominated by multiple scattering,
so a pixel sensor should have a very low material budget. Furthermore, the
harsh radiation environment at the SuperKEKB factory requires
low-noise and radiation-hard technology, and the limited space between
the beampipe and the strip detector requires a low power
consumer. Finally, a pixel detector should be operational from the
start of the Belle-II running, which is envisaged for the year 2013.

It turns out that there is only one mature technology at present
fulfilling all these requirements: The DEPFET pixel
sensor~\cite{depfet}, invented and developed at the
Max-Planck-Institute for Physics (MPI) in Munich, Germany. The name
DEPFET stands for $DEP$leted $F$ield $E$ffect $T$ransistor and combines
detection and amplification. The DEPFET principle of operation is
shown in fig.~\ref{fig:depfet}. A MOS field effect transistor is
integrated onto a fully depleted silicon substrate forming the
detector. By
means of an additional $n$-implant underneath the transistor channel a
potential minimum for electrons is created. This potential well can be
considered as an internal gate of the transistor. A particle entering
the detector creates electron-hole pairs in the fully depleted silicon
substrate. The electrons (``signal'') are collected and stored in the
internal gate. The electron charges change the potential of the
internal gate, resulting in a modulation of the transistor channel
current. After readout of the channel current, the signal electrons
are removed by a positive voltage at the clear contact. Another read
cycle then establishes the baseline which is subtracted from the
current readout before the clear. The whole readout cycle
(read-clear-read) takes about 80 ns. Because
of the small capacitance of the internal gate, the noise in the DEPFET
is very low, about 100-200 electrons are expected when operating at
Belle-II.

For a real detector a matrix of DEPFET pixels must be
constructed. Such matrices, using the SOI technique, where the sensor
wafer is bonded to a ``handling'' wafer, have already been
built and subjected to extensive beam tests. Matrices up to 512
$\times$ 512 pixels have been obtained, with pixel sizes of about 17
$\times$ 13 $\mu$m$^2$. For the DEPFET pixel detector (``PXD'') at
Belle-II typical pixel sizes would be around 50 $\times$ 50 $\mu$m$^2$.
An additional asset of the DEPFET technology is the fact that the
detector can be thinned down to about 50 $\mu$m thickness. The thinning
procedure of the handling wafer has already been demonstrated at the
MPI semiconductor laboratory.  With such a thin sensor still a signal
to noise ratio of about 40:1 can be reached. Finally, the radiation
hardness of the DEPFET matrix has been tested and found to be
satisfactory for a few years of running at SuperKEKB with full
luminosity. Further details of the DEPFET project at Belle-II can be
found elsewhere~\cite{depfetpage}.

\section{Conclusions}

Complementary to the LHC, a new generation of Super Flavor Factories
is being planned, probing the Standard Model at energy scales
beyond tens of TeV. Flavor physics continues to be the key to the
puzzle of the observed matter-antimatter asymmetry in the universe,
intimately connected to CP violation. The SuperKEKB project, together
with a substantially upgraded detector, Belle-II, will contribute to
this fascinating chapter of particle physics, taking first data  in
the year 2013, with the prospect of accumulating 50 times the presently
available data by the year 2020.

\begin{theacknowledgments}

The author would like to acknowledge the extremely stimulating and
pleasant atmosphere created by the organizers, especially Marvin
Marshak, of the CIPANP 2009 Conference. He also wishes to thank Tom
Browder, Yoshihiro Funakoshi, Katsunobu Oide, and Yoshi Sakai for
their careful reading of the manuscript and their valuable advice.

\end{theacknowledgments}

\bibliographystyle{aipproc}

\IfFileExists{\jobname.bbl}{}
 {\typeout{}
  \typeout{******************************************}
  \typeout{** Please run "bibtex \jobname" to optain}
  \typeout{** the bibliography and then re-run LaTeX}
  \typeout{** twice to fix the references!}
  \typeout{******************************************}
  \typeout{}
 }

\input{cmkbibl}

\end{document}

\endinput

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

This paper is motivated by the design and implementation of an
On-Demand Multimodal Transit System (ODMTS) for the city of Atlanta. The share of
public transit in Atlanta (about 2--3\%) is very low compared to other
American cities (e.g., about 15\% in Boston) and Atlanta is also the
8th most congested city in the world. There is thus a strong need for
a modern transit systems that leverages the train and bus
infrastructure of the city and complements it with innovative mobility
concepts.

This paper considers the design of an ODMTS for Atlanta that combines
a network of trains and buses with on-demand multimodal shuttles that
act as feeders to/from the bus/rail network and serve local
demand. ODMTS address the first/last mile problem that plagues transit
systems, while mitigating congestion on high-density corridors and
leveraging economy of scale. ODMTS and their design challenge was
introduced in \cite{MaheoEtAl2019-BendersDecompositionDesign}, which
also presents an overview of related work. The main contribution of this
paper is to generalizes prior work by
including three additional elements that are critical for ODMTS in
large cities such as Atlanta.  First, different frequencies are
allowed throughout the network.  Second, additional modes of transit
(e.g., rail) are included.  Third, a limit on the number of transfers
per passenger is introduced.  Adding a constraint to limit the number
of transfers has a significant negative impact on existing Benders
decomposition approaches as it introduces non-convexity in the
subproblem. Instead, this paper enforces the limit through
transfer-expanded graphs, i.e., layered graphs in which each layer
corresponds to a certain number of transfers.  A real-world case study
is presented for which the generalized ODMTS design problem is solved
for the city of Atlanta.  The results demonstrate that exploiting the
problem structure through transfer-expanded graphs results in
significant computational improvements.

\section{The Generalized ODMTS Design Problem}

This section presents the generalized ODMTS design problem that
enhances the model from
\cite{MaheoEtAl2019-BendersDecompositionDesign} along several
dimensions: The choice of bus frequencies, additional transportation modes and, most importantly, a constraint on the
number of transfers. The Benders decomposition approach in
\cite{MaheoEtAl2019-BendersDecompositionDesign} exploits a natural
decomposition of the ODMTS design problem.  The network design is
determined by the \emph{master problem}, while the routing of the
passengers for a given design is determined by the \emph{subproblem}.  A
major benefit of this decomposition is that the subproblem can be
solved for each trip independently. The same decomposition is used in this paper.

\subsection{The Master Problem For Network Design}

Consider a directed multigraph $G=(V,A)$, with vertices $V =
\{1,\hdots,n\}$ and arc set $A$. Let $F$ be the set of possible
frequencies, i.e., the total number of vehicles during the time
horizon, let $M$ be the set of possible transportation modes,
which may include shuttles, and let $K$ be the total number of arcs that
each passenger may travel.  By definition, $K$ is equal to the maximum
number of transfers plus one. In the multigraph $G$, each arc $a \in
A$ is uniquely defined by the quadruple $a = (i, j, m, f) \in V \times
V \times M \times F$, $i \neq j$.  Using arc $a$ means traveling from
$i$ to $j$ with mode $m$, which departs with frequency $f$.  For a
given arc $a\in A$, these elements are referred to as $i(a)$, $j(a)$,
$m(a)$, and $f(a)$, respectively.

Designing a generalized ODMTS amounts to deciding which arcs $a \in A$ are made available to passengers.
Let the binary variable $z_a \in \mathbb{B}$ be equal to one if arc $a$ is made available, and zero otherwise.
The cost of enabling arc $a$ is given by the parameter $\beta_a$.
It is assumed that $\beta_a \ge 0$ for all $a \in A$.

For a given design, a cost is incurred due to passengers traveling trough the network.
This cost $\Phi(z)$ is a function of the values of the $z$-variables that define the design.
The value of $\Phi(z)$ can be found by solving the subproblem, which is discussed in Section~\ref{sec:subproblem}.
If the subproblem is not feasible, then $\Phi(z) = \infty$.

A formulation for the master problem is presented in Figure~\ref{fig:formulation}.
For convenience, $\delta^+(i)$ is defined as the set of all arcs going out of $i \in V$.
Similarly, the set $\delta^+(i,m)$ is defined as the set of all arcs with mode $m\in M$ going out of $i \in V$.
The sets $\delta^-(i)$ and $\delta^-(i,m)$ are defined analogously for the incoming arcs.

\newcounter{masterFormulation}
\setcounter{masterFormulation}{\value{equation}}
\addtocounter{masterFormulation}{1}
\begin{figure}[t]
	\makebox[\textwidth][c]{
		\fbox{\begin{minipage}{\textwidth}
				\begin{alignat}{3}
				~ & \min \quad \mathrlap{\sum_{a \in A} \beta_a z_a + \Phi(z),} \label{eq:objectiveOveral} \tag{\arabic{masterFormulation}a} &&&\\
				~ & \textrm{s.t.} & \sum_{a \in \delta^+(i,m)} f(a) z_a - \sum_{a \in \delta^-(i,m)} f(a) z_a & = 0 & \quad & \forall i \in V, m \in M, \label{eq:masterFlowConservation} \tag{\arabic{masterFormulation}b}\\
				~ & ~ & \sum_{f \in F \vert (i,j,m,f) \in A} z_{(i,j,m,f)} & \le 1 && \forall i \in V, j \in V, m \in M, \label{eq:masterOneFrequency} \tag{\arabic{masterFormulation}c}\\
				~ & ~ & z_a & \in \mathbb{B} && \forall a \in A. \label{eq:masterVariables} \tag{\arabic{masterFormulation}d}
				\end{alignat}
	\end{minipage}}}
	\caption{Formulation for the generalized ODMTS design problem.}
	\label{fig:formulation}
	\vspace{-0.5cm}
\end{figure}
\addtocounter{equation}{1}

Objective~\eqref{eq:objectiveOveral} minimizes the cost of the design
plus the cost of routing the passengers through the network.
Constraints~\eqref{eq:masterFlowConservation} ensure that the
frequencies for each mode are balanced at each vertex.  For example,
if three buses arrive during the time horizon, then three buses
should also depart.  Constraints~\eqref{eq:masterOneFrequency} enforce that
only one frequency can be selected for a given connection and a given
mode.  Equations~\eqref{eq:masterVariables} state the integrality
requirements.

\subsection{The Subproblem: Routing Passengers Through the Network}
\label{sec:subproblem}

For a given design, the passenger trips are routed through the network
at minimum cost.  Let $T$ be the set of all passenger trips, and let
each trip $r \in T$ be defined by an origin $o(r)$, a destination
$d(r)$, and a number of passengers $p(r)$.  If trip $r \in T$ is
routed through arc $a\in A$, then a cost of $\gamma_a^r$ is incurred.
The total cost of routing all passenger trips, $\Phi(z)$, is the sum
over the costs per trip.  It is assumed that $\gamma_a^r > 0$ for
every arc $a \in A$ and trip $r \in T$, such that the optimal route is
a simple path from $o(r)$ to $d(r)$.

Solving the subproblem amounts to solving a shortest path problem from
$o(r)$ to $d(r)$ for each trip $r \in T$, with the additional
restriction that the number of arcs in the path is at most $K$.  This
problem is known as the cardinality-constrained shortest path problem
(CSP) \cite{DahlRealfsen2000-CardinalityConstrainedShortest}.  Note
that the cardinality constraint follows from the limit on the number
of transfers.  Without this limit, the subproblem is an
(unconstrained) shortest path problem (SP), as is the case in
\cite{MaheoEtAl2019-BendersDecompositionDesign}.

It is well-known that SP possesses total unimodularity and can
be solved by linear programming (LP).  Adding an additional
constraint, however, typically destroys this structure
\cite{AnejaNair1978-ConstrainedShortestPath}.
This is indeed the case when a cardinality constraint is added to the
subproblem formulation in
\cite{MaheoEtAl2019-BendersDecompositionDesign}.  As a result, the
cost function $\Phi(z)$ would change from convex to non-convex, which
negatively impacts Benders decomposition approaches (see
Section~\ref{sec:Benders}).

{\em To remedy this limitation, this paper presents a new formulation for
the subproblem that enforces the transfer limit without destroying
total unimodularity.}  This formulation uses transfer-expanded graphs,
i.e., layered graphs with a each layer for each number of transfers. Transfer-expanded graphs encode the transit constraints
directly, making it possible to use shortest-path algorithms.

\subsection{Transfer-Expanded Graphs}
\label{sec:transfergraphs}

Transfer-expanded graphs shares some similarities with time-expanded
networks, where each vertex has multiple copies for different periods
of time. This is the case, for example, for modern algorithms for evacuation
planning and scheduling
\cite{HasanHentenryck2017-ColumnGenerationAlgorithm,PillacEtAl2015-ColumnGenerationApproach,PillacEtAl2016-ConflictBasedPath}. Reference
\cite{MocciaEtAl2009-ColumnGenerationHeuristic} also uses a layered
network to solve the dynamic generalized assignment problem.  As a
result, some of the side-constraints do not need to be handled
explicitly. See \cite{BolandEtAl2017-ContinuousTimeService} for a
recent literature review on time-expanded graphs.

Let $\bar{G}^r = (\bar{V}^r, \bar{A}^r)$ be the
transfer-expanded graph for a given trip $r\in T$.  This graph
contains multiple copies of the original arcs and vertices, organized
in $K+1$ layers.  It is assumed that $K \ge 2$, as the subproblem is
trivial for $K=1$.  A vertex $\bar{v} = (i,k) \in \bar{V}^r$ in the
transfer-expanded graph is defined by a vertex $i \in V$ in the
original graph and by a layer $k \in \{1,\hdots,K+1\}$.  Similarly,
the definition of an arc is extended to $\bar{a} = (a, k, l)$, in
which $a \in A$ is the original arc, $k \in \{1,\hdots,K\}$ is the
layer of the starting vertex of $\bar{a}$ and $l \in \{2,\hdots,K+1\}$
is the layer of the ending vertex.

The transfer-expanded graph is constructed as follows.  For
convenience, Figure~\ref{fig:exampleexpanded} provides an example for
$K=3$.  First, the vertex set $\bar{V}^r$ is defined.  For the origin
and the destination of the trip, introduce the vertices $(o(r),1)$ and
$(d(r),K+1)$.  For the other vertices $i \in V\backslash\{o(r),d(r)\}$
of the original graph, add the copies $(i,k)$ for $k \in \{2, \hdots,
K\}$ to the transfer-expanded graph.  The arc set $\bar{A}^r$ is
constructed based on the arcs of the original graph, as follows:

\begin{enumerate}
\item For each arc starting in the origin, i.e., $a\in
          \delta^+(o(r))$, add the arc $(a,1,2)$ if $j(a) \neq d(r)$,
          or the arc $(a,1,K+1)$ if $j(a) = d(r)$.

\item For each arc not adjacent to the origin or the destination,
  i.e., $a \in A$ and $i(a),j(a) \notin \{o(r),d(r)\}$, add the arcs
  $(a, k, k+1)$ for all $k \in \{2, \hdots, K-1\}$.

\item For each arc ending in the destination that does not start in
  the origin, i.e., $a \in \delta^-(d(r))$, $i(a) \neq o(r)$, add the
  arcs $(a, k, K+1)$ for all $k \in \{2, \hdots, K\}$.
\end{enumerate}

\begin{figure}[!t]
	\centering
	\includegraphics[scale=0.8]{expanded.png}
	\caption{A transfer-expanded graph for $K=3$, $\lvert M
          \rvert = 1$, $\lvert F \rvert = 1$, for a complete graph as
          the original graph. The dotted arcs are removed when the
          triangle inequality holds.}
	\label{fig:exampleexpanded}
\end{figure}

\noindent
By construction, it follows that solving CSP on the original graph is
equivalent to solving SP on the transfer-expanded
graph. Figure~\ref{fig:expandedsubproblem} formulates the subproblem
as a collection of SPs on transfer-expanded graphs.  Let
$y_{\bar{a}}^r \in \mathbb{B}$ be the flow on arc $\bar{a} \in
\bar{A}^r$ of trip $r \in T$.  For convenience, define
$\bar{\delta}^+_r(\bar{v})$ to be the set of all arcs in $\bar{A}^r$
coming out of $\bar{v} \in \bar{V}^r$.  Similarly, let
$\bar{\delta}^-_r(\bar{v})$ be the set of incoming arcs.

\newcounter{subproblemFormulation}
\setcounter{subproblemFormulation}{\value{equation}}
\addtocounter{subproblemFormulation}{1}
\begin{figure}[t]
	\makebox[\textwidth][c]{
		\fbox{\begin{minipage}{\textwidth}
				\begin{alignat}{3}
				\Phi(z) = & \min \quad \mathrlap{\sum_{r \in T} \sum_{\bar{a} = (a,k,l) \in \bar{A}^r} \gamma_a^r y_{\bar{a}}^r} \label{eq:extsubObjective} \tag{\arabic{subproblemFormulation}a}\\
				~ & \textrm{s.t.} & y_{\bar{a}}^r & \le z_a & \qquad\qquad & \forall r\in T, \bar{a} = (a,k,l) \in \bar{A}^r, \label{eq:extsubArcCapacity} \tag{\arabic{subproblemFormulation}b}\\
				~ & ~ & \mathllap{\sum_{\bar{a} \in \bar{\delta}^+_r(\bar{v})} y_{\bar{a}}^r - \sum_{\bar{a} \in \bar{\delta}^-_r(\bar{v})} y_{\bar{a}}^r} & \mathrlap{= \begin{cases} 1 &\textrm{ if } \bar{v} = (o(r), 1) \\ -1 &\textrm{ if } \bar{v} = (d(r), K+1) \\ 0 &\textrm{ else } \end{cases} \quad \forall r \in T, \bar{v} \in \bar{V}^r,} \label{eq:extsubFlowConservation} \tag{\arabic{subproblemFormulation}c}\\
				~ & ~ & y_{\bar{a}}^r & \ge 0 && \forall r \in T, \bar{a} \in \bar{A}^r. \label{eq:extsubVariables} \tag{\arabic{subproblemFormulation}d}
				\end{alignat}
	\end{minipage}}}
	\caption{Formulation for the subproblem on transfer-expanded graphs.}
	\label{fig:expandedsubproblem}
	\vspace{-0.5cm}
\end{figure}
\addtocounter{equation}{1}

Objective~\eqref{eq:extsubObjective} minimizes the cost of all trips.
Constraints~\eqref{eq:extsubArcCapacity} state that passengers can
only use arcs available in the design.
Constraints~\eqref{eq:extsubFlowConservation} enforce flow
conservation, and Equations~\eqref{eq:extsubVariables} define the
variables.  Due to total unimodularity of the SPs, no integrality
conditions are required.

The main advantage of using tranfer-expanded graphs is that the limit
on the number of transfers can be enforced without destroying total
unimodularity.  A potential downside is that the number of variables
and constraints in the subproblem increases with $K$.  In public
transit, however, the number of transfers that passengers are willing
to take, and therefore the value of $K$, is typically very low.
Furthermore, a larger subproblem does not necessarily mean that the
subproblem is more difficult to solve, as algorithms may benefit from
the fact that the transfer-expanded graph is acyclic.  When the
$z$-variables are integers, for example, the acyclic subproblem for
each trip can be solved in linear time through topological sorting
\cite{CormenEtAl2009-IntroductionAlgorithms}.

Finally, it is worth pointing out that if $o(r)$ and $d(r)$ are only
served by shuttles, and shuttles satisfy the triangle inequality, then
some arcs may be removed from the transfer-expanded graph without
sacrificing optimality.  Specifically, using a shuttle on the path
$(o(r),1) \rightarrow (i, 2) \rightarrow (d(r),K+1)$ for $i\in V$ is
always dominated by using a direct shuttle from $(o(r),1)$ to
$(d(r),K+1)$.  It follows that the shuttle arcs between $(i,2)$ and
$(d(r),K+1)$ may be removed for all $i\in V$, as also indicated in
Figure~\ref{fig:exampleexpanded}.  For $K\le 3$, it then follows that
the transfer-expanded graph does not require more edges than the
original graph.

\section{Benders Decomposition}
\label{sec:Benders}

Following \cite{MaheoEtAl2019-BendersDecompositionDesign}, a Benders
decomposition approach is presented for the generalized ODMTS design
problem.  The goal is to solve the master
problem~(\arabic{masterFormulation}), which is complicated by the fact
that $\Phi(z)$ is defined implicitly.  To apply Benders decomposition,
replace $\Phi(z)$ in Objective~\eqref{eq:objectiveOveral} by a new
variable $\theta \in \mathbb{R}$, and add the constraint $\theta \ge
\Phi(z)$.  Note that this does not change the problem, as $\theta =
\Phi(z)$ in any optimal solution. In Benders decomposition, the constraint $\theta \ge \Phi(z)$ is enforced through \emph{Benders cuts}.
For subproblem~(\arabic{subproblemFormulation}), these cuts are
\begin{equation}
\label{eq:extbenderscut}
\theta \ge \Phi(\bar{z}) + \sum_{r\in T} \sum_{a \in A} \sum_{k=1}^K \lambda_a^{rk}(\bar{z}) (z_a - \bar{z}_a),
\end{equation}
with $\lambda_a^{rk}(z)$ the dual values of Constraints~\eqref{eq:extsubArcCapacity} and $\bar{z}$ any feasible solution to the LP relaxation of the master problem \cite{Benders1962-PartitioningProceduresSolving}.
For the case study in this paper, the subproblem is always feasible.
If this assumption is not satisfied, \emph{Benders feasibility cuts}, which are similar to \eqref{eq:extbenderscut}, may also be included \cite{Benders1962-PartitioningProceduresSolving}.

The Benders decomposition approach is implemented in C++ and Gurobi
8.1.1.  The master problem is the main model, and the Benders
cuts~\eqref{eq:extbenderscut} are separated in both the MIP solution
callback (in case the $z$-variables are integer) and in the MIP node
callback (in case the $z$-variables are fractional).  The subproblem
for each trip is also solved with Gurobi, and dual simplex is used to
ensure that the basis remains feasible when the subproblem is solved
for different values of $z$.  To prevent excessive calls to the
subproblem, feasibility heuristics are disabled.  The number of cut
separation rounds in the root node is set to the maximum value to get
the best possible bound.  Finally, the $2\epsilon$-trick is used to
stabilize the master problem
\cite{FischettiEtAl2017-RedesigningBendersDecomposition}.  This
stabilization uses $\epsilon = 0.00001$ and the trivial core point
obtained by assigning $z_a = \frac{1}{4}$ to every bus arc.

Without transfer-expanded graphs, the subproblem is not totally unimodular and
$\Phi(z)$ is not convex (see Section~\ref{sec:subproblem}).  In that
case, Benders decomposition cannot be applied directly.  Instead,
$\theta \ge \Phi(z)$ may be enforced by adding \emph{combinatorial
  Benders cuts} in the MIP solution callback and Benders cuts for the LP relaxation of the
subproblem in both callbacks \cite{CodatoFischetti2006-CombinatorialBendersCuts,LaporteEtAl2002-IntegerLShaped}.
However, it is well-known that relying too much on combinatorial
Benders cuts may result in slow algorithmic progress, and many cuts
may be necessary to find the optimal solution.

\section{Atlanta as a Case Study}

The generalized ODMTS design problem was solved for the city of Atlanta.  In Atlanta, the Metropolitan Atlanta Rapid Transit Authority (MARTA) operates two modes: bus and rail. The case study added on-demand shuttles and the bus
system was redesigned accordingly. More precisely, define the three modes $M=\{S,B,R\}$ for shuttle, bus, and
rail respectively.  Shuttle arcs are introduced to connect from
origins to hubs and from hubs to destinations, as well as to serve
the local demand. The corresponding $z_a$ variables are fixed to one, as
shuttles are always available. Following \cite{MaheoEtAl2019-BendersDecompositionDesign}, the cost of using a
shuttle is a weighted sum of cost and convenience, controlled by the
parameter $\alpha \in [0,1]$.  Let $d_a$ and $t_a$ be the travel
distance and the travel time of arc $a \in A$, respectively.  The
parameter $c^S$ is the cost of using a shuttle per person per unit of
distance.  The cost of traversing arc $a \in A$ for trip $r \in T$ is
then defined as $\gamma_a^r = p(r)\left((1-\alpha) c^S d_a + \alpha
t_a\right)$.

Bus arcs are defined between the potential hub locations and between
each hub and the three nearest rail stations.  The cost of enabling
bus arc $a \in A$ is given by $\beta_a = (1-a)c^Bf(a)d_a$.  That is,
the distance is multiplied by the cost per unit distance and the
number of buses over the time horizon.  The cost of traversing a bus
arc is given by $\gamma_a^r = \alpha \left(t_a + L + \frac{H}{2f(a)}
\right)$.  Here $L$ is the fixed time required for a transfer, $H$ is
the time horizon, and $\frac{H}{2f(a)}$ is the expected waiting time
before the next bus arrives, which depends on the frequency. Rail arcs are defined between all rail stations that are connected by the same rail line. The costs of traversing an arc is defined in the same way as for the buses.
For each rail arc $a \in A$, the variable $z_a$ is fixed to one, which makes the cost of enabling an arc irrelevant.

The case study uses the following data and parameters to create a realistic instance and evaluate the computational benefit of transfer-expanded graphs. It uses passenger trip data provided by MARTA for March 16, 2018, between 6am and 10am. Connecting trips have been chained together to obtain origin and destination pairs. This resulted in 2588 unique trips, with 7167 passengers in total. There are 5563 bus stops and rail stations in total, and their locations were also provided by MARTA. Eleven hubs were selected manually on the map.
For the distances $d_a$, great-circle distances are used. To estimate travel times $t_a$, the distances are divided by a constant speed of 30 mph. The cost parameters are set to $c^S = 5$ and $c^B = 1$. The fixed transfer time is chosen to be five minutes, i.e., $L=5$ minutes, and the time horizon is set to four hours, i.e., $H=240$ minutes. To balance cost and convenience, $\alpha = 0.5$ is used. The rail frequency is assumed to be fixed to six per hour, i.e., $f(a) = 6 \times 4 = 24$, and bus frequencies are determined by the model to be either three per hour or six per hour.
At most two transfers are allowed, i.e., $K=3$.

Figure~\ref{fig:atlanta} presents the result of solving the generalized ODMTS problem using transfer-expanded graphs.
In total, it took 122 seconds to obtain the optimal solution and prove optimality, with a minimum objective value of 131,905. Without transfer-expanded graphs, i.e., when adding combinatorial Benders cuts, it was not possible to obtain an optimal solution in reasonable time. Instead, the evaluation considered a relaxation in which the combinatorial Benders cuts were ignored and only the Benders cuts for the LP relaxation of the subproblem were added. Solving this relaxation to optimality took 3.8 hours. Keep in mind that this relaxation explores routes that may require many transfers. To evaluate the quality of the design obtained by the relaxation, the passengers were routed through the transfer-expanded formulation with the $z$-variables fixed to their values found in the relaxation. The result is presented in Figure~\ref{fig:relaxed} and has an objective value of 131,965. Solving the relaxed problem results in a smaller public transit network because the relaxation does not completely enforce the transfer limit.

In summary, the main benefit of the transfer-expanded formulation is the significant computational benefits it provides in capturing the transfer limit. Without transfer-expanded paths, it can be optimal to fractionally select long paths that do not adhere to this constraints. These longer fractional paths likely play a role in the difference of computational performance.

\begin{figure}[t!]
	\centering
	\subfloat[Based on transfer-expanded graphs.\label{fig:atlanta}]{%
		\centering
		\includegraphics[width=0.49\textwidth]{atlanta.png}
	}
	\hfill
	\subfloat[Based on relaxed problem.\label{fig:relaxed}]{%
		\centering
		\includegraphics[width=0.49\textwidth]{relaxed.png}
	}
	\caption{Network designs for Atlanta showing shuttles (thin lines), rail (thick lines), and buses (arrows, orange/light for low frequency and purple/dark for high frequency).}
	\vspace{-0.5cm}
\end{figure}

\section{Conclusion}

This paper presented a generalization of the ODMTS design problem that
introduces three critical elements in practice: different
frequencies, additional transit modes, and a limit on the number of
transfers.  Transfer-expanded graphs are introduced to handle the
transfer limit without negatively impacting existing Benders
decomposition approaches.  The Atlanta case study demonstrates that
this approach is very effective, as transfer-expanded graphs
significantly improve computational performance. Exploiting the problem
structure through transfer-expanded graphs
opens the door to designing increasingly realistic networks in the
future.  One possible extension is to incorporate the capacity of the
on-demand shuttles.  As capacity of these shuttles is typically small,
expanded networks could also be used to model capacity efficiently.

\section{Acknowledgments}
This research is partly supported by NSF Leap HI proposal NSF-1854684.

\clearpage

\bibliographystyle{splncs04}