%% $Id: infocom05.tex,v 1.10 2004/07/23 20:57:37 gokhale Exp $
\documentclass[twocolumn,11pt]{IEEEtran}

\usepackage{epsf}
\usepackage{psfrag}
\usepackage{subfigure}
\usepackage{epsfig}

\def\tenrm{\fontsize{10}{12}\normalfont\rmfamily\selectfont}
\def\BibTeX{{\rmfamily B\kern-.05em{\scshape i\kern-.025em b}\kern-.08em \TeX}}

\newcommand{\Figdir}{./}
\newcommand{\figdir}[1]{\renewcommand{\Figdir}{#1/}}
\newcommand{\Idrawdir}{./}
\newcommand{\idrawdir}[1]{\renewcommand{\Idrawdir}{#1/}}

\newcommand{\textfigure}[3]{
\begin{figure}[ht]
\begin{center}
\begin{tabular}{|l|} \hline
#3 \hline
\end{tabular}
\caption{#2}
\label{fig:#1}
\end{center}
\vspace*{-0.5cm}
\end{figure}
}

\newcommand{\ls}[1]
{\dimen0=\fontdimen6\the\font
\lineskip=#1\dimen0
\advance\lineskip.5\fontdimen5\the\font
\advance\lineskip-\dimen0
\lineskiplimit=.9\lineskip
\baselineskip=\lineskip
\advance\baselineskip\dimen0
\normallineskip\lineskip
\normallineskiplimit\lineskiplimit
\normalbaselineskip\baselineskip
\ignorespaces
}

\newcommand{\commentout}[1]{}

\newcommand{\ix}{\hspace*{1em}}
\newcommand{\hs}{\hspace*{2em}}
\newcommand{\et}{{\it et. al.\ }}
\newcommand{\ets}{{\it et. al.}'s\ }

 
\newtheorem{theorem}{Theorem}[section] 
\newtheorem{example}{Example}[section] 
\newtheorem{proposition}[theorem]{Proposition} 
\newtheorem{corollary}[theorem]{Corollary} 
\newtheorem{lemma}[theorem]{Lemma} \newtheorem{claim}[theorem]{Claim} 
\newtheorem{fact}[theorem]{Fact} 
\newtheorem{definitionx}[theorem]{Definition} 
\newtheorem{algorithm}[theorem]{Algorithm} 
\newenvironment{definition}{\begin{definitionx}\rm}{\end{definitionx}} 
\newcommand{\qedsymb}{\hfill $\square$} 


\newcommand{\twoup}[6]{   
\begin{figure}[t]  
\renewcommand{\subfigtopskip}{3pt}
\renewcommand{\subfigcapskip}{3pt}
%% \hspace*{-.6cm}
\subfigure[#2]{\includegraphics[angle=-90,width=.45\textwidth]{\Figdir #1}}
%\hspace*{-.6cm}
\subfigure[#4]{\includegraphics[angle=-90,width=.45\textwidth]{\Figdir #3}}
%\vspace*{-.7cm}
\caption{#5}
\label{#6}
\end{figure}
}


\begin{document}
\title{Network Simulation via Hybrid System Modeling: A Time-Stepped Approach}
\author{Amogh Kavimandan$^\dag$, Wonsuck Lee$^\ddag$, Marina Thottan$^\ddag$, 
        Anirudha Gokhale$^\dag$, and Ramesh Viswanathan$^\ddag$\\
        $^\dag$Department of EECS, Vanderbilt University, Nashville, TN 37235\\
        $^\ddag$Bell Laboratories, Lucent Technologies, Murray Hill, NJ 07920}

\markboth{IEEE Infocom 2005}{Kavimandan et. al.}
\maketitle

\begin{abstract}
The ever increasing complexity of networks dramatically increases the
challenges faced by service providers to analyze network behavior and
(re)provision resources to support multiple complex distributed
applications. Robust and accurate simulation tools are pivotal to this
cause. The current state of the art in network modeling and
simulation tools, namely, packet-based discrete event simulations and
fluid models, are limited in that they are either applicable only to
particular phases of the network's dynamics or do not scale well to
capture the behavior of immense or complex networks. The recently
proposed {\it hybrid system model} for data communication networks
shows promise in addressing these challenges.  Based on the hybrid
system modeling paradigm, this paper provides novel contributions to
TCP network modeling and analysis/simulation of large scale
networks. The salient distinguishing features of our approach include
a new finite state machine model for a drop-tail queue, a new model
for fast recovery/fast retransmit mode, a revised sending rate model,
and an embedded time-out mode transition mechanism all of which employ
a time-stepped solution method to solve the hybrid system network
models. Our simulation results are consistent with well-known
simulators, such as {\bf ns-2}.

%Second, we describe how our hybrid systems models can be composed to
%provide multi-layer network simulations that can analyze very large
%topologies -- something that is not feasible with state of the art
%packet-based discrete event simulations. Finally, results of our
%experiments comparing our approach with the well known ns-2 simulator.
\end{abstract}

\begin{keywords}
Data Communication Network, Network Simulation, TCP Traffic Modeling,  
Fluid Modeling, Hybrid Systems.
\end{keywords}

\section{Introduction}
%\PARstart{A}{pplications} and services for next-generation distributed
%systems must be flexible, reusable, robust, and capable of providing
%scalable, low-latency quality of service (QoS) to delay-sensitive
%applications.  In addition, communication software must allow
%bandwidth-sensitive applications to transfer data efficiently over
%high-speed networks.  Robustness, flexibility, and reusability are
%essential to respond rapidly to changing application requirements that
%span an increasingly wide range of media types and access
%patterns~\cite{Clark:90}.  Examples of such applications include {\em
%internet collaboration tools} for {\em telemedicine} and {\em
%scientific applications}, {\em distributed mission training},
%emergency response systems, and tele-immersion systems. Supporting
%these critical applications requires adaptive (re)provisioning of 
%resources from the network layer up to the
%application layer. This paper focuses on the challenges of visualizing
%complex network behavior and our solutions so that communication
%service providers are better equipped to support the QoS needs of
%these applications. 
%\subsection{Shortcomings of Existing Tools of the Trade}
%\label{shortcomings}
\PARstart{C}{ommunication} service providers use a variety of tools
including visualization, analysis and simulations to track the
performance of their networks and (re)provision existing and new
applications. For example, the Network Simulator {\bf ns-2}~\cite{ns2}
is the most widely used packet-based simulation tool to understand
network behavior. Simulators like {\bf ns-2} are based on discrete
event simulation of per-packet behavior in the network. Although these
simulators provide accurate analysis of network behavior, they have a
number of shortcomings when applied to the networks of today and the
future.

It is difficult to study heterogeneous and multilayered networks where
they are made up of several layers starting with the physical
layer on which are overlaid data link layers like ATM and MPLS
followed by the internet protocols like TCP and IP. Existing
techniques to analyze network behavior are typically restricted to
individual layers. For example, simulators like {\bf ns-2} are
particularly useful to analyze IP traffic. Thus, it is tedious and in
many cases infeasible, with existing techniques, to determine the
effect of disruptions at one layer on other layers of the network.

Another challenge includes networks of networks where the ever
expanding size of the internet has given rise to a large collection of
interconnected but independently administered networks that include
both public and private networks. Analyzing the behavior of these
networks requires accurately modeling the structure, policies and
behavior of the networks while also scaling to the large size of the
networks. Current state of the art in network analysis do not provide
these capabilities.

As the scale of the network increases it becomes computationally
infeasible to use sequential packet-based discrete event simulators to
analyze the networks. Moreover, the storage requirements for the trace
data generated by these simulations become cumbersome to manage and
mine.

%\begin{itemize}
%\item {\em Heterogeneous, multilayered networks} -- where networks are made
%      up of several layers starting with the physical layer on which
%      are overlaid data link layers like ATM and MPLS followed by the
%      internet protocols like TCP and IP. Existing techniques to analyze
%      network behavior are typically restricted to individual
%      layers. For example, simulators like {\bf ns-2} are particularly
%      useful to analyze IP traffic. Thus, it is tedious and in many
%      cases infeasible, with existing techniques, to determine the
%      effect of disruptions at one layer on other layers of the
%      network. 

%\item {\em Networks of networks} -- where the ever expanding size of the
%      internet has given rise to a large collection of interconnected
%      but independently administered networks that include both public
%      and private networks. Analyzing the behavior of these networks
%      requires accurately modeling the structure, policies and
%      behavior of the networks while also scaling to the large size of the
%      networks. Current state of the art in network analysis do not
%      provide these capabilities.

%\item {\em Computational complexity and scalability of discrete event simulators} --
%      as the scale of the network increases it becomes computationally
%      infeasible to use sequential packet-based discrete event
%      simulators to analyze the networks. Moreover, the storage
%      requirements for the trace data generated by these simulations
%      become cumbersome to manage and mine.

%\item {\em Impact of topology changes} -- service providers are
%      interested in understanding the behavior of the networks of
%      networks as and when there are changes in topology due to node
%      or link failures or congestion. This information is critical to
%      them to determine how best to reprovision resources so that
%      application service level agreements (SLAs) are honored, while
%      also ensuring that the return on investments (ROI) is kept
%      optimal.

%\end{itemize}

Addressing the challenges outlined above requires a new approach to
analyzing network behavior where the computational requirements must
be kept manageable. 


%while allowing {\em cross-layer} network analysis for both normal and
%disruptive network behavior. This paper provides a novel approach to
%the behavioral analysis of multi-layer, networks of networks using
%hybrid systems modeling and analysis.

Hybrid system models combine the strengths of
continuous and discrete models. A system is modeled as a discrete
event system at a coarse granularity where a discrete state transition
occurs when certain guard conditions are met. Within each state,
however, a system is modeled to evolve according to some continuous
dynamics.

% (both parallel and product form),
% I need to verify this claim - Andy
%
% What does parallel and product mean --- is this asynchronous and
% synchronous parallel composition?  Parallel and product suggests to
% me that product is not supposed to be for concurrent execution of
% two hybrid systems. - Ramesh
% For now I have removed the contentitious phrase :-) - Andy

%Hybrid systems are closed under composition, 
%which allows composition of hybrid models of network layers and the 
%analysis of whole system behavior.

\section{Survey of Work on Network Modeling \& Simulation} % Wonsuck  
Simulation of data communication networks is instrumental to study the
dynamics of the network under various conditions. Simulation studies
are important and relevant to i) designing a new network, ii)
evaluating the existing network, iii) performing stress test of a new
communication protocol, and iv) testing network readiness to support a
new service. Network modeling and simulation paradigms can be
categorized into three mainstream approaches: packet-level models,
flow-level models, and hybrid models.

Packet-level simulation has been the most widely accepted and used
simulation methodology by the network simulation
community. Packet-level simulators are based on event-logic driven
simulation paradigm. For example, {\bf ns-2}~\cite{ns2}, {\bf
SSFNet}~\cite{ssfnet} and {\bf JSim-INET}~\cite{jsim,inet} fall into
this category. Packet-level simulators keep track of all the
individual packets in the network. Therefore, they offer highly
accurate and detailed results; however, the computational complexity
increases rapidly as the size of a network becomes large.

During the past decade fluid-level simulations have been studied
extensively and continue to be researched actively. In the framework
of a fluid models, a chunk of packets is modeled as a single
continuous artifact, namely, a fluid. The dynamics of the network is
then represented by a system of ordinary differential equations, which
describe fluid rate changes inside the network. Most of the research
efforts in this area~\cite{Bonald, Misra99a, Misra99b,Altman} have
focused on modeling TCP protocols~\cite{Jacobson}.  In this framework,
the arrival and the loss of packets are often modeled as stochastic
processes~\cite{Altman, Misra99a}.  The traffic source and queue
dynamics constitute a system of equations which governs the TCP
controlled data transfer over an IP network.

Recently, hybridization of the packet and the flow modeling approaches
have been studied.  Guo \et proposed the Time-stepped Hybrid
Simulation (TSHS)~\cite{Guo00a} method to deal with the scalability
issue faced by traditional packet-level discrete-event simulation
methods. TSHS considers a chunk, grouped from the packets that are in
the same time-step, as a unit entity. The packets in a chunk are
assumed to be evenly spaced within the time-step. In order to identify
any event, every node in the network is checked at every
time-step. TSHS is capable of offering flexible choice in simulation
fidelity based on the simulation abstraction level. With a slightly
different objective, Gu \et studied the integration of a fluid model
with a packet-level simulation to maintain both the performance
advantage as well as providing details for selected packet
traffic~\cite{Gu04}. They proposed models that can account for the
interaction between flow-based background TCP traffic and foreground
packet-based traffic flows. Another approach called Time-Driven Fluid
Simulation (TDFS)~\cite{Yan99} models traffic nodes in a network as
fluid servers which process workload continuously. In TDFS models,
traffic sources may employ either continuous or discrete event
models. Using TDFS framework, Yan \et studied a single class of
traffic. They also proved that the discretization error does not
depend on the network size but is determined by the length of the
discretization interval used. TDFS can also be easily parallelized to
improve computational performance. All the above mentioned hybrid
approaches show promising avenues for network simulation but have not
yet been used as an alternative to {\bf ns-2}.

A newly proposed hybrid system modeling~\cite{Bohacek03}, which is
different from the aforementioned hybridization approaches, describes
continuous dynamic behavior of a network within an organized finite
state machine (FSM) formalism. Bohacek \et proposed such a framework
for data communication networks and studied various communication
protocols using hybrid systems~\cite{Schaft}. Bohacek \ets work forms
the groundwork for what we are reporting in this paper.


\section{Hybrid Systems Modeling of TCP Networks} % wonsuck
The hybrid system modeling approach for data communication network
simulation is motivated by the drawbacks of the discrete-event driven
packet level simulations and the flow-level simulations. Though the
packet level simulators provide the most detailed network simulation,
they often suffer from scalability issues. It is also difficult to
analyze network performance related to network parameters, such as a
certain combinations of link bandwidth in neighboring nodes and impact
of queuing policies etc.  The flow-level simulation models, on the
other hand, statistically average quantities for better computational
performance and are capable of providing information on network
parameters and network performance relationship. However, the network
dynamics for its entire duration is difficult to model naturally in
the fluid modeling framework because discrete phase changes may be
caused instantaneously by external events, such as a packet drop.

The hybrid system modeling approach~\cite{Bohacek03} proposed by
Bohacek \et effectively reduces the complexity of network simulations
by using flow-level models in each state of the FSM description. Logic
based events trigger transitions between states to capture the
discrete nature of queue dynamics.

In this section, we describe the details of a hybrid system TCP
network model we have developed.  First, a mathematical model of
network is presented. For a TCP network model, we focus on TCP
New-Reno source, drop-tail queuing policy, and their interactions
under FSM models. The models that differ from Bohacek \ets work will
be highlighted. Our solution method differs from Bohacek \ets in that
we use the time-stepped solution method to solve the constructed
hybrid system. This approach resembles the solution techniques used in
TSHS~\cite{Guo00a}.

\subsection{Mathematical Model of a Network} % Wonsuck
Mathematically a network is a directed graph with positive real
numbers assigned to each edge. Thus a communication network is a graph
$(\mathcal{N},\mathcal{L})$ where a set of vertices $\mathcal{N}$
denotes nodes and a set of edges $\mathcal{L}$ denotes a collection of
links connecting two nodes. For instance, two nodes $u \in
\mathcal{N}$ and $v \in \mathcal{N}$ are connected by a link
$l=l_{u,v} \in \mathcal{L}$.  Three real positive numbers are assigned
to each link $l \in \mathcal{L}$, these are: bandwidth (link-rate),
propagation delay, and maximum queue size. Let $B^l$, $\tau^l$, and
$q^{l}_{\rm max}$ denote the bandwidth, the propagation delay, and the
maximum queue size of a link $l$, respectively. These three parameters
characterize each link in a network. Before we discuss the
relationships and the constraints imposed by these parameters, we
first define what we mean by a flow.

Let $\mathcal{F}$ represent a set of {\it flows}. Each
flow $f \in \mathcal{F}$ is associated with a source node $u_s
\in\mathcal{N}$ and a destination node $u_d \in\mathcal{N}$ where $f$
flows from $u_s$ to $u_d$. We assume that flow $f$ is generated and
enters node $u_s$ with an incoming rate $a_f$ at a prescribed time
$t=t^0_f$. Once the network is loaded with data flows, all or part of the
links become path ways for data flows. A flow $f$ on a link $l$ can be
quantified by a link-in-rate (a sending rate), $s_f^l$, and a
link-out-rate (an arrival rate), $a_f^l$.  Bohacek \et named $s_f^l$
as {\it l-link/f-flow} {\it transmission rate} and $a_f^l$ as
{\it l-link/f-flow} {\it arrival rate}~\cite{Bohacek03}.

Note that we assume the source and the destination nodes are unchanged
over time for a given flow $f$.  However, intermediate nodes, that
flow $f$ passes through, may dynamically change in time but at a much
slower scale compared to the data transmission rates. Hence network
dynamics can be well defined on top of any flow path changes within
network.

As discussed earlier our network model is characterized by the three
parameters assigned on a link $l$. Bandwidth $B^l$ imposes an upper
bound on flow sending rate as follows:
\begin{equation}
\sum_{f\in\mathcal{F}} s^l_f \leq B^l, \qquad\forall\  l \in \mathcal{L}~.
\end{equation}
Here $s^l_f = s^l_f(t)$ is a dynamically evolving quantity.

A queue is a temporary storage for flows/packets and is utilized if
the bandwidth of the out-bound link-rate is not enough to process the
flow-rate from the in-bound link. Note that, in our model, a queue is
associated with a link $l\in\mathcal{L}$. $q^l_f$ denotes the size of
the flow $f$ in the queue of a link $l$. Therefore, given the maximum
capacity of the queue $q^l_{\rm max}$ of the link $l$, the following
inequality should hold for all times.
\begin{equation}
\sum_{f\in\mathcal{F}} q^l_f \leq q^l_{\rm max}, 
\qquad\forall\  l \in \mathcal{L}~.
\end{equation}
Here $q^l_f$ is also a function of time, {\it i.e.}, $q^l_f=q^l_f(t)$.

In general the transmission delay is composed of link propagation
delay, processing delay, and queuing delay. Associated with each link
$l$, we only consider a fixed propagation delay $\tau^l$ and a queuing
delay $q^l/B^l$ where $q^l = \sum_{f\in\mathcal{F}} q^l_f$. For the
sake of simplicity, we assume the processing delay is minimal compared
to the other two terms.

Round-trip-time (RTT) is a measure of the current delay on a network
or time elapsed to receive an acknowledgment packet (ACK) from the
destination node.  To compute RTT for a flow $f$, $RTT_f$, the queuing
delay and the propagation delay of the original data packet plus the
propagation delay of the ACK is considered. In other words, we ignore
the queuing delay of the ACK from the destination node to the source
node.  If $\mathcal{L}_f$ is an ordered set of links that flow $f$
takes on its way to the destination, $RTT_f$ becomes
\begin{equation}
RTT_f = \sum_{f\in\mathcal{L}_f}(2\tau^l + \frac{q^l}{B^l})~.
\label{rtt_eqn}
\end{equation}
Here $RTT_f$ and $q^l$ are functions of the time variable $t$ and the
coefficient 2 in front of $\tau^l$ accounts for propagation delay of
the original packet and the ACK. We note that $RTT_f(t)$ is the
estimation of RTT at certain time $t$. In a real system a packet sent
at time $t$ may experience queuing delay different from $q^l(t)/B^l$
at the link $l$ since the packet will exit link $l$ at a certain time
past $t$.

%\subsection{The TCP Protocol}
%\label{sec:TCP}
%In order to study transient behavior of a network, modeling TCP
%protocols has been our first goal since not only is a significant
%portion of internet traffic controlled by the TCP protocol but also it
%potentially possesses rich dynamics due to the concept of varying the
%sending rate in response to the network conditions. We first briefly
%introduce TCP protocol for the readability and completeness of the
%paper.

%TCP protocol~\cite{Jacobson} ensures fair management of network
%resources among geographically spread users. In essence, the TCP
%protocol is implemented over the millions of hosts in the network so
%that a network is somewhat self-adjustable or self-manageable. In this
%subsection we only overview TCP-New Reno policy, which is general
%enough to cover TCP-Tahoe, TCP-Reno, TCP-SACK and TCP-Vegas with minor
%changes in some of the modes.

%The modes of TCP-New Reno policy include slow-start,
%fast-recovery/fast-retransmit (fr/fr), congestion-avoidance and
%time-out. The TCP-New Reno algorithm chooses a certain mode based on
%the limited information it holds and has received from the destination
%host.

%Data transfer starts in the slow-start mode with congestion window
%size, $w$, equal to 1. The source limits the number of sending
%packets below $w$ minus the number of outstanding (unacknowledged)
%packets.  Once the transmitted packet arrives at the destination host,
%the receiver returns a small acknowledgment packet (ACK) with
%information containing the sequence number of the next expected
%packet. If the ACK is successfully received by the sender, the sender
%increases $w$ by ACK so that it has more room to send a new packet. If
%the packet is lost or the order of the packets are altered during
%transmission, the sender may receive duplicate acknowledgments
%(dupACKs) or may not receive any ACKs.

%If the sender receives triple dupACKs, the source enters fr/fr mode.
%Initial congestion window size in fr/fr mode is given by $w^{-}/2+3$
%where $w^{-}$ is the congestion window just before fr/fr mode is
%entered. A variable called slow-start threshold ($ssthr$) is set to
%$w^{-}/2$ which is the maximum congestion window size allowed in the
%next slow-start mode, if there is any. In TCP-New Reno policy, the
%fr/fr mode ends if all the lost packets are successfully
%retransmitted. Subsequently, the congestion avoidance mode follows.

%In the congestion avoidance mode, unlike in slow-start mode, every ACK
%increases the congestion window size by a small portion of the
%congestion window, usually it takes a form of $w^{\rm new} = w^{\rm
%old} + 1/w^{\rm old}$. Normally, a source honors the maximum
%advertised congestion window value satisfying the following
%inequality: $w \leq w^{\rm max}$.  From the congestion avoidance mode,
%depending on number of dupACKs, the subsequent mode is either fr/fr or
%time-out mode which is discussed next.

%If the sender has not received any ACKs, it will not be able to send
%any new packets and the time-out timer will eventually expire. The
%transmission effectively ceases during this waiting
%period. After the time-out period is over a new slow-start mode is
%re-initiated with $w=1$ and $ssthr = w^{-}/2$ where $w^{-}$ is the
%congestion window after the last ACK was received.


TCP protocol~\cite{Jacobson} ensures fair management of network
resources. Modeling TCP protocols has been our first goal since not
only is a significant portion of internet traffic controlled by the
TCP protocol but also it potentially possesses rich dynamics due to
the concept of varying the sending rate in response to the network
conditions.  In the following subsection we focus on TCP-New Reno
policy, other TPC variants can be modeled in similar manner with minor
changes in some of the modes.

\subsection{Source Dynamics}
\label{sec:source_dynamics}
In the following subsection, we revisit Bohacek \ets hybrid TCP source
model~\cite{Bohacek03}. A new model for fr/fr mode is introduced along
with a more accurate sending rate formula. Transition to time-out mode
is naturally embedded without further modeling effort.  We start with
the continuous time (CT) domain model of slow-start mode.  Note that
our primary unknown in CT domain of source dynamics is the congestion
window size.  The sending rate is a derived quantity from the
congestion window size variable.

\subsubsection{Slow-Start} \ \\
The congestion window size $w_f$ of a given flow $f \in\mathcal{F}$
is governed by the following ordinary differential equation:
\begin{equation}
w_f'(t) = \frac{\log m_{ss}}{RTT_f(t)} w_f(t)~,
\label{ss_cwnd_eqn}
\end{equation}
where $\prime$ denotes the time derivative $d/dt$, $m_{ss}$ is a
multiplicative constant such that $w_f$ is being multiplied by $m_{ss}$
every RTT, and $RTT_f$ is computed from Eqn.~(\ref{rtt_eqn}). $w_f$
is initialized to 1 whenever the source enters slow-start mode. Hence,
in slow-start CT domain we solve the initial value ODE problem. Derivation
of Eqn.~(\ref{ss_cwnd_eqn}) is given in the Appendix section.

In TCP flow control $w_f$ is used to limit the number of packets to be
sent out from the sender. Because the sender is able to send more
packets only when the ACKs arrive, effectively, only $w_f$ number of
packets are sent in $RTT_f$ time window. With the assumption that the
packets (or flow $f$) are sent over $RTT_f$ period with evenly spaced
pattern, the instantaneous sending rate $s_f$ should be
\begin{equation}
s_f(t) = \frac{w_f(t)}{RTT_f(t)}~.
\label{sendrate_eqn1}
\end{equation}

Here we should note that the evenly spaced packets over $RTT_f$ time
period is a rather strong assumption. It may be a good model when
traffic flow is up and running for some time so that ACKs become
spread over time due to a certain random process. However, when
traffic is at its early stage, like the very first slow-start mode,
TCP tends to burst out traffic. In such cases the sending rate model,
Eqn.~(\ref{sendrate_eqn1}) might be poor. We revisit this question in
the discussion section.

Bohacek \et used an artificial coefficient $\beta$ to adjust the sending
rate formula as follows: 
\begin{equation}
s_f(t) = \frac{\beta w_f(t)}{RTT_f(t)}~,
\label{sendrate_eqn2}
\end{equation}
where $\beta = 1.45$ is obtained from the trace comparison against
{\bf ns-2} simulation. However, we argue that it might be an artifact
from the hybrid simulation paradigm employed in Bohacek \ets
paper~\cite{Bohacek03}.  In the discussion section we will show how a
CT-FSM ({\it i.e.} hybrid) model may introduce this error by deducing
the value of $\beta$ based on certain assumptions.
% Wonsuck, did we say what these assumptions were - Andy

Though the congestion window size evolves with Eqn~(\ref{ss_cwnd_eqn})
within the slow-start CT domain, $w_f$ should not grow over the
advertised maximum window size $w_f^{\rm adv}$, which is assigned to
the source host by the destination host at the initial hand-shake
stage. Therefore we use Eqn.~(\ref{sendrate_eqn1}) with a slight
modification. The model we use for instantaneous flow sending rate is
\begin{equation}
s_f(t) = \frac{\max \left\{{w_f(t), w_f^{\rm adv}}\right\} }{RTT_f(t)}.
\label{sendrate_eqn3}
\end{equation}

\subsubsection{Congestion Avoidance} \ \\
Congestion avoidance mode increases congestion window size much more
cautiously.  During the congestion avoidance mode every ACK increases
congestion window size by a small portion of the current window size,
\begin{equation}
w_f^{\rm new}  = w_f^{\rm old} + \frac{L}{w_f^{\rm old}}~,
\end{equation}
where $L$ is a real positive number (normally $L=1$). Since $w_f$
packets are sent every RTT, if the congestion avoidance mode is
maintained during this time window, $w_f$ packets will be
returned. Therefore $w_f$ is increased by 1 every RTT if $L=1$. The
ODE model that governs the congestion window size evolution in the
congestion avoidance mode is
\begin{equation}
w_f'(t) = \frac{L}{RTT_f(t)}~.
\label{ca_cwnd_eqn}
\end{equation}
It is straight forward to deduce Eqn.~(\ref{ca_cwnd_eqn}), see the Appendix.

The model shows a linear increase of congestion window size in
congestion avoidance mode as long as no drop is detected and $w_f \leq
w_f^{\rm max}$ is satisfied.  Occurrence of dupACKs or no ACK ends the
congestion avoidance mode and is followed by fr/fr or time-outs.

\subsubsection{Fast Recovery/Fast Retransmit} \ \\
In TCP-New Reno policy fr/fr mode is initiated whenever the sender
receives triple dupACKs. These dupACKs do
not increase the congestion window size in real TCP flow
control. Therefore, it is not accurate if the state remains at either
slow-start or congestion avoidance mode since they are governed by
their own dynamic discipline which evolves the congestion window size
in time.

In order to model congestion window evolution correctly over the state
transition, the source should enter fr/fr once a drop is
detected. However the actual transition should be delayed by the right
amount of time to account for elapsed time between the drop event and
the arrival of dupACKs in the real system.

Since each dupACK increases $w_f$ by 1, the speed of window size
increments is the same as that in the slow-start mode, {\it i.e.}, an
exponential increase, unless there is a drop. Hence
Eqn.~(\ref{ss_cwnd_eqn}) might well be applied as long as no drop
occurs. However, in a real packet network a packet drop is usually
followed by another drop because it takes a certain amount of time
to send buffered packets out onto the out-bound link from a congested
queue. In that event, a hybrid system may experience too many state
transitions in a short time epoch and the model becomes inefficient.

To overcome this modeling difficulty within the hybrid system
framework, we propose an algebraic model for fr/fr mode. Along with it
we provide a formula for the duration of fr/fr mode.

Let $w_f^{-}$ be the congestion window size of flow $f$ just before
the source enters fr/fr mode and $t_0^{\rm fr}$ be the time when the
source mode becomes fr/fr. Suppose $a_f^D(t)$ is the instantaneous
arrival rate of flow $f$ at destination node $D$ at time $t$,
$\zeta(t)$ is size of outstanding flow at time $t$ and $ack(t)$ is the
acknowledged flow size since $t=t_0^{\rm fr}$. Then we have following
relations,
\begin{eqnarray}
ack(t)  &=& \int_{t_0^{\rm fr}}^{t} a_f^D(t'-RTT_f/2) \ dt'~, \nonumber\\
\zeta(t) &=& (w_f^{-} - w_f(t_0^{\rm fr}-RTT_f)) + \eta(t) - H(t)~,\nonumber\\
& &
\end{eqnarray}
for $t \ge t_0^{\rm fr}$. Here $(w_f^{-} - w_f(t_0^{\rm fr}-RTT_f))$
represents the initial outstanding flow size when the source enters
fr/fr, $\eta(t)$ denotes size of flow newly transmitted since the
initiation of fr/fr mode. If there is no drop in $[t-RTT_f(t),t]$,
$H(t)$ is equal to $ack(t)$.  Otherwise, if the source learns about
the drop at $t_{\rm drop} \le t$, $H$ returns the value of $ack(t')$ 
provided no drop is detected in time period of $[t'-RTT_f(t'),t']$ for
$t'\le t_{\rm drop}$. Hence, function $H$ does not count dupACKs in the
events of drop during fr/fr mode. We remark that $\zeta$ approximates
the outstanding flow size for the drop-tail queuing policy.

In fr/fr mode every ACK (either normal ACK or dupACK) increases the
congestion window by 1 so the model for congestion window evolution is
\begin{equation}
w_f(t) = w_f^{-}/2+ ack(t)~.
\end{equation}

The instantaneous sending rate is different from
Eqn.~(\ref{sendrate_eqn1}) since we do not allow the drop event to
trigger state transition. Our flow sending rate model is
\begin{equation}
s_f(t) 
= \frac{\max\left\{{(w_f(t)-\zeta(t)), w_f^{\rm adv}}\right\}}{RTT_f(t)}.
\label{sendrate_eqn4}
\end{equation}

Let $n_f^{\rm drop}$ be the cumulated flow size that are lost during
the transmission. The increment of $n_f^{\rm drop}$ is then given by
\begin{equation}
n_f^{\rm drop} (t) = \int_{t_0^{\rm fr}}^{t} d_f(t') \ dt'~,
\end{equation}
where $d_f$ is a drop rate which will be discussed in detail in the
next section, Sec.~\ref{sec:queue_dynamics}.  Since TCP New-Reno only
retransmits one probably-lost packet at one RTT, the decrement of
$n_f^{\rm drop}$ is 1 at every RTT unless all the flow sent were
lost. In the event of total loss $n_f^{\rm drop}$ should not be
decreased. More precisely, for every $RTT_f$, if $ack(t) \neq
ack(t-RTT_f)$, $n_f^{\rm drop}$ should be decreased as follows:
\begin{equation}
n_f^{\rm drop} (t) = n_f^{\rm drop} (t-RTT_f) - 1~.
\end{equation}
Here $ack(t) \neq ack(t-RTT_f)$ ensures that not all the flows were
dropped during time interval of $[t-RTT_f,t]$.

In our TCP New-Reno implementation, the source exits fr/fr mode if
$n_f^{\rm drop}$ becomes 0. Other TCP variants should differ in fr/fr
dynamic behavior and the duration of fr/fr. For example, TCP Reno
policy exits as soon as the fast retransmit of the first lost packet
is successful. If two packets were dropped, another fr/fr mode follows
immediately after the recovery (fr/fr mode) of the first lost packet.
Therefore, multiple packet losses might lead the source mode to
time-out by multiple steps of halving the congestion window size.

\subsubsection{Time-out} \ \\
In our implementation, whenever a drop is detected, the source enters
fr/fr. In fr/fr, if ACKs cease to return, the congestion window stays
at a constant value and the number of outstanding packets is not
decreased resulting in no new packets being sent for some time period.

Let $t_{\rm out}$ be a timer which is set to 0 when the source enters
fr/fr. $t_{\rm out}$ keeps on increasing with the simulation clock. If
$t_{\rm out}$ exceeds the given time-out period, say, $RTO$, then the
source enters slow-start mode with $w_f = 1$ and $ssthr =
w^{-}_f/2$. Note that if time-out occurs in the simulation, the fr/fr
state entered exhibits no effect on congestion window dynamics because
no new packets were sent and no ACKs were received by the sender.

\subsection{Queue Dynamics}
\label{sec:queue_dynamics}
Suppose $\omega$ and $\nu$ are in-bound (west) and out-bound (east)
links, respectively, in the traffic path of a flow $f$.  Let
$\mathcal{F}_\alpha$ be a set of such flows sharing the links $\omega$
and $\nu$.  We have the following equation describing the queue
changing rate
\begin{equation}
\frac{d}{dt} q^{\omega}_{f} = a^{\omega}_f - s^{\nu}_f - d^{\omega}_f~,
\label{q_eqn}
\end{equation} 
where $q^{\omega}_f$ is the size of flow $f$ buffered at
$q^{\omega}$. $a^{\omega}_f$ is flow arrival rate from the link
$\omega$. $s^{\nu}_f$ denote flow sending rate into the link $\nu$.

Suppose there is no packet loss during transmission over the links in
a network, then incoming flow into a link $\omega$ must all exit from
that link, that is, $a^{\omega}_f(t) = s^{\omega}_f(t-\tau_{\omega})$
where $\tau_{\omega}$ is the link propagation delay. Here
$s^{\omega}_f(t-\tau_{\omega})$ is readily available at time $t$.

Let $q_{\rm max}^{\omega}$ represent maximum queue size associated
with the in-bound link $\omega$. $B^{\nu}$ ($B^{\omega}$) is a given
fixed bandwidth of the link $\nu$ ($\omega$).

Bohacek \et determined the values of $s^{\nu}_f$ and $d^{\omega}_f$
based on queue state. In their model, the FSM for a queue consists of
three states, they are:
\begin{itemize}
\item Empty queue ($q^{\omega} = 0$)
\begin{equation}
d^{\omega}_f = 0, \quad
s^{\nu}_f = \left\{ 
 \begin{array}{l}
   a^{\omega}_f \sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}}
   \le B^{\omega} \\
    
   a^{\omega}_fB^{\omega}/(\sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}})
 \end{array} \right.
\label{q_empty}
\end{equation}
\item Queue neither empty nor full ($0 < q^{\omega} < q^{\omega}_{\max}$ or
 $q^{\omega}=q^{\omega}_{\max}$ but 
 $\sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}} \leq B^{\omega}$)
\begin{equation}
d^{\omega}_f = 0, \quad s^{\nu}_f = 
\frac{q^{\omega}_{f}}{\sum_{\bar{f}\in\mathcal{F}_\alpha}q^{\omega}_{f}} B^{\nu}
\label{q_notfull}
\end{equation}
\item Queue full and still filling ($q^{\omega}=q^{\omega}_{\max}$
and  $\sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}} > B^{\omega}$)
\begin{equation}
d^{\omega}_f = \frac{a^{\omega}_f (\sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}} - B_{\omega})}{\sum_{\bar{f} \in\mathcal{F}_\alpha} a^{\omega}_{\bar{f}}},
\quad s^{\nu}_f = 
\frac{q^{\omega}_{f}}{\sum_{\bar{f}  \in\mathcal{F}_\alpha}q^{\omega}_{f}} B^{\nu}
\label{q_full}
\end{equation}
\end{itemize}
There are two assumptions made to obtain the above queue model. They are,
i) packet arrival rate does not change rapidly and is assumed to be
constant in a short time interval, namely, {\it Arrival Uniformity} ii)
packets of different flows are uniformly distributed, {\it Queue
Uniformity}.

