%        File: VTthesis_template.tex
%     Created: Thu Mar 24 11:00 AM 2016 EDT
% Last Change: Wed Apr 06 02:37 PM 2016 EDT
%      Author: Alan M Lattimer, VT
%
% Instructions for using the document are embedded in the document.  Additionally, there
% are comments at the end of the file that give some suggestions on writing your
% thesis.  To get started you should add VTthesis.cls to your path or keep it in
% the same directory as this file.  In addition to the standard formatting options
% the following options are defined for the VTthesis class:
%  proposal, prelim, doublespace, draft
\documentclass[doublespace]{VTthesis}
% Using the following header instead will create a draft copy of your thesis
%\documentclass[doublespace,draft]{VTthesis}

% The lipsum package is just included to put dummy text in the document in order to demonstrate 
% page headers and table of contents behavior.  You should remove it once you begin 
% writing your actual thesis or dissertation
\usepackage{lipsum}
\usepackage{pgfplots}
\pgfplotsset{compat=1.14}
\usepackage{algorithm}
\usepackage{algorithmic}
\newcommand\tab[1][1cm]{\hspace*{#1}}

% Title of your thesis
\title{Online Message Delay Prediction for Model Predictive Control over CAN}

% You should include 3-5 keywords, separated by commas
\keywords{Model Predictive Control (MPC), Controller Area Network (CAN), online prediction, delay}

% Your name including middle initial
\author{Amith Kaushal Bangalore Narendranath Rao}

% Change this to your program, e.g. Physics, Civil Engineering, etc.
\program{Computer Engineering} 

% Change this to your degree, e.g. Masters of Science, Masters of Art, etc.
\degree{Master of Science} 

% This should be your defense date
\submitdate{June 16, 2017} 

% Committee members. Only have five readers and one chair available
% only use the ones you need and don't include the ones you don't need.
% You can also declare a Co-advisor. If you do, the principal and co-advisors
% will be listed as co-advisors on the title page.  Per the VT ETD standards, 
% you should not include titles or educational qualifications such as PhD or Dr.
% You should however include middle initials if possible.
\principaladvisor{Haibo Zeng}
%\coadvisor{Vicente Esparza}
\firstreader{Pratap Tokekar}
\secondreader{Feng Guo}
%\thirdreader{James Smith}
%\fourthreader{Fourth Committee Member}
%\fifthreader{Fifth Committee Member}

% The dedication and acknowledgement pages are optional.  Uncomment to 
% include them in your thesis or dissertation
\dedication{To my parents, my sister, and my grandparents for the unconditional love and support.}
\acknowledge{ 
I would firstly like to express my heartfelt gratitude to my advisor/committee chair, Dr. Haibo Zeng for guiding me throughout my time here; from the time when I took his class, to when I served as a TA for a class he taught and finally as a Research Assistant at CESCA, and steering me in the right direction. He has always been more than accommodating with respect to my personal and academic requirements and has helped me grow towards being a professional as well.
\newline
I want to thank the experts on my review committee, Dr. Pratap Tokekar and Dr. Feng Guo, for taking time out from their busy schedules to read, review and provide valuable comments.
\newline
I would like to extend my special thanks to Dr. Anton Cervin (Lund University) for clarifying queries about TrueTime and Dr. Liuping Wang (RMIT University) for the discussions on MPC which proved invaluable in my research. 
\newline
I thank my fellow researchers in the Lab who had to put up with me switching workstations every two months to accommodate my research requirements. They were all very helpful and particularly nice to me.
\newline
Further, I thank my parents and my sister, for always believing in my abilities more than I ever have, thereby driving me to reach heights greater than my imagination. 
\newline
Being 16,000 miles away from my family back in Bangalore, I would like to thank my close friends who have become my family here in Blacksburg over these two years. They have shared my struggles and their joys, laughed and cried with me, partied and studied with me, and I can safely say that the journey would have been very different without them. I would like to make a special mention: Shriya - from my family here, who gave me a very practical insight into Control Theory, which was not my area of focus. Because of those sessions, my research now includes the advanced Model Predictive Control. 
\newline
Lastly, I would like to thank the entirety of Virginia Tech for deeming me worthy of a Master's Degree and fully taking care of my tuition right from Day 1. It made a big difference in my life. The Hokie community, in all its diversity has given me so many opportunities to grow as an individual and develop a well rounded personality. I cannot thank my stars enough for it.  
\newline
\newline
Amith Kaushal Bangalore Narendranath Rao
}

% The abstract is required and should be <=250 words for thesis, <=350 words for disseration
\abstract{Today's Cyber-Physical Systems (CPS) are typically distributed over several computing nodes communicating by way of shared buses such as Controller Area Network (CAN). Their control performance gets degraded due to variable delays (jitters) incurred by messages on the shared CAN bus due to contention and network overhead. This work presents a novel online delay prediction approach that predicts the message delay at runtime based on real-time traffic information on CAN. It leverages the proposed method to improve control quality, by compensating for the message delay using the Model Predictive Control (MPC) algorithm in designing the controller. By simulating an automotive Cruise Control system and a DC Motor plant in a CAN environment, it goes on to demonstrate that the delay prediction is accurate, and that the MPC design which takes the message delay into consideration, performs considerably better. It also implements the proposed method on an 8-bit 16MHz ATmega328P microcontroller and measures the execution time overhead. The results clearly indicate that the method is computationally feasible for online usage.}

% The general audience abstract is required. There are currently no word limits.
\abstractgenaud{In today's world, most complicated systems such as automobiles employ a decentralized modular architecture with several nodes communicating with each other over a shared medium. The Controller Area Network (CAN) is the most widely accepted standard as far as automobiles are concerned. The performance of such systems gets degraded due to the variable delays (jitters) incurred by messages on the CAN. These delays can be caused by messages of higher importance delaying bus access to the messages of lower importance, or due to other network related issues. This work presents a novel approach that predicts the message delays in real-time based on the traffic information on CAN. This approach leverages the proposed method to improve the control quality by compensating for the message delay using an advanced controller algorithm called Model Predictive Control (MPC). By simulating an automotive Cruise Control system and a DC motor plant in a CAN environment, this work goes on to demonstrate that the delay prediction is accurate, and that the MPC design which takes the  message delay into consideration, performs considerably better. It also implements the proposed approach on a low end microcontroller (8bit, 16MHz ATmega328P) and measures the time taken for predicting the delay for each message (execution overhead). The obtained results clearly indicate that the method is computationally feasible for use in a real-time scenario.  }

\begin{document}
	% The following lines set up the front matter of your thesis or dissertation and is required to ensure proper
	% formatting per the VT ETD standards. 
  \frontmatter
  \maketitle
  \tableofcontents

	% The list of figures and tables are now optional per the official ETD standards.  Unless you have a very
	% good reason for removing them, you should leave these lists in the document.  Comment them out to remove
	% them.
	\listoffigures
	\listoftables

	% The following sets up the document for the main part of the thesis or dissertation. Do not comment out or
	% remove this line.
	\mainmatter

	%now go ahead and start writing your thesis
	\chapter{Introduction} \label{ch:introduction}
		Cyber-Physical Systems (CPS) are often deployed over multiple computing platforms, where the communication among distributed computing nodes is supported by shared network buses. For example, today's automobiles consist of over 50 or even 100 computing nodes and dozens of in-vehicle communication buses. In this work, we consider CPS communicating over Controller Area Network (CAN). CAN is the most popular communication protocol in automotive, but it is also applied to many other CPS application domains such as factory and plant controls, robotics, medical devices, and avionics.

CAN brings certain advantages that fits the need of many CPS systems, such as low cost, flexibility, and built-in error handling capability. However, it also comes with a set of issues that makes the design and implementation of CPS particularly challenging. One of the most significant and relevant challenges for this work is the variable timing delays introduced into the feedback control loop, as communication among sensors, actuators, and controllers is supported by CAN. Specifically, computing nodes are unsynchronized in CAN, such that each node is sending messages according to their own clocks. This, coupled with the fact that CAN messages are arbitrated based on their ID (priority), makes the timing delay of CAN message strongly affected by the (random) queuing times of messages from other nodes. Such timing variation is known to cause significant degradation in control performance and may even lead to instability~\cite{c19}. To be effectively able to deal with these challenges, we will require two things:

a) a procedure that predicts these delays in real-time thereby reducing unpredictability and 

b) a robust control algorithm which can handle timing issues in real-time.

This work is aimed at providing a feasible way to achieve both of these requirements in a single design on a Controller Area Network.

To satisfy requirement (a), we propose a timing model that effectively predicts the delays on a CAN bus in real-time. The approach will be two-staged: Before deployment into the real world, we perform an offline analysis on the CAN trace of messages and determine the nominal timing properties of the messages involved. We then develop an online strategy that uses the results of the offline analysis to successfully predict the delay (response time) of every arriving message in real-time. The challenges in coming up with an online strategy include the variable timing (jitters) in CAN. As stated earlier, the contention for the medium access causes this variability and unpredictability which cannot be ignored. In addition, all the methods cannot accommodate the {\em blackbox}-based supply chain that is typical in industries such as automotive. More specifically, since the automotive supply chain values Intellectual Property rights and their protection, nodes are often provided as blackboxes. The timing characteristics of messages sent by other nodes on the network are unknown, including message period and configuration in the software stack that is responsible for message transmission. This provides the motivation for an online approach.

Now, we move on to requirement (b). Delay compensation is an appealing idea to counter these issues in controller design for CPS~\cite{c19}. However, the current current approaches have significant drawbacks. They either perform joint {\em offline} design on timing and control, hence inapplicable for online employment~\cite{c22,c8}; or rely on precise knowledge on the timing information of software tasks in the remote sender node, which is ill-suited for CAN where nodes are unsynchronized~\cite{c2}. 

To fully leverage the benefit of delay compensation in controller design, the control algorithm shall offer a way to embed timing related state and control constraints. Conventional controllers, such as Proportional-Integral-Derivative (PID) and Linear-Quadratic Regulator (LQR) lack this key capability to anticipate future events and adjust control actions accordingly: the control value is computed (often in closed-form) based on a fixed plant model~\cite{c2,c6,c14}. In particular, the random nature of timing delay of CAN messages will be impossible to capture in such a controller design. Therefore, this work considers Model Predictive Control (MPC), an advanced method that is becoming popular in chemical plants, process industries, and more recently, robotics and automotive~\cite{c13}. This algorithm predicts optimal control efforts for a small window (prediction horizon) into the future, while implementing the current control effort. Different from PID and LQR, it uses the current plant measurements and the current dynamic state of the process to dynamically adjust the control action. This  conveniently allows incorporation of delays in the control loop at each time step~\cite{c2,c6,c14}, which is vital in the approach towards delay compensation.


		\section{Contributions} \label{se:contributions}
% 			\subsection{A sub-section} \label{ss:this_subsection}
% 				\lipsum[1-4]
% 		\section{Another Section} \label{se:another_section}

\begin{figure}[htp]
	\centering
	\includegraphics[scale =0.65]{CANSetup.pdf}
	\caption{CAN based Control Systems with Feedback Loops with a proposed Delay Predictor for Compensation}\label{fig:CAN Setup}
\end{figure}
			This work considers the problem of online delay compensation for MPC design, and claims the following contributions: 
\begin{itemize}
\item We propose a novel approach that effectively predicts the message delay on a CAN bus in real-time. The approach contains two stages. The first is to analyze a collected set of traces (message reception events) and determine the nominal timing properties of the messages, without assuming detailed knowledge on the sender nodes. The second stage uses the results of the first stage to successfully predict the delay of every relevant message in real-time. 

\item We also ensure its real-time feasibility by implementing the method and measuring its worst-case execution time (WCET) on a low-end 8-bit microcontroller. 

\item We then include the timing information as constraints into the MPC controller design of a Cruise Control system~\cite{c20} and a DC servo~\cite{TrueTime}, to demonstrate the control performance improvement in terms of the rise time and maximum overshoot of the plant output. 

The entire setup is represented by the block diagram shown in Figure \ref{fig:CAN Setup}.
\end{itemize}

	\chapter{Review of Literature} \label{ch:lit_review}
    With respect to controller design, the literature does offer quite a few works that perform offline timing and control analysis together (e.g.,~\cite{AntaRTSS2009,c15,c16,c19,c22}). However, as the name suggests, offline methods lead to overly conservative designs, and they are unsuitable to be employed in real-time as well as in conjunction with the MPC design~\cite{c2,c13}.These are usually performed using computer simulations for all time stamps based on the employed scheduling algorithm. As far as Model predictive control is concerned, these offline computer simulations will prove too expensive. Online methods must design a schedule of events and apply the control strategy at each sampling instant for the \textit{finite} horizon. This saves a lot of computation as compared to the offline strategies as we are optimizing only for a finite window and not for all the samples.

The timing delay of CAN messages has been studied in several contexts. Since CAN is typically used to support applications with hard real-time constraints, its worst-case delay (i.e., the maximum delay from its readiness to the completion of transmission) has been studied in a variety of scenarios to make sure that the message deadline is met~\cite{c1,c18,DavisECRTS2011,LiuSIES2013}. However, this analysis only provides an upper bound on the message delay. A  stochastic analysis framework~\cite{ZengTII2009} is provided based on Markov chain model for computing CAN message delay distribution. In another work~\cite{c7}, the same authors build a probability mix model to predict the distribution of CAN message delay based on simulation data and regression techniques. \cite{c9} uses real message traces to analyze message delay distribution. However, these analysis methods~\cite{ZengTII2009,c7,c9} still only provide a statistical distribution on message delay in an offline fashion, but MPC controller design shall adjust the control action based on the actual delay of each message instance. Furthermore, they rely on strong assumptions on the CAN system, such as a particular implementation strategy in the CAN hardware controller that is often unrealistic~\cite{DiNataleETFA2013}.

Ever since MPC became popular with the chemical plants and process industries~\cite{c13}, significant research has gone into this area of control theory. However, there is still very little work which studies MPC operating over CAN with online delay analysis. \cite{c2} is probably the only exception. It has a few noticeable drawbacks. First, it depends on the availability of the delays for all software tasks in the CAN network. This information is very difficult to predict for tasks on the same node, not to mention those from other remote (blackbox) nodes. Also, it requires the apriori knowledge of the message transmission times. However in CAN, the message transmission time depends on the actual content of the message because of the possible bit stuffing. Hence, it is random and possibly difficult to accurately predict. There are also other works on MPC design over different network buses such as Ethernet~\cite{c5,c21}, but CAN is very different due to the timing characteristic (deterministic arbitration in CAN vs. random backoff mechanism in Ethernet), and the approach is not transferable.
	\chapter{Model Predictive Control} \label{ch:mpc}
    	\section{Overview} \label{se:mpc_overview}
        Simple lower order systems are often well controlled by simple PID controllers. When significant time delays are introduced in these systems, or when we have to deal with a system of a very large order, it is considered very complex for a PID controller. Such systems require the additional complexity and the predictive ability of the Model Predictive Control in a controller in order to be controlled. 

The main design objective of Model Predictive Control is to optimize the trajectory of the future manipulated variable $u$ (control effort) over a finite time window in order to optimize the future performance output $y$ of the plant.  Unlike offline procedures, the MPC algorithm performs predictions on-the-fly on the dependent variables based on the changes in the independent variables of the system. The prediction is based on the past moves, the current measurements, the current states of the variables and the constraints imposed on the output - all in an effort to hold the dependent variables as close as possible to the target (or a reference trajectory) in the near future. At each step, the algorithm uses an optimization cost function $J$ over the prediction horizon ($H$),
to calculate the optimum control efforts. This horizon keeps being shifted forward at every new step while the algorithm re-calculates the new control effort. Because of this behavior, the horizon is also referred to as the \textit{receding control horizon} in the literature.

This can be understood of in terms of a day-to-day planning activity. Say we have to prepare for an exam that is in 2 months. Initially, we optimally plan for the first 5 days and implement the first day's plan. Depending on how much of it we were able to implement and how well we did, we re-plan for the next 5 days. This process is iterative until we reach the target. At each step, we plan for the upcoming 5 steps while implementing the current step, based on the prior progress, current state, while trying to stay as close to the optimal plan. 

    	\section{Design Parameters} \label{se:mpc_design}
        The key design parameters used in tuning the control performance are as follows:
        \begin{itemize}
\item Sample Time($t_s$): This is the control interval duration and it is usually the parameter that needs to be \textit{decided} based on the application and held constant as the other parameters are tuned. 

It is important to recognize that as $t_s$ decreases, rejection of disturbances gets better at the cost of the computation effort. It is vital to find the optimal balance of performance and the computational effort needed.


\item Prediction Horizon($H$):  The prediction horizon refers to  the number of future control intervals($t_s$) that the controller will predict by evaluation while optimizing its manipulated variable(MV) at the current control interval.

Unfavorable plant characteristics combined with a small value of $H$ can generate an internally unstable controller. The remedy is to increase $H$ if not large, but, if $H$ is already large, we'll need to reconsider our $t_s$ or re think other parameters. 

\item Control Horizon($m$): The control horizon refers to the number of MV moves to be optimized at a given control interval. It lies between 1 and $H$. 

Smaller values of $m$ are preferred since it leads to fewer variables in the Quadratic Program that needs to be solved at each control interval. This implies faster computation. The MATLAB MPC Toolbox uses the KWIK algorithm to solve the QP \cite{c4}.

\end{itemize}
 
        \section{Delay Compensation in MPC Design} \label{delay_comp}
    The basic idea, as we know, in MPC is its iterative, finite-horizon optimization of the plant model and the control action. At each time $t$ it samples the current plant state, predicts the changes in the other dependent state variables using online calculation, and computes a cost minimizing control strategy (via a numerical minimization algorithm) for a (relatively short) time horizon of length $H$, i.e., $[t,t+H]$. It then implements the first step of the control strategy. For the next step, all the predictions are re-calculated over the prediction horizon, yielding a new control effort. This happens over and over again whenever the controller is triggered. 

For simplicity, we use a single-input single-output linear time-invariant  system to illustrate the MPC design, but it can easily be extended to more sophisticated control systems. The control system can be described by the equations 
\begin{equation}\label{eqn:control}
\dot{x}(t) = Ax(t) + Bu(t);
\end{equation}
\begin{equation}\label{eqn:control_2}
y(t) = Cx(t)
\end{equation}
where $u(t)$ is the control effort at time $t$, $y(t)$ is the output of the system, $x(t)$ is the plant state, and $A$, $B$, and $C$ are constant.

In the control implementation, the sensor senses the status of the plant $x(t)$ and triggers the controller. The controller uses the sensed information and generates a control signal $u(t)$ every time it is triggered. 
Since the sensed data and the control signal may need to transmit over the CAN bus, there exists a lag from the moment when the sensor measures the plant states, to the moment when the actuator executes the control signal. 

The computed control signal $u(t)$ is applied to the plant until the next sensor message triggers the controller~\cite{c2}. Thus, the control effort $u(t)$ in Equation~(\ref{eqn:control}) is a piecewise constant function: 
\begin{equation}\label{eqn:horizon}
u(t) = \mu[k],  t \in [a_{s,k} + \delta_k, a_{s,k+1} + \delta_{k+1}]
\end{equation}
where $a_{s,k}$ is the $k$-th sampling instant of the sensor, $\mu[k]$ is a constant, $\delta_k$ is the delay that exists between the time of sampling $a_{s,k}$ and the reception time $t_{u,k}$ of $\mu[k]$, i.e., $\delta_k = t_{u,k} - a_{s,k}$.

In MPC, the prediction time horizon for the $k$-th sampling is $[a_{s,k},a_{s,k}+H]$. The goal is to find a control command $u(t)$ that minimizes the difference between the predicted plant output $y(t)$ and the reference trajectory $r(t)$ over the whole prediction horizon. In the control design, a cost function models the impact of various parameters in the effort of the controller to track the reference signal $r(t)$. Weights (or penalties) are assigned to the involved variables to penalize the controller for violating the associated constraints. The designed cost function is then optimized to output the optimal control effort $u_n(t)$. A typical and suitable example cost function is~\cite{c2,c4} 
\begin{equation}\label{eqn:mpc}
\begin{aligned}
& J(x(t),u(t)) \\
= & \int_{a_{s,k}}^{a_{s,k}+H} \left \{(r(s)-y(s))^T \cdot Q_1 \cdot  (r(s)-y(s)) 
+ u(s)^T \cdot Q_2 \cdot u(s) \right \} d s\\
+ & x^T(a_{s,k} + H) \cdot Q_3 \cdot x(a_{s,k} +H) 
\end{aligned}
\end{equation}
where $Q_1$, $Q_2$ and $Q_3$ are non-negative constants (weights/penalties). 

The objective is to minimize the above cost function subject to the constraints put forth in the plant model in Equations~(\ref{eqn:control})--(\ref{eqn:horizon}). Thus the design of the MPC control requires the values for $\delta_k$ in the prediction horizon.

However, in CAN, numerous messages from various unsynchronized computing nodes contend to access the medium. Time-variant message delays are introduced, since the clocks of computing nodes are randomly drifting, and consequently the time that messages are ready for bus contention is random. Therefore, it is challenging to obtain the $\delta_k$ predictions accurately and in an online fashion, which the next chapter will address.


    \chapter{Online Delay Prediction on CAN} \label{ch:online_pred}
   The purpose of having a timing model is to explicitly define and discuss the kind of timing properties that the CAN bus can inject into the system. This will help in predicting timing faults or non-ideal timing behaviors for all the CAN messages and will help in overcoming them. In this chapter, we will discuss the CAN timing model, and the  algorithm which predicts message delays and outputs it to the MPC for delay compensation. Our approach is based on two of simple characteristics in CAN.


The first is that the CAN arbitration protocol is both {\em priority-based} and {\em non-preemptive}. Priority-based means that when multiple nodes with ready messages try to access the shared bus medium, the node sending the message with lowest identifier is always the winner of the arbitration round and is transmitted next. Hence, in CAN, each message is assigned with a unique identifier (which is referred to as unique purity), and the lower the identifier, the higher the priority. Non-preemptive means that a message being transmitted can not be preempted by higher priority messages that are made available after the transmission has started. 

In addition, we observe that sensor data are periodically sampled in control systems . Hence, in most implementations of CAN, there is a software task in the middleware layer, called {\em TxTask}, which is periodically activated and responsible for packing messages from signal data and queuing them for transmission~\cite{c8}. This is supported by the automotive AUTOSAR/OSEK standard and implemented by the commercial tools from Vector~\cite{vector}. In other words, this TxTask synchronizes the queuing of all messages from the same node. A representation of the CAN system architecture is shown in Figure~\ref{fig:sysarch}. 
%
\begin{figure}
  \centering
  \includegraphics[scale = 1]{archCAN.pdf}
  \caption{CAN System architecture.}
  \label{fig:sysarch}
\end{figure}

For the purpose of timing analysis, each periodic message stream $m_i$ is characterized by the tuple
%
\[
m_i = \{ id_i, L_i, T_i\}
\]

\noindent
where $L_i$ is the data content in bits, $id_i$ is the CAN identifier, and $T_i$ is its period. CAN messages can only contain a data content that is an integer multiple of 8 bits. The actual message length needs to account for the protocol bits in addition to the data content, including stuffed bits (34 of the 46 protocol bits are subject to stuffing) that are dependent on the data content. 

For each instance $M_{i,j}$ of the message $m_i$, the arrival time (i.e., the time the signal data is sampled and the message is ready for packing by the TxTask) is denoted as $a_{i,j}$, the queuing time (i.e., the time that the message is copied to the queue in the driver and ready for bus arbitration) is denoted as $b_{i,j}$, the start time (i.e., the time the message starts transmission) is denoted as $s_{i,j}$, and finally, the reception time (i.e., the time it finishes transmission and received by the receivers) is denoted as $t_{i,j}$. The delay of the message, defined as the difference between the message reception and the message arrival, is denoted as $d_{i,j}=t_{i,j} - a_{i,j}$.

The procedure we employ is two-staged. The first stage collects a set of traces (i.e., message reception events),  performs a trace analysis, and generates a lookup table mapping message IDs to their corresponding nominal periods (indicates possible times of arrival). This is required since the nodes may come from several suppliers and provided as blackboxes. Hence, the nodes' timing properties (like the nominal periods of the messages) are unknown before hand. Therefore, we perform the trace analysis to understand the timing behaviors of the node and attribute periods to messages. 

The second stage proposes an algorithm that stores, accesses and continuously updates the lookup table generated by the earlier stage while predicting delays at the event of reception of every message. The iterative outputs of this stage are used by the MPC for delay compensation. 

    	\section{Stage 1: Trace Analysis} \label{se:trace_analysis}
        
        The analysis starts with the collection of a sufficient amount of trace data (typically 20 minutes). We emphasize that this analysis only needs to be performed once assuming the message nominal periods remain the same.

\begin{table}[t!]
\centering
\scalebox{0.79}{
\noindent
\begin{tabular} {l|l|l|l|l|l|l}
\hline \hline
Reception  &  ID & Data Content  & Arrival & Start & Delay & REF \\\hline
1145332 &  1C1   & 03 45 03 4C & 1144241 & 1145150 & 1091 & \\ 
1148042 &  0C1   & 20 04 7C 82 20 11 BF 37 & 1147679 & 1147802 & 363 & YES \\ 
1148282 &  0C5   & 20 13 FB 69 20 0E 9F B8 & 1147679 & 1148042 & 603 & \\
1148532 &  0F9   & 00 00 40 00 00 00 03 FF & 1148201 & 1148282 & 331 & \\
1148765 &  199   & CF FF 0E 70 F1 8D 00 FF & 1148201 & 1148532 & 564 & \\
1148999 &  1E5   & 46 05 6C E0 00 FA 91 00 & 1147679 & 1148765 & 1320 & \\
1149185 &  2F9   & C8 01 0F 00 00 & 1147679 & 1148999 & 1560 & \\
1149357 &  348   & 00 00 00 00 & 1147679 & 1149185 & 1678 & \\
1149531 &  34A   & 00 00 00 00 & 1147679 & 1149357 & 1852 & \\
1151167 &  0F1   & 1C 02 00 40 & 1150810 & 1150981 & 357 & YES \\
\hline\hline
\end{tabular}
}
\vspace{1ex}
  \caption{Trace segment from a production vehicle, with computed arrival time, start time, and delay (time unit: $\mu s$).}
  \label{tab:fulltrace}
\end{table}

In the message trace, each line depicts an event of successful message reception on the bus. The trace $\mathbb{T}$ consists of a finite sequence of events $e_i$, $\mathbb{T}=\{e_0, e_1, \ldots, e_n \}$, where each event  $e_i$ is a tuple that includes the event time stamp (i.e., the time the message successfully transmitted on the bus and received by all other nodes), the CAN ID of the message to which the event is referred, and the data content which implicitly gives the corresponding total number of bits (including stuff-bits). By studying the trace $\mathbb{T}$ we can extract an estimate of information about a (possibly blackbox) node, including:
\begin{itemize}
\item The reconstruction of the possible message arrival time and delay starting from its reception time stamp;
\item The detection of true message periods;
\item The analysis of the scheduling delays of TxTask.
\end{itemize}

As an example, a segment of the trace collected from a production vehicle and the analyzed message delay are provided in Table~\ref{tab:fulltrace}, where the first three columns are the reception time stamp, message ID, and the data content. The bus speed is 500kbps, hence, the message transmission time is between 94$\mu s$ (microseconds) to 270$\mu s$. 



\begin{figure*}[t!]
  \centering
  \includegraphics[scale = 0.68]{activationtime.pdf}
  \caption{From message arrival to message reception.}
  \label{fig:msgtoreceive}
\end{figure*}

The reconstruction of the message arrival times is based on their reception times. This is done by backtracking in time and subtracting the factors that contribute to the total message delay, as illustrated in Figure~\ref{fig:msgtoreceive}. Specifically, the message delay $d_{i,j} = t_{i,j}-a_{i,j}$ of a message instance $M_{i,j}$ is composed of three components: 
\begin{itemize}
\item The scheduling delay of the TxTask, i.e., $(b_{i,j}-a_{i,j})$. At the end of the TxTask execution, the message is queued to the hardware buffer and ready for arbitration. 
\item The queuing delay $(s_{i,j}-b_{i,j})$, i.e., the time from its queuing to its start of transmission. 
\item The message transmission time $(t_{i,j}-s_{i,j})$, i.e., the time between its start of transmission to its reception.
\end{itemize}

The message transmission time is easy to calculate from the trace, as it is simply the quotient of the the number of bits in the message (including stuff-bits that are added depending on the data content) divided by th bus speed (how many bits are transmitted each second). Since CAN is non-preemptive, the start time is derived by subtracting the message transmission time from $t_{i,j}$.


The queuing delay and the TxTask scheduling delay depend on the other tasks and messages in the system, and are more challenging to find out. We observe that the queuing delay is generally unknown (as how much a message is delayed by messages of other unsynchronized nodes is unknown), except for a few {\em reference events} for which it can be neglected. Specifically, for each event in the trace, if the calculated start time does not overlap (i.e., separated by a safe distance to account for possible imprecision in the recorded time) with the reception time of the previous event, a bus idle time interval between the two events is identified. The message that is transmitted after the idle time is labeled as a \textit{reference event} with zero queuing delay.

\textit{Base rate messages}: Messages from the same node are almost always enqueued with harmonic periods by a dedicated task (on middleware). This supports the assumption we made clear earlier. Therefore, every instance of tasks with harmonic periods on a node is bound to be launched on the bus with another task, thus delaying one of them based on priority. Thus we define a base rate message as that which runs on each node with a period = ({\em gcd}) of the periods of the messages and has the highest local priority. Hence, the base rate message acts as a reference of the arrival times of all the other messages since every message is launched with one of the instances of the base-rate message. This relation between the arrival times is used for the reconstruction of arrival times of all messages as we will show.

First we describe the procedure for the message events in phase with the base rate messages: For each message $m_i$, we will denote all the trace events  as $e_{i,0}, e_{i,2}, \ldots , e_{i,n}$.  We then identify the first and last reference events in the trace, $e_{i,k}$ and $e_{i,j}$, with $j>k$ with the corresponding time stamps $t_{i,k}$ and $t_{i,j}$. See for example, Table~\ref{tab:fulltrace} where the reference events are marked with a label {\bf YES} in the last column. It is important to note that not all the trace events from $m_i$ are reference events of message receptions and therefore it is not necessarily $k=0$ or $j=n$. The time difference between these two reference events is $\Delta t_{i} = t_{i,j} - t_{i,k}$. The number of message receptions between the two reference events is $j-k$ and the actual message period is defined as
%
\begin{equation}
T_i = \frac{\Delta t_i}{j-k}
\end{equation}
In this way, we can minimize the impact of variable queuing delays for the messages (as all reference events suffer no queuing delay). The nominal period $\tilde{T}_i$ can be taken as the closest round integer (due to engineering practice). 

Consider for instance, the first reference event for a message is at time 3500$\mu s$, and the last event is at time 2993550$\mu s$ which is the 300-th message event. Hence, the average true period will be (2993550-3500)/299 = 10000.167$\mu s$. This is typically different from the message nominal period due to clock drift.  Hence, we can consider its nominal period to be the closest round integer 10000$\mu s$. 

\begin{figure*}
  \centering
  \includegraphics[width=0.7\linewidth]{activations.pdf}
\vspace{-1ex}
  \caption{Finding the reference event with the smallest TxTask scheduling delay.}\label{fig:queuing}
\vspace{-1ex}
\end{figure*}

We now discuss the steps to calculate the arrival times of the message $m_i$. {\em These steps are essentially the same for the second stage}. However, the second stage differs from the first stage in that it initially use the nominal period as the true period, and then continuously updates the true message period.

Since the reference transmission events of $m_i$ do not suffer queuing delays, the arrival time $a_{i,q}$ of reference event at $t_{i,q}$ can be estimated by subtracting its transmission time from $t_{i,q}$. However, this calculation does not take into account the scheduling delay of the middleware TxTask,  which cannot be neglected since this will be inherited by the messages as their queuing jitters. 

We notice that the arrival times must come at a periodic base due to the periodic activation of TxTask. Hence, starting from the first reference event at $t_{i,k}$, with the value of true period $T_i$ estimated, the arrival times are also estimated using $\tilde a_{i,q} = a_{i,k} + T_i \times (q-k)$ (shown as dotted lines in Figure~\ref{fig:queuing}). 

Then, we find the reference event with the smallest TxTask scheduling delay and use it as the base for arrival time calculation. More specifically, for all the {\em reference events} $e_{i,q}$ for message $m_i$, the value $\tilde l_{i,q} = a_{i,q} - \tilde a_{i,q}$ is computed, and the one $e_{i,q^*}$ with the smallest $\tilde l_{i,q}$ indicates the case suffering the smallest TxTask scheduling delay. The minimum value
\begin{equation}
J_{a^*_i} = \min_q \{ \tilde l_{i,q} \}
\end{equation}
(which shall be negative) is then applied as a correction to obtain the estimated arrival times $a_{i,q} = \tilde a_{i,q}+ J_{a^*_i}$. 

In Figure~\ref{fig:queuing}, a minimum value $J_{a^*_i}$ is obtained for the event indicated by the arrow, and the arrival times are shifted left to the solid lines, yielding the final result. As in the figure, the corrected arrival times $a_{i,q}$ are equal to $\tilde a_{i,q}$ for the event $e_{i,q^*}$, and for the other events it is $a_{i,q} < \tilde a_{i,q}$. The $(t_{i,q} - a_{i,q})$ values for all the reference events from the same node is assumed as the queuing jitter, i.e., the scheduling delay of the TxTask.

For message events that arrive at with a phase difference with respect to the base rate messages, we first find the phase associated with each of them and then assign the corresponding arrival times.  The following procedure detects the phases of such messages on a node in reference to the base rate message and then assigns arrival times based off of the results. The method is hinged on the initial assumption that the majority of the message instances are transmitted with a response time lower than the period of the base message. This assumption is valid in most cases and exceptions can be detected a posteriori and corrected.
The procedure is formally described in Algorithm~\ref{alg:phasedetect}, and an example is provided in Figure~\ref{fig:phasedetect}. An index {\tt baseIx} is assigned to each instance of the base message, starting with 0. An array {\tt mul[]} is defined to store the period ratio of all the other messages and the base rate message, one entry per message.

\begin{algorithm}
{\footnotesize
\caption{Calculate phase $\phi_i$ of $e_{i,j}$}
\label{alg:phasedetect}
\begin{algorithmic}[1]

\FOR{each trace event $e = e_{i,j}$}
   \STATE  Message m = e.GetMessage();
   \STATE  ECU c = m.GetEcu();
   \STATE  Message bm = e.GetBaseRateMessage();

   \IF{bm == m}
       \STATE c.SetBaseCount(index);
   \ELSE
       \STATE ph = c.GetBaseCount() \% m.GetMultiplier();
       \STATE m.phase\_count[ph]++;
   \ENDIF
\ENDFOR

\FOR{each message $m = m_i$}
   \STATE  m.phase = Index\_of\_max(m.phase\_count);
\ENDFOR

\end{algorithmic}
}
\end{algorithm}
\begin{figure}
  \centering
  \includegraphics[width=1\columnwidth]{phasedetect.pdf}
\vspace{-1mm}
  \caption{Detection of the message phase}
  \label{fig:phasedetect}
\vspace{-1mm}
\end{figure}
The minimum response time of message $m_i$ can be calculated as the minimum queuing delay caused due to the higher priority messages on the same node which are always queued before $m_i$, added to the transmission time of $m_i$ itself. Next, the phase of the message should is adjusted such that the analyzed response time in the trace should always be larger than its calculated minimum response time.
 An array of counters {\tt phase\_count[]} of size {\tt mul[i]} is defined for each message.
 For every instance of $m_i$ with period larger than the base period, the algorithm finds the value of {\tt baseIx} for the latest transmission of the base message and computes {\tt val=baseIx \% mul[i]}. Then, the counter in {\tt phase\_count[]} which matches {\tt val} is incremented. By the end of this, the \textit{index} of the counter array with the largest value is output as  the phase $\phi_i$ of the message. 


\begin{algorithm}
{\footnotesize
\caption{Calculate arrival time of $e_{i,j}$}
\label{alg:activation}
\begin{algorithmic}[1]

\FOR{each trace event $e = e_{i,j}$}
   \STATE  Message m = e.GetMessage();
   \STATE  ECU c = m.GetEcu();
   \STATE  Message bm = e.GetBaseRateMessage();

   \IF{bm == m}
       \STATE c.SetBaseCount(index);
       \STATE c.SetLastBaseArrival(e.GetArrivalTime());
   \ELSE
       \STATE mcnt = m.GetBaseDynCnt();
       \STATE bcnt = c.GetBaseCount();
       \STATE m.SetBaseDynCnt(mcnt + m.GetMultiplier());

       \STATE $\Delta$ = bcnt-mcnt;
       \IF {$\Delta < 0$}
           \STATE $\Delta$ = $\Delta$ + m.GetMultiplier();
       \ENDIF
       \STATE e.SetArrivalTime(c.GetLastBaseArrival() - $\Delta \times $bm.GetTPeriod();
   \ENDIF
\ENDFOR
\end{algorithmic}
}
\end{algorithm}

 In the example illustrated by Figure ~\ref{fig:phasedetect}, four transmission events return a phase count equal to 0 and only one returns a phase count of 1. Therefore, the algorithm outputs $\phi_i=0$ and assigns the corresponding arrival time by back-tracking until the first event of the base rate message with index $k \times mul + \phi_i$, where $k \in 0,1,2,\ldots $ is encountered. If such an event does not exist, i.e., the base rate transmission of index $\phi_i$ occurs after the message transmission event in question (the one that the algorithm is considering), then the arrival time of the base rate event with index $\phi_i$ is decremented by the period of the message and assigned as the message arrival time. This is formally described as Algorithm~\ref{alg:activation}. 


        \section{Stage 2: Delay Prediction} \label{se:delay_pred}
         Before the delay compensation is in effect, the system is first subject to the trace analysis described in the above subsection. This will populate a lookup table with message IDs and the nominal message periods. This step only needs to be carried out once, unless the nominal periods of messages change. 

The delay prediction stage calculates the message delay for each message $m_i$ following the same principle of Section~\ref{se:trace_analysis}, except that 
\begin{itemize}
\item the message true period shall be updated each time it encounters a reference event of the message, to adjust it according to the actual clock speed of the nodes; 
\item the reference event $q^*$ with the smallest TxTask scheduling delay shall be updated whenever a smaller $J_{a^*_i}$ is observed.
\end{itemize}

We will need to keep updating the $T_i$ values with time. As and when a new message arrives into the predictor, its nominal time period $\tilde{T}_i$ is looked up from the table. The difference between the nominal arrival time and the actual time of arrival of the message is returned as the delay $\delta_k$. These delays are aggregated and stored in an array $delaytracker$, one index per message id. Another array ($dynCount$) maintains the current count indicating the number of times each message has arrived up to that point. The look up table id updated after the occurrence of a Reference event. The aggregated delay value for every message is divided by the number of occurrences in that reference window and added to the $\tilde{T}_i$ corresponding to the $msgID$ in the look-up-table. The aggregate delay array is then reset. Thus, the predictions of arrival times get closer and closer to the accurate value, thereby reducing the delay in every feedback loop with time. 


\begin{algorithm}[t!]
{\footnotesize
\caption{Calculate and predict $\delta_k$ for every incoming message $m_{i,j}$}
\label{alg:onlinepred}
\begin{algorithmic}[1]

\FOR{each incoming message $m = m_{i,j}$}
  \STATE r = \textbf{nextReferenceInstance};
  \IF{(m == r)}
  		\STATE \textbf{updateTPs}(delaytracker,dynCount);
        \STATE \textbf{reset}.dynCount;
        \STATE \textbf{reset}.delaytracker; 
   \ENDIF
   \STATE delay=\textbf{calculateDelay}(m,count,delaytracker,dynCount);

//\textit{delay is given as the input to the MPC}
\ENDFOR

//Functions: 
\STATE \textbf{calculateDelay}:
	\newline \tab index = find(msgId == id);
	\newline \tab cnt[index] = cnt[index] +1;
	\newline \tab DynCount[index] = DynCount[index] + 1;
	\newline \tab expected\_arr = cnt[index] x $\tilde{T}$[index];
	\newline \tab $delay$ = time - expected\_arr;
	\newline \tab delayTracker(index)=delayTracker(index)+delay;
	\newline \textbf{end}
    
\STATE \textbf{updateTPs}:
\newline \tab  \textbf{for} each $\tilde{T}_i$ entry in the lookuptable \textbf{do}
	\newline \tab \tab offset=\{DelayTracker[i]/(DynCount[i]+1)\};
    \newline \tab \tab $\tilde{T}_i$ = $\tilde{T}_i$ + offset ;
	\newline \tab \textbf{end for}
\newline \textbf{end} 
%save('datalookup.mat','TPnom','msgId');

\end{algorithmic}
}
\end{algorithm}

The procedure is formalized in Algorithm~\ref{alg:onlinepred}.
In the main loop of execution, we maintain two functions: $calculateDelay$ and $updateTPs$. For every new message $msgIn$ that arrives into the predictor, we call the $calculateDelay$ function which returns the delay estimate for that particular message. We also check if the event is a reference event (message transmission with no queuing delay). If so, we call $updateTPs$ by passing the arrays $delaytracker$ and counter array $dynCount$. 

The function $calculateDelay$ looks up the message's Time Period, calculates the next prediction for the arrival time, and returns the difference as the delay. It also maintains the counters $count$ and $DynCount$ by incrementing the count in the index corresponding to the message ID. 

The function $updateTPs$, updates the $T$ column in the lookup table for all the messages. This is done by averaging the aggregated delay (stored in $delayTracker$) for each message over the number of times each message arrived in that window (stored in DynCount).The window is refreshed in the main loop soon after $updateTPs$ is called upon. 

Table \ref{tab:Example} shows the initial lookup Table (result of the offline trace analysis, say) for a toy example with 3 messages in the lookup. The experiment is conducted on a hardware which implements the online algorithm, and Table \ref{tab:Example_updated} shows the updated lookup table after 38 message events out of which 2 were reference message events. The algorithm considers the higher message ID to have the higher priority (0x02 in this case). Therefore, as it can be seen from the table, the message 0x00 acts as the base rate message for this particular ECU. The second and fourth instances of 0x02 are identified as the reference instances of the highest priority message on this ECU. Therefore, the table \ref{tab:Example_updated} shows the lookup table after two updates (once for every reference event). Figure \ref{fig: CANresults} also shows a snapshot the Arduino Serial port console. We can see the delay values being predicted for every message and the lookup updated upon detection of a reference event.

\begin{figure}[htp]
\begin{center}
  \includegraphics[clip, trim = 00cm 6.5cm 0cm 2cm,scale = 01]{CANresults.pdf}
\end{center}
\vspace{-1ex}
\caption{Serial Port Console showing online delay predictions and LookUp Updations}
\label{fig: CANresults}
\end{figure}


\begin{table}
\centering
{\footnotesize
\begin{tabular} {|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|
                 |@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|}
\hline
\textbf{Message ID} & \textbf{T}\\ \hline
0x02 &	 40\\\hline
0x01 &	 20\\\hline
0x00 &   10\\\hline

\end{tabular}
}
\vspace{-1mm}
  \caption{Initial lookup table for the Example with N = 3 }
  \label{tab:Example}
\vspace{-1mm}
\end{table}


\begin{table}
\centering
{\footnotesize
\begin{tabular} {|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|
                 |@{\hspace{5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|@{\hspace{0.5mm}}l@{\hspace{0.5mm}}|}
\hline
\textbf{Message ID} & \textbf{T}\\ \hline
0x02 &	 40.34\\\hline
0x01 &	 20.25\\\hline
0x00 &   9.63\\\hline

\end{tabular}
}
\vspace{-1mm}
  \caption{Updated lookup table for the Example with N = 3 after 2 reference events}
  \label{tab:Example_updated}
\vspace{-1mm}
\end{table}

Now we are fully equipped to solve the problem posed in Section~\ref{delay_comp} since we have the needed information for calculating $\delta_k$ within the prediction
horizon. Specifically, $\delta_k$ is the difference between the arrival of the sensor message instance $a_{s,k}$ and the reception of the actuator message $t_{u,k}$. Hence, the controller, once receiving the sensor message, can calculate the arrival time $a_{s,k}$ using the approach developed in this section. Once the actuator message is transmitted, the time stamp  $t_{u,k}$ is recorded, and the delay $\delta_k$ for the control loop is calculated as $t_{u,k}-a_{s,k}$. 
These values then will be  used for predicting the future delays.

    \chapter{Experimental Results} \label{ch:results}
    
    In this section, we present our experimental results. Specifically, we verify that our online delay prediction algorithm is in fact capable of performing the computations in real-time, as it has small runtime overhead. We also motivate the need of online delay prediction by analyzing the error incurred by clock drifts on delay calculation. We then use two example control systems to demonstrate the effectiveness of the delay compensation. 




    	\section{Measurement of Online Algorithm Timing
Overhead} \label{se:overhead}
To ensure that our online delay prediction algorithm is actually feasible, we implement it and measure the timing overhead. We use an Arduino Uno board which houses a low-end, 8-bit, 16 MHz ATmega328P controller. We used a CAN Shield to provide CAN capabilities to the microcontroller. This CAN  Shield adopts a MCP2515 CAN Bus controller with SPI interface and MCP2551 CAN transceiver~\cite{c23}. The experimental setup we used is as follows:

We use two Arduinos A and B; each connected to a different serial port on the same machine or one serial port on two different machines. A and B are connected via CAN. A sends over messages with different message IDs every 50 milliseconds on the CAN. B receives these messages and represents the delay predictor box. This means that B stores the lookup table that results from the offline trace analysis and performs the computation of the delays in real-time as and when the message is received over CAN from A, updates the lookup, etc. It is helpful to note that the delay prediction algorithm on B does not deal with growing data structures for a given system. 

\begin{figure}[htp]
    \centering
\begin{tikzpicture}
\begin{axis}[
%    title={WCET vs. Number of Messages},
    xlabel={Number of messages},
    ylabel={WCET [$\mu$s]},
    xmin=0, xmax= 8,
    ymin=0, ymax=450,
    xtick={0,1,2,3,4,5,6,7,8},
    ytick={0,50,100,150,200,250,300,350,400,450},
    ymajorgrids=true,
    grid style=dashed,
]
 
\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (0,0)(2,148)(3,220)(4,268)(5,320)(6,372)(7,420)
    };

\end{axis}
\end{tikzpicture}
\vspace{-1ex}
\caption{WCET trend with increasing number of messages.}
    \label{fig:WCET}
\end{figure}


We measure the worst-case execution time (WCET) of stage 2 (i.e., the stage after the nominal period is constructed). As in Figure~\ref{fig:WCET}, the WCET of the proposed method increases almost linearly with the number of messages. As a general trend, it can be observed that every message adds around 48$\mu$s to the WCET. This is significantly smaller than the typical period of control loop (in the range of tens of milliseconds to several seconds), and thus the algorithm is feasible for real-time delay prediction. 
		\section{Necessity of Online Delay Prediction} \label{se:necessity}
        
        We demonstrate the need of online delay prediction, due to the long operation time of the controller and the unsynchronized nature of CAN nodes. We show that the delay prediction shall  be adjusted in an online fashion because the clock drift is in general not constant and the true message period changes all the time. 

\begin{table}[htp]
\centering
{
\begin{tabular} {l|l|l}
\textbf{Trace Time (s)} & \textbf{\# Received Instances} & \textbf{True Period ($\mu s$)} \\ \hline
0-1000 &	1000153 	& 9984.71  \\ \hline
1000-2000&	1000150 	& 9985.08 \\
\end{tabular}
}
\vspace{1ex}
\caption{Evidence of clock drift on the same node: a message on a production vehicle with different number of instances received in the same length of time.}
\label{tab:clkdrift}
\end{table}

\begin{figure}[htp]
\begin{center}
  \includegraphics[scale = 1]{clkdrift.pdf}
\end{center}
\vspace{-1ex}
\caption{Clock drifts over time for two CAN nodes in the production vehicle.}
\label{fig:drift_hyb}
\end{figure}

Table~\ref{tab:clkdrift} shows the effect of the clock drift on a message of a production vehicle. The trace is recorded for more than 7000 seconds. As exemplified in the table, there are about three more message instances received in the first 1000 seconds than in the second 1000 seconds. Clearly, the constant clock speed assumption, which translates in constant message periods, cannot be sustained. Otherwise, the predicted message delay will be constantly growing, which reaches about three times of its period (period is approximately 10 ms) i.e., about 30 milliseconds, as opposed to the actual delay that is always below 1 millisecond after the system operates for 2000 seconds. 

Figure~\ref{fig:drift_hyb} shows the measured clock drifts of two different nodes in the same vehicle. The clock drift is calculated as the ratio between the message true period and its nominal period. It is evident that the clocks of different nodes are generally drifting at different speeds all the times. Hence, even if the nominal message periods are known apriori, it is still necessary to continuously compute the true periods in an online fashion to avoid accumulating the error. 


        \section{MPC Performance} \label{se:mpc_performance}
        
        Next, to simulate the control system and evaluate the performance of the controller delay compensation, we use the Model Predictive Control Toolbox on MATLAB in a CAN environment provided by TrueTime~\cite{c19}. It is a simulator software for real-time control systems that works with MATLAB and Simulink. We apply the delay prediction method to two example systems. One is a cruise control system described in~\cite{c20}, the other is a DC servo~\cite{TrueTime}. 

In the cruise control system~\cite{c20}, the vehicle mass ($m$) is set to 1000kg, the damping coefficient ($b$) is 50Ns/m, and the nominal control force is 500N. It is assumed that rolling resistance and air drag are proportional to the car's speed. The control system equation in the state-space form is
%
\begin{equation}
[\dot{v}] = \begin{bmatrix}
\dfrac{-b}{m}
\end{bmatrix} [v] + \begin{bmatrix}
\dfrac{1}{m}
\end{bmatrix} [u]
\end{equation}
\begin{equation}
y = \begin{bmatrix}
1
\end{bmatrix} [v] 
\end{equation}
where $v$ is the velocity of the car which is also the output state.

The DC Servo, found in the Truetime Library~\cite{TrueTime}, has the following state equations:
\begin{equation}
[\dot{x}] = \begin{bmatrix}
-50 & 0\\ 
1 & 0
\end{bmatrix} [x] + \begin{bmatrix}
1 & 0
\end{bmatrix} [u]
\end{equation}
\begin{equation}
[y] = \begin{bmatrix}
0 & 1000 \\ 
\end{bmatrix} [x]
\end{equation}
 




We inject delays in the system for both  the actuator message  (from controller to actuator) and the sensor message (from sensor to controller). To match the actual implementation of CAN systems, the message trace used to inject delay in the control system is recorded from a CAN bus in a production vehicle. 

\begin{figure}[htp]
    \centering
        \includegraphics[clip, trim = 01cm 2cm 3cm 0cm,width=1\columnwidth]{sample.eps}
    \caption{A sample output explaining the concepts of rise time and overshoot.}
    \label{fig:sample}
\end{figure}

Figure~\ref{fig:sample} shows a sample output of the cruise control plant for a particular MPC design with online delay prediction. The y-axis shows the plant output $y(t)$ against time $t$ on the x-axis.  The black solid line represents the plant output and the red dashed line denotes the reference trajectory. We are interested in comparing the rise time (i.e., the time needed to meet the changing trajectory),  and the maximum overshoot (the amount the control output exceeding its reference) against MPC designs for (i) no delay prediction and (ii) online prediction timing model.  

\begin{figure}[t!]
	\centering
\begin{subfigure}{.5\textwidth}    
\begin{tikzpicture}
\begin{axis}[
%    title={Rise Time ($T_r$) for various constraints on Control Effort $u(t)$},
    xlabel={+/- Constraints},
    ylabel={Rise Time [s]},
    xmin=400, xmax= 1100,
    ymin=1, ymax=3.4,
    xtick={400,500,600,700,800,900,1000,1100},
    ytick={0,0.4,0.8,1.2,1.6,2,2.4,2.8,3.2,3.6},
    legend pos=north west,
    ymajorgrids=true,
    grid style=dashed,
]
 
\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (400,3.0333)(500,2.4927)(600,2.1579)(700,1.9601)(800,1.8586)(900,1.8084)(1100,1.7905)
    };

 
 \addplot[
    color=red,
    mark=triangle,
    ]
    coordinates {
    (400,2.9253)(500,2.3844)(600,2.0283)(700,1.8107)(800,1.6879)(900,1.6229)(1100,1.5976)
    };
    \legend{No Prediction,Online Prediction}
\end{axis}
\end{tikzpicture}
 \caption{Rise Time ($T_r$) for various constraints on control effort $u(t)$.}
  \label{fig:cruise_rise}
\end{subfigure}
\linebreak
\linebreak
\begin{subfigure}{.5\textwidth} 
\begin{tikzpicture}
\begin{axis}[
%    title={Overshoot for various constraints on Control Effort $u(t)$},
    xlabel={+/- Constraints},
    ylabel={Overshoot},
    xmin=400, xmax= 1100,
    ymin=0, ymax=0.6,
    xtick={400,500,600,700,800,900,1000,1100},
    ytick={0,0.1,0.2,0.3,0.4,0.5},
    legend pos=north west,
    ymajorgrids=true,
    grid style=dashed,
]
 \addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (400,0.084)(500,0.089)(600,0.0984)(700,0.1146)(800,0.1296)(900,0.1412)(1100,0.1581)
    };
    
 \addplot[
    color=red,
    mark=triangle,
    ]
    coordinates {
    (400,0.1754)(500,0.2141)(600,0.2424)(700,0.2812)(800,0.3114)(900,0.3497)(1100,0.3987)
    };
    

    \legend{No Prediction,Online Prediction}
\end{axis}
\end{tikzpicture}
\caption{Overshoot for various constraints on control effort $u(t)$.}
  \label{fig:cruise_overshoot}
\end{subfigure}
\caption{MPC control performance for cruise control system.}
\label{fig:cruise_performance}
\end{figure}

\begin{figure}[t!]
	\centering
\begin{subfigure}{.5\textwidth}    
\begin{tikzpicture}
\begin{axis}[
%    title={Rise Time ($T_r$) for various constraints on Control Effort $u(t)$},
    xlabel={+/- Constraints},
    ylabel={Rise Time [s]},
    xmin=0.04, xmax= 0.09,
    ymin=0, ymax=2,
    xtick={0.04,0.05,0.06,0.07,0.08,0.09,0.1},
    ytick={0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2},
    legend pos=north west,
    ymajorgrids=true,
    grid style=dashed,
]
 
\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (0.04,1.8961)(0.05,1.5376)(0.06,1.2793)(0.07,1.1251)(0.08,1.052)(0.09,1.0501)
    };

 \addplot[
    color=red,
    mark=triangle,
    ]
    coordinates {
    (0.04,1.5189)(0.05,1.2516)(0.06,1.0929)(0.07,0.9564)(0.08,0.8854)(0.09,0.8835)
    };
    \legend{No Prediction,Online Prediction}
\end{axis}
\end{tikzpicture}
 \caption{Rise Time ($T_r$) for various constraints on Control Effort $u(t)$.}
  \label{fig:servo_rise}
\end{subfigure}
\linebreak 
\linebreak
\linebreak
\begin{subfigure}{.5\textwidth} 
\begin{tikzpicture}
\begin{axis}[
%    title={Overshoot for various constraints on Control Effort $u(t)$},
    xlabel={+/- Constraints},
    ylabel={Overshoot},
    xmin=0.04, xmax= 0.09,
    ymin=0, ymax=0.35,
    xtick={0.04,0.05,0.06,0.07,0.08,0.09,0.1},
    ytick={0,0.05,0.1,0.15,0.2,0.25,0.3,0.35},
    legend pos=north west,
    ymajorgrids=true,
    grid style=dashed,
]
\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (0.04,0.007)(0.05,0.014)(0.06,0.0266)(0.07,0.0418)(0.08,0.0522)(0.09,0.0524)
    }; 
\addplot[
    color=red,
    mark=triangle,
    ]
    coordinates {
    (0.04,0.117)(0.05,0.1821)(0.06,0.2268)(0.07,0.2752)(0.08,0.3286)(0.09,0.3295)
    };
    

    \legend{No Prediction,Online Prediction}
\end{axis}
\end{tikzpicture}
\caption{Overshoot for various constraints on control effort $u(t)$.}
  \label{fig:servo_overshoot}
\end{subfigure}
\caption{MPC control performance for DC servo.}
\label{fig:servo_performance}
\end{figure}


The observed rise time and overshoot are plotted in Figures~\ref{fig:cruise_performance} and~\ref{fig:servo_performance} respectively. Here the x-axis represents the applied constraints on the control effort in the MPC Design. For example, in Figure \ref{fig:cruise_performance}, 800 represents a constraint of +/- 800N on the Force in the cruise control. It is important to note that a higher value on the x-axis indicates a less constrained control design. 

From the figures ~\ref{fig:cruise_performance} and~\ref{fig:servo_performance}, we can see that having the online prediction timing model reduces the rise time at the cost of maximum overshoot. As per the application, the designer can decide on a trade-off between the rise-time and the corresponding overshoot to decide on the appropriate constraint on the control effort. The other observation is that, as we go on reducing the constraint on the the control effort, the overshoot increases, thereby decreasing the rise time.  Similar trend is observed with the DC Servo plant (plotted in Figures~\ref{fig:servo_rise} and~\ref{fig:servo_overshoot}), where the x-axis has different scales since the plant is different. 

\begin{figure}[t!]
\centering
\begin{subfigure}{1\textwidth}
  \includegraphics[scale = 0.25]{inaccurate_pred_double.eps}
  \centering
  \caption{Inaccurate prediction with worst-case delay.}
  \label{fig:inaccurate CC}
\end{subfigure}%
\linebreak
\begin{subfigure}{1\textwidth}
 \centering
  \includegraphics[scale = 0.25]{accurate_pred_double.eps}
  \caption{Prediction with our method.}
  \label{fig:accurate CC}
\end{subfigure}
\caption{MPC control output with two different prediction methods of $\delta_k$ for cruise control system.} 
\label{fig:comparisonCC}
\end{figure}

\begin{figure}[t!]
\centering
\begin{subfigure}{1\textwidth}
  \includegraphics[scale = 0.3]{servo_inaccurate_double.eps}
  \centering
  \caption{Inaccurate prediction with worst-case delay.}
  \label{fig:inaccurate servo}
\end{subfigure}%
\linebreak
\begin{subfigure}{1\textwidth}
 \centering
  \includegraphics[scale = 0.3]{servo_accurate_double.eps}
  \caption{Prediction with our method.}
  \label{fig:accurate servo}
\end{subfigure}
\caption{MPC control output with two different prediction methods of $\delta_k$ for DC servo.}
\label{fig:comparisonServo}
\end{figure}

Finally, Figures~\ref{fig:comparisonCC} and~\ref{fig:comparisonServo} respectively show a sample plant output of the Cruise Control and the DC Servo  plant while tracking a square trajectory. The y-axis shows the plant output $y(t)$ against time $t$ on the x-axis. In the subfigure (a), the output is based on inaccurate predictions of $\delta_k$ using the worst-case message delay, and the subfigure (b) shows the output with the online delay prediction. The DC Servo output is evidently better with accurate predictions than the one with inaccurate predictions. We can see the output growing in amplitude with increasing overshoots above the reference trajectory. We also confirmed that the worst case predictions have rendered the system \textit{unstable} by using the "\textbf{isstable}" command on MATLAB. The same system is held \textit{stable} and tracks the reference better when our prediction algorithm is used. The cruise control plant output is stable in both cases, but the performance and reference tracking is noticeably better in Figure~\ref{fig:accurate CC} compared to Figure~\ref{fig:inaccurate CC}. With inaccurate worst-case delay predictions, the rise time is higher than with accurate predictions. This also confirms the results presented in Figures~\ref{fig:cruise_performance} and~\ref{fig:servo_performance}.

%	\chapter{Discussion} \label{ch:discussion}
	\chapter{Conclusions} \label{ch:conclusions}
    
    This work studies the problem of online delay prediction for control systems design in CAN-based Cyber-Physical Systems. The main contribution of this work is the timing model with online delay prediction algorithm for messages scheduled on a CAN bus. It assumes no detailed knowledge on the other computing nodes. It leverages the characteristic of CAN, and finds a special kind of message reception events, to backtrack from the message reception to the message arrival and calculate the message delay. The work then demonstrates that the online delay prediction is in fact feasible to be deployed in a real-time scenario by implementing it and measuring the Worst-Case Execution Time (WCET). Such a delay prediction method is then applied to the design of the model predictive control algorithm. The work analyzes the rise time and maximum overshoot of the control output, to illustrate the benefit of the proposed method. 
    



	% This is the standard bibtex file. Do not include the .bib extension in <bib_file_name>.
	% Uncomment the following lines to include your bibliography 
	\bibliography{onlinempc}
	\bibliographystyle{plainnat}   

	% This formats the chapter name to appendix to properly define the headers
	\appendix

	% Add your appendices here.  You must leave the appendices enclosed in the appendices 
	% environment in order for the table of contents to be correct.
% 	\begin{appendices}
% 		\chapter{First Appendix} \label{app:appendix_one}
% 			\section{Section one} \label{ase:app_one_sect_1}
% 				\lipsum[1-3]
% 			\section{Section two} \label{ase:app_one_sect_2}
% 				\lipsum[1-3]
% 		\chapter{Second Appendix} \label{app:appendix_two}
% 			\lipsum[2]
% 	\end{appendices}

\end{document}


%****************************************************************************
% Below are some general suggestions for writing you dissertation
%
% 1. Label everything with a meaningful prefix so that you
%    can refer back to sections, tables, figures, equations, etc.
%    Usage \label{<prefix>:<label_name>} where some suggested
%    prefixes are:
%			ch   : Chapter
%     se   : Section
%     ss   : Subsection
%     sss  : Sub-subsection
%			app  : Appendix
%     ase  : Appendix section
%     tab  : Tables
%     fig  : Figures
%     sfig : Sub-figures
%     eq   : Equations
%
% 2. The VTthesis class provides for natbib citations.  You should 
%    in one or more *.bib bibtex files and save them in your working 
%    directory. Suppose you have two bib files: some_refs.bib and 
%    other_refs.bib.  Then your bibliography line to include them
%    will be:
%      \bibliography{some_refs,other_refs}
%    where multiple files are separated by commas.  In the body of 
%    your work you can cite your references using natbib citations.
%    Examples:
%      Citation                     Output
%      -------------------------------------------------------
%      \cite{doe_title_2016}        [18]
%      \citet{doe_title_2016}       Doe et al. [18]
%      \citet*{doe_title_2016}      Doe, Jones, and Smith [18]
%
%    For a complete list of options see
%      https://www.ctan.org/pkg/natbib?lang=en
%
% 3. Here is a sample table.  Notice that caption is centered at the top. Also
%    notice that we use booktabs formatting.  You should not use vertical lines
%    in your tables.
% 
%				\begin{table}[htb]
%					\centering
%					\caption{Approximate computation times in hh:mm:ss for full order versus reduced order models.}
%					\begin{tabular}{ccc}
%						\toprule
%						& \multicolumn{2}{c}{Computation Time}\\
%						\cmidrule(r){2-3}
%						$\overline{U}_{in}$ m/s & Full Model & ROM \\
%						\midrule
%						0.90 & 2:00:00 & 2:08:00\\
%						0.88 & 2:00:00 & 0:00:03\\
%						0.92 & 2:00:00 & 0:00:03\\
%						\midrule
%						Total & 6:00:00 & 2:08:06\\
%						\bottomrule
%					\end{tabular}
%					\label{tab:time_rom}
%				\end{table}
% 
% 4. Below are some sample figures.  Notice the caption is centered below the
%    figure.
%    a. Single centered figure:
%					\begin{figure}[htb]
%						\centering
%						\includegraphics[scale=0.5]{my_figure.eps}
%						\caption{Average outlet velocity magnitude given an average input velocity magnitude of 0.88 m/s.} 
%						\label{fig:output_rom}
%					\end{figure}
%    b. Two by two grid of figures with subcaptions
%					\begin{figure}[htb]
%						\centering
%						\begin{subfigure}[h]{0.45\textwidth}
%							\centering
%							\includegraphics[scale=0.4]{figure_1_1.eps}
%							\caption{Subcpation number one}
%							\label{sfig:first_subfig}
%						\end{subfigure}
%						\begin{subfigure}[h]{0.45\textwidth}
%							\centering
%							\includegraphics[scale=0.4]{figure_1_2.png}
%							\caption{Subcpation number two}
%							\label{sfig:second_subfig}
%						\end{subfigure}
%
%						\begin{subfigure}[h]{0.45\textwidth}
%							\centering
%							\includegraphics[scale=0.4]{figure_2_1.pdf}
%							\caption{Subcpation number three}
%							\label{sfig:third_subfig}
%						\end{subfigure}
%						\begin{subfigure}[h]{0.45\textwidth}
%							\centering
%							\includegraphics[scale=0.4]{figure_2_2.eps}
%							\caption{Subcpation number four}
%							\label{sfig:fourth_subfig}
%						\end{subfigure}
%						\caption{Here is my main caption describing the relationship between the 4 subimages}
%						\label{fig:main_figure}
%					\end{figure}
%
% 5. You should place your images and figures in a subdirectory.  You can 
%    then tell Latex to look for all your images there using:
%      \graphicspath{{figures/}}
%    Note: you must have both sets of curly braces and include the final
%          directory /
% 
%----------------------------------------------------------------------------
%
% The following is a list of definitions and packages provided by VTthesis:
%
% A. The following packages are provided by the VTthesis class:
%      amsmath, amsthm, amssymb, enumerate, natbib, hyperref, graphicx, 
%      tikz (with shapes and arrows libraries), caption, subcaption.
%      listings, verbatim
%
% B. The following theorem environments are defined by VTthesis:
%      theorem, proposition, lemma, corollary, cojecture
% 
% C. The following definition environments are defined by VTthesis
%      definition, example, remark, algorithm
%
%----------------------------------------------------------------------------
%
%  I hope this template file and the VTthesis class will keep you from having 
%  to worry about the formatting and allow to focus on the actual writing.
%  Good luck and happy writing.
%    Alan Lattimer, VT, 2016
%
%****************************************************************************