The governing equation~(\ref{q_eqn}) of the queue size and the finite
queue state description (\ref{q_empty}) - (\ref{q_full}) is a
mathematically well-defined hybrid system.  However, the number of
state transitions can be numerous specially if $q^{\omega} \approx
q^{\omega}_{\max}$ to detect zero-crossing for the state
transition~\cite{Schaft} in the hybrid system. Furthermore, a large
system of ODE for $q^{\omega}_{f}$ for all $f\in\mathcal{F}_\alpha$
needs to be solved.

In fact, to find $q^{\omega}_{f}(t)$ from the previously known queue
size $q^{\omega}_{f}(t_0)$ for $t_0 < t$, all that is required is the
net increase of queue size. If we do not impose $q^{\omega}_{f} \leq
q^{\omega}_{\max}$, we have
\begin{equation}
q^{\omega}_{f}(t) = q^{\omega}_{f}(t_0) 
+ \int_{t_0}^{t} (a^{\omega}_f(t') - s^{\nu}_f(t')) \ dt'
\label{q_inteqn}
\end{equation}
It is simply an integral form of Eqn.~(\ref{q_eqn}) without the drop term.
Here $s^{\nu}_f$ is the sending rate into the link $\nu$, {\it i.e.}, 
the sending rate from the $q^{\omega}$. 

We now discuss our queue model. Unlike Bohacek \ets model, the
states are not determined by the queue size $q^{\omega}_{f}$.  Rather they
depend on the bandwidth of the out-bound link, $B^{\nu}$, and the size
of the arriving and departing packet. For the flow sending
$s^{\nu}_f$, we have,
\begin{equation}
s^{\nu}_f = \left\{ 
 \begin{array}{l}
  \frac{q^{\omega}_{f}}{\sum_{\bar{f} \in\mathcal{F}_\alpha}q^{\omega}_{f}},
  \quad \mathrm{if} \ \ 
  \sum_{\bar{f} \in\mathcal{F}_\alpha}q^{\omega}_{f} \geq B^{\omega} \\ \\

  a^{\omega}_{f} +q^{\omega}_{f}, 
  \ \mathrm{if} \ \ 
  \sum_{\bar{f} \in\mathcal{F}_\alpha} (a^{\omega}_{f} +q^{\omega}_{f})
   \leq B^{\omega} \\ \\

  q^{\omega}_{f} + 
  \frac{a^{\omega}_{f}}{\sum_{\bar{f} \in\mathcal{F}_\alpha}a^{\omega}_{f}}
  (B^{\omega} - \sum_{\bar{f} \in\mathcal{F}_\alpha}q^{\omega}_{f}), 
  \ \mathrm{o.w.}
  \end{array} \right.
\label{sendrate_eqn}
\end{equation}
Equations~(\ref{q_inteqn}) and~(\ref{sendrate_eqn}) completely
describe the queue size and the flow sending rate changes in time.
However, the sum of solutions $\sum q^{\omega}_{f}$ may be greater
than $q^{\omega}_{\max}$ which should be understood as buffer overflow
or flow drop.

The drop model should describe which flows suffer losses at the time a
drop occurs. For example, one can drop a larger flow first or vice
versa. A probability distribution $p(f)$ may be used to select a flow
to drop, for instance, $p(f)$ can be
\begin{equation}
p(f) = 
\frac{a^{\omega}_f}{\sum_{\bar{f} \in\mathcal{F}_\alpha}a^{\omega}_{\bar{f}}}~.
\end{equation}

Another question is, if a flow is decided to be dropped, what is the
right size of the flow that should be dropped. We propose the
following model. Let $D^{\omega}_f(t)$ denote the size of dropped flow
$f$, we model $D^{\omega}_f(t)$ as,
\begin{equation}
D^{\omega}_f (t) 
= \max{ \left\{ 0, 
\Bigg\lceil \frac{a^{\omega}_{f}\left( \sum_{\bar{f} \in\mathcal{F}_\alpha}
q^{\omega}_{\bar{f}}(t)-q^{\omega}_{\max}\right)}{\sum_{\bar{f} \in\mathcal{F}_\alpha}a^{\omega}_{\bar{f}}} \Bigg\rceil \right\} }
\label{drop_formula}
\end{equation}

Eqns.~(\ref{q_inteqn})~-~(\ref{sendrate_eqn}) are not entirely
different from the Eqn.~(\ref{q_eqn}) and Eqns.~(\ref{q_empty}) -
(\ref{q_full}).  Eqns.~(\ref{q_inteqn})~-~(\ref{sendrate_eqn}) take an
integral form and $s^{\nu}_f$ is based on packet arrival rate
$a^{\omega}_f$ and queue size $q^{\omega}_{f}$ instead of the queue
state. Note that Bohacek \ets drop rate $d^{\omega}_{f}$ is determined
by queue state only, however, in our model, the value of the drop
size~(\ref{drop_formula}) is a deduced quantity.

\subsection{Solution Techniques}
In this section we describe the solution method and the algorithm used
for our network simulation approach.  We have created the simulator
without relying on any existing simulation frameworks, such as,
Stateflow/Simulink~\cite{simulink}, Modelica~\cite{modelica},
Ptolemy~\cite{ptolemy,ptolemy_paper} etc. There are certain
disadvantages, mostly the amount of work required throughout the
development cycle, however, we possess full control of the whole
system including the solution method of the hybrid system, numerical
scheme to be used in CT domain, and scope of component
model. Furthermore, it is easy to identify any computational
bottleneck and the sources of inaccuracy in the model specific to the
implementation.

First, we take a closer look at the system we have constructed in the
previous sections from CT-FSM point of view. In general FSM is not a
timed model, however, its composition with CT domain requires
detection of any events that may be triggered on the time axis. Each
state in FSM is being refined during the continuous evolution of the
sources and the queues in our model. At each time abscissa, which may
be a priori chosen or dynamically updated, active states are
evaluated. If the state transition condition is enabled, then the FSM
makes corresponding transitions according to the condition and the
system continues to evolve. The details of our solution technique
are described next.

Assume that CT-FSM composition refers a master clock for their
time.  Consider one tick of that master clock to be the time step
$\Delta t$ of our hybrid system simulation. It is trivial that a
smaller $\Delta t$ introduces small error in state transition. In our
simulation, we do not allow FSM state transitions during an epoch of
length $\Delta t$.  Now, CT domain simulation can be performed as
accurately as possible, say, over an interval $[n\Delta t, (n+1)\Delta
t]$ where $n=0,1,\dots$.  At each time point, $n\Delta t$, current
states are evaluated, checked against transition enablers, and changes
are made to the states, as necessary. The tasks described above can be
summarized in the following algorithm.

%\begin{figure}[ht]
\begin{center}
{\bf Algorithm for the Hybrid Simulation}
\end{center}
{\it\small
Input:\\
- A network $(\mathcal{N},\mathcal{L})$,\\
- Bandwidth of all the links $l\in\mathcal{L}$,\\
- Maximum queue size at all the queues,\\
- Set of source nodes $\mathcal{S}$,\\
- Set of destination nodes $\mathcal{D}$,\\ 
- Source wake up time and interval,\\
- Time-step $h$,\\
- End of simulation time $T_{end}$. \\

\begin{tabbing}
Construct $\mathcal{L}_f$ for each $f\in\mathcal{F}$ \\
Set $t = 0$ \\
{\bf while} $t\leq T_{end}$ {\bf do} \\
\quad $t = t + h$ \\
\quad Wake-up sources as scheduled \\
\quad Solve equations for all $s\in\mathcal{S}$ and $q\in\mathcal{N}$ \\
\quad {\bf for} all $s\in\mathcal{S}$ {\bf do} \\
\quad\quad {\bf switch} mode(s): \\
\quad\quad {\bf case} (mode(s)=slow start):\\
\quad\quad\quad {\bf if} cwnd(s) $\geq$ {\it ssthr} {\bf then} \\
\quad\quad\quad\quad $\mathrm{mode(s)} \stackrel{t+\mathrm{delay}}{\longrightarrow}$ c.a. \\
\quad\quad\quad {\bf else if} flow $f$ from $s$ suffer drop {\bf then} \\ 
\quad\quad\quad\quad $\mathrm{mode(s)} \stackrel{t+\mathrm{delay}}{\longrightarrow}$ fr/fr \\
\quad\quad\quad {\bf else} \\
\quad\quad\quad\quad remain in slow start mode\\
\quad\quad {\bf case} (mode(s)=fr/fr):\\
\quad\quad\quad {\bf if} $n^{\rm drop}_f=0$ {\bf then} \\
\quad\quad\quad\quad $\mathrm{mode(s)} \stackrel{t}{\longrightarrow}$ c.a.\\
\quad\quad\quad {\bf else} \\
\quad\quad\quad\quad remain in fr/fr \\
\quad\quad {\bf case} (mode(s)=c.a.):\\
\quad\quad\quad {\bf if} flow $f$ from $s$ suffer drop {\bf then} \\
\quad\quad\quad\quad $\mathrm{mode(s)} \stackrel{t+\mathrm{delay}}{\longrightarrow}$ fr/fr\\
\quad\quad\quad {\bf else} \\
\quad\quad\quad\quad remain in c.a. \\
\quad {\bf end for}\\
\quad {\bf for} all queue {\bf do} \\
\quad\quad evaluate and set queue state for the next iteration \\
\quad {\bf end for}\\
{\bf end while}
\end{tabbing}
}
%\caption[]{Algorithm}
%\label{algorithm}
%\end{figure}

Now we comment on a numerical scheme for 
Eqns.~(\ref{ss_cwnd_eqn})~-~(\ref{ca_cwnd_eqn}) which describes 
dynamics of the congestion window evolution.  We have, in their general form,
\[
w' = g(t,w), \quad w(0) = w_0
\]
where $w_0$ is given.

It is well known that Euler's method~\cite{Conte} is not practical 
due to the requirements of small step size for reasonably accurate
numerical solutions. Runge-Kutta methods do not require 
higher order derivative by evaluating $g(t,w)$ at selected points on each
interval of computation.  Higher order methods can be obtained by paying
extra for the evaluation of $g$. For example, Runge-Kutta method of order 2
employs the following recursive formula:
\[
\begin{array}{l}
 w_{n+1} = w_n +\frac{1}{2}(k_1+k_2), \\
 \qquad\mathrm{where}\quad k_1 = hg(t_n, w_n), \\
 \qquad\qquad\quad\  k_2 = hg(t_n+h, w_n+k_1)~.
\end{array}
\]
Here $t_n=nh$ and $w_n = w(t_n)$. The local error of this method is of
order $h^3$ while Euler's method is of order $h^2$. Henceforth, we can
use larger step size or save computation time with higher order
numerical scheme.  As one might expect, in terms of computational
efficiency, we found that the larger the hybrid simulation time step,
the better it is with a higher order method.

Finally, we remark that the proposed solution approach is somewhat
similar to the Time-stepped Hybrid Simulation (TSHS)~\cite{Guo00a} and
shares the idea of packet smoothing, however, the concept of hybrid
system of CT-FSM composition is fundamentally different from TSHS.

\section{Hybrid Model \& Simulation Validation}% Marina, Andy, Amogh

\subsection{Simulation Set Up \& Environment} % Andy/Amogh
Our simulation studies are comprised of multiple experiments using the
packet-based {\bf ns-2} simulator and the hybrid model-driven
simulator that we have developed. Since we are interested in comparing
the two simulators we had to use the same traffic characteristics for
both simulators. Furthermore the specific choice of the traffic pattern
is dictated by how well it models network traffic. Experimental studies
have shown that individual sources of traffic can be approximated by
ON/OFF processes, where the On and Off periods have subexponential
distributions. Within the ON state, a source transmits at peak rate. We
use this model for our experiments with both the {\bf ns-2} and the
hybrid simulations. Packets are sent at a fixed rate during ON state
while no packets are sent during OFF state. Both on and off periods are
taken from a Pareto distribution with constant size packets. These
traffic sources can be used to generate aggregate network traffic that
exhibits the long range dependency seen in the IP network traffic.

%Our simulation testbed comprises multiple experiments using the
%packet-based {\bf ns-2} simulator and the hybrid model-driven
%simulator we have developed. Since we are interested in internet level
%measurements, it was important that we chose the right traffic
%characteristics. Moreover, for our results comparing both types of
%simulators we had to use the same traffic characteristics for both 
%simulators. Choice of the right traffic is thus crucial for our
%experiments.

%Recent studies have revealed that several types of network traffic
%exhibit multiple time-scale behavior, where considerable
%fluctuations in the traffic rate are observed over several time scales
%to be above or under the average rate. A Poisson or Markovian traffic
%model for the arrival process would have a characteristic burst length
%that smoothes out by averaging over a long enough time
%scale. Unfortunately, this model is inappropriate since real network
%traffic indicates that significant traffic burstiness is still present
%over a wide range of time scale. Experimental studies have shown that
%such behavior can be modeled using subexponential distributions over
%multiple time scales. A subexponential distribution has the property
%to decay slower than any exponential distribution and is shown to fit
%many network traffic parameters including inter arrival time for
%requests, file sizes and duration of sessions.

%Experimental studies have shown that individual sources of traffic can
%be approximated by ON/OFF processes, where the On and Off periods have
%subexponential distributions. Within the ON state, a source transmits
%at peak rate. We use this model for our experiments with both ns-2 and
%hybrid simulations. Packets are sent at a fixed rate during ON state
%while no packets are sent during OFF state. Both on and off periods
%are taken from a Pareto distribution with constant size packets. These
%sources can be used to generate aggregate traffic that exhibits long
%range dependency.

\subsection{Simulation Validation \& Results} % Marina
The first goal of our experimental set up was to baseline the
performance of our hybrid model using well known packet level
simulators such as {\bf ns-2}. For this purpose we used two topologies. A
simple topology as shown in Figure~\ref{simpletop} with a single TCP New-Reno
source and a dumbbell topology as shown in Figure~\ref{dumbtop} with 4
TCP New-Reno sources.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/simpletop.eps}
%\psfig{figure=figs/simpletop.pdf, height = 3}
\end{center}
\caption{A Simple Topology with 1 TCP New-Reno Source.}
\label{simpletop}
\end{figure}
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/dumbbelltop.eps}
%\psfig{figure=figs/dumbbelltop.pdf, height = 3.5in}
\end{center}
\caption{A Dumbbell topology with 4 TCP New-Reno Sources.}
\label{dumbtop}
\end{figure}
For the {\bf ns-2} simulations we used an identical TCP source model
with no background traffic. The simple topology is used to evaluate
the behavior of the queue dynamics as well as TCP source behavior. It
consists of a single edge and two network nodes. The link propagation
delay was 2ms. In the dumbbell topology, we have 4 TCP
sources that are aggregated at the queue for router 1 and each of
these sources have link delays of 1, 2, 3, and 4ms each. All four TCP
sources are destined for the same destination node {\bf d}. The queue
size is set to 15 packets with each packet consisting of 512
bytes. The run time for each of the simulations was 1000ms. The TCP
flows were triggered to start at 5ms. The metrics used for comparing
the hybrid simulator with the {\bf ns-2} simulator are, the queue
sizes and the congestion window sizes of the TCP source.

\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/cwndcomp.ps}
%\psfig{figure=figs/cwndcomp.pdf, height = 3.5in}
\end{center}
\caption{Comparison of congestion window evolution from the hybrid
model and {\bf ns2} simulation. TCP New-Reno Sources.} 
\label{cwndcomp}
\end{figure}
Figure~\ref{cwndcomp} shows a comparison of the congestion window
evolution of {\bf ns-2} versus the hybrid model. Observe that the
behavior of the two simulators is more closely matched in the
congestion avoidance phase of the simulation while there is some
mismatch in the slow start phase as well as the fr/fr phases. In the
case of the hybrid model the exponential increase in the slow start
phase reaches a value of 55 while in {\bf ns-2} the maximum value of the
congestion window in the slow start phase is only about 27. Also the
increase in congestion window size in the hybrid model is smooth due
to the sending rate model that is used. In the hybrid model we use a
time averaged sending rate for the packets. The plausible explanation
for the deviation in the behavior of {\bf ns-2} is that {\bf ns-2}
employs a burst model sending rate in the initial slow start stages of
the TCP flow. Also in the fr/fr mode we observe that the hybrid model
is more graceful in evolving the congestion window while in {\bf ns-2} there
is an abrupt change in the window size. 

The graceful congestion window
evolution in the hybrid simulation is attributed to the wait time that
is incorporated into the model after a drop occurs to account for the
arrival of dupACKs. Furthermore, in the congestion avoidance phase we
observe that the saw tooth behavior is more well behaved in the hybrid
model as compared to ns. This smooth behavior is attributed to the
lack of randomness in the arrival rates at the individual queues. In
{\bf ns-2} the randomness arises from the sending rate model that is
used by the TCP sources.
% Marina - so is this still alright i.e., in ns-2 there is a
% randomness but not in hybrid - so what is real? - Andy

\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/QSizeComp.ps}
%\psfig{figure=figs/QSizeComp.pdf, height = 3.5in}
\end{center}
\caption{Comparison of queue size changes using the hybrid model and
{\bf ns-2}. Queue size changes are evaluated using different time step values
(0.25, 0.5, 1.0) in the hybrid simulator.} 
\label{Qsizecomp}
\end{figure}
Figure~\ref{Qsizecomp} shows the queue size changes using different
time steps in the hybrid simulator as well as in {\bf ns-2}. The queue sizes
computed using the finer time steps in the hybrid model are much
closer to the actual behavior. However using a finer time step
increases the computational complexity. The queue size changes
obtained by the hybrid simulator using a time step of 0.25 provides
the closest match to the queue size change observed in {\bf ns-2}.

In the dumbbell topology we simulate 4 TCP flows each with links of
different delay values but having the same destination node. We do not
assume any background traffic.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/CwndHybDumbComp.ps}
%\psfig{figure=figs/CwndHybDumbComp.pdf, height = 3.5in}
\end{center}
\caption{Congestion window sizes evolution over time: 4 sources in the dumbbell topology using the hybrid simulation.}
\label{CwndHybDumbComp}
\end{figure} 
Figure~\ref{CwndHybDumbComp} shows the evolution of the congestion
window of the 4 different sources. Note that as expected in steady
state all of the TCP sources are synchronized to a saw tooth pattern
with a period of 100ms. The {\bf ns-2} simulation for a similar topology is
shown in Figure~\ref{CwndNsDumbComp}.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/CwndNsDumbComp.ps}
%\psfig{figure=figs/CwndNsDumbComp.pdf, height = 3.5in}
\end{center}
\caption{Congestion window sizes evolution over time: 4 sources in the dumbbell topology using {\bf ns-2}.}
\label{CwndNsDumbComp}
\end{figure}
Note that in the steady state a good match is obtained for congestion
window evolution in the hybrid as well as the {\bf ns-2} models. The average
values of the steady state congestion window sizes obtained using the
hybrid model and {\bf ns-2} are in Table~\ref{tab:cwnd}.
\begin{table}[!h]
\begin{center}
 \begin{tabular}{|c||c|c|c|c|}\hline
 \ & ${\rm source}_1$ & ${\rm source}_2$ & ${\rm source}_3$ & ${\rm source}_4$ 
 \\ \hline \hline
 hybrid & 26.92 & 19.92 & 15.89 & 13.84 \\ \hline
 ns-2   & 26.92 & 19.82 & 16.02 & 13.94 \\ \hline
 \end{tabular}
\end{center}
\caption{Dumbbell topology, steady state ($t \geq$ 500ms) average congestion window size.}
\label{tab:cwnd}
\end{table}
Figure~\ref{QSizeDumb} compares the queue size changes between the
hybrid and the {\bf ns-2} simulations. As the simulation time progresses and
reaches a steady state we obtain a close match between the queue
behavior of the {\bf ns-2} and the hybrid models.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/QSizeDumb.ps}
%\psfig{figure=figs/QSizeDumb.pdf, height = 3.5in}
\end{center}
\caption{Comparison of the queue size changes: the hybrid model {\sl vs.}
{\bf ns-2}.}
\label{QSizeDumb}
\end{figure} 

\section{Simulation of Large Network Topologies}
In order to study the effectiveness of the hybrid simulator we
evaluated its performance using two different network
topologies. These two topologies differ significantly in terms of the
number of nodes and their degree of connectivity. Similar topologies
were simulated in {\bf ns-2} and we compare the behavior of queue
sizes under different TCP source configurations.

\subsection{Network Topology} 
Figure~\ref{UStop}, called to be {\it US Topology}, shows a network of 8
nodes connected by duplex-links of 100 Mbps bandwidth and 2ms
link-delay. Links 2-3 and 8-7 have higher, 1000 Mbps, bandwidth with 2ms
link-delay also.  For all the nodes, maximum queue size is set to 50
packets. Queue sizes are traced for the link 6-7 (L1), the link 4-5
(L2), and the link 2-3 (L3). We have used 53, 86, and 129 TCP sources
where destination nodes are spread over the network.

A network, named {\it KL Topology}, shown in Figure~\ref{KLtop},
consists of 15 nodes connected by links of 100 Mbps and 1000 Mbps
bandwidths. The propagation delay is 4ms for all the the links in the
network.  Maximum queue size used is 75 packets for all the nodes.
Total 160 sources were simulated in this topology.
% Amogh - Figures and Description
% figures from Amogh added - Andy
\begin{figure}[htb]
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/us_top.eps}
%\psfig{figure=figs/us_top.pdf, height = 3.5in}
\end{center}
\caption{US Topology loaded with 53 sources: A dashed-numbers ({\it
e.g.} 1-8-7-6) denotes the path of a flow. The parenthesized number
represents how many flows are in that path. }
\label{UStop}
\end{figure}

\begin{figure}[htb]
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/kl_top.eps}
%\psfig{figure=figs/kl_top.pdf, height = 3.5in}
\end{center}
\caption{KL Topology: D denotes the destination node of the queue for which 
we present the simulation results.} 
\label{KLtop}
\end{figure}

\subsection{Results} 

\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/Q3SizeComp.ps}
%\psfig{figure=figs/Q3SizeComp.pdf, height = 3.5in}
\end{center}
\caption{Comparison of queue size changes using the hybrid model and
{\bf ns-2} for US topology. Queue size changes are evaluated using different number of sources.  Legend: $\cdots$ {\bf ns-2} and -- hybrid} 
\label{Q3SizeComp}
\end{figure}
Figure~\ref{Q3SizeComp} shows the change in queue (Q2) size under three
different source configurations (53, 86 and 129 sources each) for both
the {\bf ns-2} as well as the hybrid simulator.  In our simulation
scenarios, of all the queues in the network this particular queue was
the one that had the most discrepancy with the {\bf ns-2} simulation.
However note that even this worst case deviation from {\bf
ns-2} is not very significant as can be seen from Table~\ref{tab:q3}.
\begin{table}[!h]
\begin{center}
 \begin{tabular}{|c||c|c|c|}\hline
\  & ${\rm Q}_1$ & ${\rm Q}_2$ & ${\rm Q}_3$\\ \hline\hline
53  connections & (23.7:24.8)  & (21.0:22.3)  & (8.7:10.5) \\ \hline
86  connections & (24.4:23.8)  & (20.5:22.5)  & (11.8:13.0) \\ \hline
129 connections & (28.5:28.5)  & (13.9:16.4)  & (21.0:20.6) \\ \hline
 \end{tabular}
\end{center}
\caption{US topology average queue length for Q1, Q2, Q3 written for (hybrid:ns-2)}
\label{tab:q3}
\end{table}

Figure~\ref{Q3SizeKL} shows the queue size variation between {\bf ns-2}
and the hybrid simulator for {KL Topology} loaded with 160 TCP sources.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/QSizeKL.ps}
%\psfig{figure=figs/QSizeKL.pdf, height = 3.5in}
\end{center}
\caption{Comparison of queue size changes using the hybrid model and
ns for KL topology. See Figure~\ref{KLtop} for queue location. Legend:
$\cdots$ {\bf ns-2} and -- hybrid }
\label{Q3SizeKL}
\end{figure}
Once again we observe a reasonably good match between our simulator
and {\bf ns-2}. The average queue size for the hybrid simulator was 25
and that of {\bf ns-2} was 24.

\section{Discussion \& Conclusion} % Everyone
In this section we discuss the sending rate model~(\ref{sendrate_eqn1})
and the hybrid model implementation and complexity issues. We conclude 
the section with a brief summary of the paper.
 
\subsection{Transport Layer \& Application Layer} % Wonsuck
In Sec.~\ref{sec:source_dynamics} we pointed out that the sending rate
model~(\ref{sendrate_eqn1}) might be inappropriate due to the bursty
nature of the TCP protocol and the dependency on packet availability
from applications. For example, in {\bf ns-2} simulation the queue
overflows much earlier than the hybrid simulation (see
Figure~\ref{cwndcomp}).  It is because the
application sent all the packets in a much shorter period of time
than the sending rate modeling equation ~(\ref{sendrate_eqn1}).

The model in equation~(\ref{sendrate_eqn1}) sends out packets (flow)
over a time interval of length RTT. More specifically, if congestion
window size is 30 and RTT is 6ms then 5 packets are sent per 1ms in
our model. {\bf ns-2} is more likely to send all the 30 packets almost
instantaneously. Hence if the bandwidth of the out-bound link is not
capable of processing all the incoming traffic, the queue
overflows. However, with the model in equation~\ref{sendrate_eqn1},
the queue has enough time to process the flows that are arriving
spread out over time.

The model in equation~(\ref{sendrate_eqn1}) is quite appropriate at
the asymptotic stage, when the packets are more evenly spread out.
One can observe that in Figure~\ref{cwndcomp}, the plots from the hybrid
simulation and {\bf ns-2} match well for $t\ge q$ 300ms.

%A solution to this problem cannot be presented in transport layer only,
%because the available packets are determined by application and an
%operating system (OS) running on a host. At the least, the application
%layer model should be available to correctly simulate real network
%along with transport layer model.  The application layer only prepare
%the packets depending the chosen application and the transport layer
%process the given information without any modification.

\subsection{Implementation Issues} % Wonsuck
Our current hybrid system simulator employs a fixed time-step
approach.  More specifically, with the notion of a global clock, time
marches by a fixed time-step toward the end of simulation time. During
the time advancement all the sources and the queues are computed and
checked against the transition enablers. The checks are carried out
regardless of the states and the events at any time interval. Hence
the simulation advancing method might be inefficient. More intelligent
time marching techniques are necessary for further speed up.

\subsection{Conclusions} % Marina & Andy
In this paper, we have demonstrated the strength of the hybrid system
modeling and simulation method for TCP network. New models,
enhancements, and revisions were proposed and studied. We have used
the time-stepped solution method to solve the hybrid system model of a TCP
network. The extensive experiments showed good agreements with the
results from Network Simulator {\bf ns-2}. However, we have not yet
studied the stability and the convergence issues of the numerical
methods we have employed. The robust simulation tool which we have developed
can be used to analyze complex network behavior so that 
service providers are better equipped to support, for instance, 
low-latency quality of service (QoS) needs of their applications.


\setcounter{subsection}{0}
\section*{Appendix} % Wonsuck
\subsection{Equation for Congestion Window Size}
Consider slow-start mode in TCP source dynamics. In slow-start mode,
the congestion window size increases exponentially in time.  Suppose
that the window size is multiplied by $m_{ss}\in R^{+}$ for every RTT,
where $R^{+}$ is the set of real positive number. Usually $m_{ss}=2$
since most of of the TCP protocols double the congestion window size
for every RTT.  Let $RTT = n\Delta t$ for $\Delta t\in R^{+}$ and
$n=0,1, \dots$~.

We assume that congestion window $w$ uniformly increases over $RTT$
time period $[t,t+RTT]$, that is, $w(t+\Delta t) = \sqrt[n]{m_{ss}}
w(t)$.  Then we have,
\begin{equation}
\frac{w(t+\Delta t)-w(t)}{\Delta t} 
= w(t) \frac{(m_{ss}^{\Delta t/RTT} - 1)}{\Delta t}~.
\label{w_diff}
\end{equation}

Take $\lim_{\Delta t\to 0}$ on both sides of Eqn.~(\ref{w_diff}) and
applying L'Hospital's rule. We obtain,
\[
w'(t) = \ln(m_{ss})\frac{w(t)}{RTT}~. 
\]

The governing equation for congestion avoidance mode,
Eqn.~(\ref{ca_cwnd_eqn}) can be obtained in a similar manner.

\subsection{Correction Coefficient $\beta$}
Bohacek \et used a correction factor $\beta=1.45$ for the sending rate
formula~(\ref{sendrate_eqn2}) to best match their simulation results
to the ones from {\bf ns-2}. Without the specifics of the numerical
method used in Bohacek \ets work, we conjecture that the presence of
$\beta$ is necessary if time points for event checking and triggering
is RTT apart.

Let $m_{ss}=2$ for the following discussion. Consider slow-start mode
which is continued from $t=t_0$. Then congestion window $w$ becomes
$2^n$ at $t=t_0+nRTT$. If a fixed $w$ is used throughout the time
interval $[nRTT, (n+1)RTT]$, the total number of packets sent during
this time period is exactly $w$ provided that
Eqn.~(\ref{sendrate_eqn2}) is used for the packet sending process.
\begin{figure}
\begin{center}
\epsfxsize=7.5cm \epsffile{figs/beta_explain.eps}
%\psfig{figure=figs/beta_explain.pdf, height = 3.5in}
\end{center}
\caption{Schematic description congestion window size changes for the
TCP controlled flow in slow-start mode. With infinitesimal packets the
congestion window size supposedly follows the solid curve. The
step-wise increase may represent what happens in a real network. The
dashed line indicates a fixed value of window size during RTT time
period.}
\label{fig:beta}
\end{figure}

However, in reality, the congestion window size increases during this time
interval due to arriving ACKs, therefore the source is able to send more
packets.  It can be explained schematically in
Figure~\ref{fig:beta}. The dashed line represents fixed
window size ($w=2^n$) over $[nRTT, (n+1)RTT]$, the step-wise variation may
represent what happens in reality, and the solid curve is a plot of
$y=2^{x/RTT}$.

With the sending rate model~(\ref{sendrate_eqn2}), {\it i.e.} $r =
w/RTT$, the area between the dashed line and the time axis represents
the size of flow sent in $[nRTT, (n+1)RTT]$, where $n=4$ in the
figure. It is always smaller than the area below the step-wise
function provided the slow-start mode persists.  Suppose that the
congestion window size variation can be closely approximated by 
$y=2^{x/RTT}$ curve, {\it i.e.} with infinitesimal $\Delta t$, the area
under the $y=2^{x/RTT}$ curve for $[nRTT, (n+1)RTT]$ is
\[
\int_{nRTT}^{(n+1)RTT} 2^{x/RTT} dx = 2^n\frac{RTT}{\ln(2)}~.
\]
Note that $(\ln(2))^{-1} \approx 1.4424$ which is close to the value of
$\beta = 1.45$ reported in~\cite{Bohacek03}.

%\section*{Acknowledgments}
%The authors would like to thank Lawrence Cowsar and Sanjoy Paul both
%at Bell Labs for their help in bringing this project to fruition.

\begin{thebibliography}{1}
\bibitem{Clark:90}
D.D. Clark and D. L. Tennenhouse,
\newblock {``Architectural Considerations for a New Generation of Protocols'',}
\newblock{ACM SIGCOMM '90 (Philadelphia, PA, September 1990)}

\bibitem{ns2}
The VINT Project on NS2,
\newblock {``The Network Simulator'', ns2: http://www.isi.edu/nsnam/ns}

\bibitem{ssfnet}
SSFNet,
\newblock {``The Scalable Simulation Framework'', 
http://www.ssfnet.org/homePage.html}

\bibitem{jsim}
J-Sim,
\newblock {``A component-based simulation environment'', 
http://www.j-sim.org}

\bibitem{inet}
INET,
\newblock {``Network simulation in the J-Sim environment'', 
http://www.j-sim.org/drcl.inet/index.html}

\bibitem{Misra99a}
V. Misra, W-B. Gong, and D. Towsley,
\newblock {``Stochastic Differential Equation Modeling and Analysis of TCP Windowsize Behavior'',}
\newblock{Performance '99 (Istanbul, Turkey, October 1999)}

\bibitem{Misra99b}
V. Misra, W-B. Gong, and D. Towsley,
\newblock {``A Fluid-based Analysis of a Network of AQM Routers Supporting TCP Flows with an Application to RED'',}
\newblock{Proceedings of ACM SIGCOMM'00, (Stockholm, Sweden, September 2000)}

\bibitem{Altman}
E. Altman, K. Avrachenkov, and C. Barakat,
\newblock{``A Stochastic Model of TCP/IP with Stationary Random Losses'',}
\newblock{ACM SIGCOMM Computer Communication Review, 30: 231-242, 2000}

\bibitem{Bonald}
T. Bonald,
\newblock{``Comparison of TCP Reno and TCP Vegas via fluid approximation'',}
\newblock{Technical Report RR-3563, INRIA, 1998}

\bibitem{Jacobson}
V. Jacobson,
\newblock{``Congestion Avoidance and Control'',}
\newblock{Computer Communication Review, vol. 18, 1988}

\bibitem{Guo00a}
Y. Guo, W-B. Gong, and D. Towsley,
\newblock{``Time-stepped Hybrid Simulation (TSHS) for Large Scale Networks'',}
\newblock{Proceedings of IEEE Infocom 2000, 441-450, 2000}

\bibitem{Gu04}
Y. Gu, Y. Li, and D. Towsley,
\newblock{``On Integrating Fluid Models with Packet Simulation'',}
\newblock{Proceedings of IEEE Infocom 2004, 2004}

\bibitem{Yan99}
A. Yan and W-B. Gong,
\newblock{``Time-Driven Fluid Simulation for High-Speed Networks'',}
\newblock{IEEE Trans. on Information Theory, vol. 45, no. 5, July 1999}

\bibitem{Bohacek03}
S. Bohacek, J.P. Hespanha, J. Lee, and K. Obraczka,
\newblock{``A Hybrid Systems Modeling Framework for Fast and Accurate 
Simulation of Data Communication Networks'',}
\newblock{ACM SIGMETRICS'03, June 2003}

\bibitem{Schaft}
A.J. Van Der Schaft and J.M. Schumacher,
\newblock{An Introduction to Hybrid Dynamical Systems (Lecture Notes in Control and Information Sciences, 251)}
\newblock{Springer-Verlag Telos, 1999}

\bibitem{Conte}
S.D. Conte and C. Boor,
\newblock{``Elementary Numerical Analysis (An Algorithmic Approach)''}
\newblock{McGraw-Hill, 1981}

\bibitem{simulink}
The MathWorks Inc.,
\newblock{``Simulink'', http://www.mathworks.com/products/simulink}

\bibitem{modelica}
The Modelica Association,
\newblock{``Modelica: '', http://www.modelica.org}

\bibitem{ptolemy}
Dept. of EECS, University of California at Berkeley, 
\newblock{``The Ptolemy Project'', http://ptolemy.eecs.berkeley.edu}

\bibitem{ptolemy_paper}
E. Lee 
\newblock{``Overview of the Ptolemy Project''},
\newblock{Technical Memorandum No. UCB/ERL M03/25, University of California, Berkeley, CA, 94720}

\end{thebibliography}

\end{document}