\documentclass[11pt]{article}   	% GP.tex is 10pt font; this is 11pt
\usepackage{graphicx}				
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{arydshln}  % needed for partitioned matrices
\usepackage{color}
\usepackage{bbm}
\usepackage{upgreek}
\usepackage{subfig}
\usepackage{booktabs}
\usepackage{wrapfig}
\usepackage{lipsum}
%%%%%%%%
\DeclareCaptionSubType[alph]{figure}
\captionsetup[subfigure]{labelformat=simple,justification=centerlast}

\renewcommand\thesubfigure{(\alph{subfigure})}
\makeatletter
\renewcommand\p@subfigure{\thefigure}
\makeatother
%%%%%%%%%


% the following two define argmin and argmax. NOTE: package ams math is also needed

\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}

\newcommand{\R}{\mathbb{R}}  % the REAL symbol
\newcommand{\E}{\mathbb{E}}  % the EXPECTATION symbol

\newcommand{\x}{{\text{\bf{x}}}}

\title{}
\author{}
\date{}							% Activate to display a given date or no date

\begin{document}
\maketitle
%\section{}
%\subsection{}

\section{Research Description}
\subsection{Introduction}

Imagine an autonomous cyber-physical system (CPS) deployed into a remote area. This CPS has a robot arm used to collect rock samples. The robot arm is positioned via a sequence of unit step inputs. A proportional-integration-derivative (PID) controller is used to produce smooth robot arm movements. The PID controller gains $K_p, K_I$ and $K_d$ were tuned prior to deployment. The CPS communicates with a rear area base station over a radio link.

Now suppose a fault has occurred in the robotic arm. The exact nature of this fault is unknown, but the PID controller gains no longer produce acceptable arm movement. Fault recovery must be initiated to restore, as much as possible, the previous behavior. No redundant hardware is available so the only viable fault recovery method is to readjust the PID controller gains. Unfortunately, the CPS does not have the computational resources to compute new gain values. It does, however, have the capability of setting the gains to any value received from the base station over the radio link.

A duplicate CPS is located at the base station. Under ideal circumstances the fault could be duplicated at the base station and the new controller gains could be determined using say Zeigler-Nichols rules. Unfortunately the exact nature of the fault is unknown so this approach won't work.


The base station can send the CPS candidate gain values and then direct the CPS to apply a unit step input and report response times, steady-state errors, and so forth. In other words, the base station works with the CPS to get observation data needed for the Bayesian optimization.

Once Bayesian optimization has determined the updated PID controller gain values, the CPS can reconfigure the PID controller and recommence operation.
 
\subsection{CPS Research Focus}

Autonomous CPSs operating in remote environments must be fault-tolerant. Fault tolerance involves two tasks: fault detection and isolation (FDI) and fault recovery (FR). FDI is not addressed in this project. 

Several fault recovery methods exist. Replacing a fault system with redundant hardware is one method but in many cases is impractical because there is no room for the extra hardware. In these cases reconfiguration of the faulty hardware may be the only viable option. That said, it is difficult to know exactly how to reconfigure a faulty system if the type of fault is not known. In this research effort Bayesian optimization is used to reconfigure a PID controller. To the best of our knowledge Bayesian approaches have been used for FDI (e.g. see~\cite{nn}), but not FR, which is the focus of this research effort.

\subsection{Background}

\subsubsection{Bayesian Optimization}
This section gives a brief overview of Bayesian optimization. Further details can be found in the book by Rasmussen and Williams~\cite{gpbook} or in the review article by Shahriari et.~al~\cite{bayes}.

The problem of interest is to find the global minimum (or maximum) of some objective function $f$
\[
\x^* = \argmin_{\x\in\mathcal{X}} f(\x)
\]

\noindent where $\mathcal{X}$ is the problem's design space. In general $\mathcal{X}\subset\R^d$. This is known as a black box optimization problem whenever the analytic form of the function $f$ is unknown. It may be possible to estimate the gradient of $f$ numerically or a grid search could be used to find the minimum. But such approaches are impractical whenever sampling $f$ is costly. For example, suppose $\x$ is a set of hyperparameters for a deep neural network and $f(\x)$ is a regression output. Each $\x$ must undergo testing, which would be costly if it takes two hours to complete. 

Every sample $\x$ from $\mathcal{X}$ produces a real number evaluation $y=f(\x)$. That is, $f: \R^d\rightarrow\R$. These evaluations (sometimes called observations) are often corrupted by noise so that the observed value is
\[
y = f(\x) + \epsilon\hspace{0.49cm} ;\ \epsilon\sim\mathcal{N}(0,\sigma)
\]
\noindent and $\E\left[ y\mid f(\x)\right]  =f(\x)$.  A global minimizer should ideally need only a few samples (especially whenever drawing samples from $f$ is costly.)

\textit{Bayesian optimization} is a sequential process that can solve these kind of black box optimization problems. It creates a probabilistic model of $f$ which is much cheaper to evaluate. Observations from previous experiments serves as training data for the model.

The probabilistic model used is a \textit{Gaussian process} (GP). The GP is a probability distribution over functions where any finite collection of function values is a multi-variate Gaussian. A GP is completely described by its mean function $\mu(\x)$ and covariance function $k(\x,\x')$, also called a \textit{kernel}. 

\begin{equation}
f(\x)\sim\text{GP}(\mu(\x), k(\x,\x'))
\end{equation}


\noindent More specifically, at test points $\x_1, \ldots, \x_n$

\begin{equation}
%\left[ \begin{matrix}y_1\\ \vdots \\ y_n\end{matrix} \right] = 
\left[ \begin{matrix}f(\x_{1})\\ \vdots \\ f(\x_{n})\end{matrix} \right]  \sim \mathcal{N}\left( \left[ \begin{matrix}\mu (\x_{1})\\ \vdots \\ \mu (\x_{n})\end{matrix} \right]  ,\left[ \begin{matrix}k(\x_{1},\x_{1})&\cdots &k(\x_{1},\x_{n})\\ \vdots &&\vdots \\ k(\x_{n},\x_{1})&\cdots &k(\x_{n},\x_{n})\end{matrix} \right]  \right)  
\label{eqn:0}
\end{equation}

\vspace{2ex}
\noindent where $\mathcal{N}(\mu, C)$ denotes a multi-variate normal distribution with mean $\mu$ and covariance matrix $C$. The squared exponential kernel
\begin{equation}
k\left( x,x^{\prime }\right)  =\exp\left( -\frac{(x-x^{\prime })^{2}}{2l^{2}} \right)  
\end{equation}

\noindent is popular. Figure~\ref{fig:f1} shows five sample functions with 1000 test points from a GP with a squared exponential kernel. Other kernel functions are available. The squared exponential kernel function tends to produce overly smooth functions which may be unsuitable for some problems. In such cases a more appropriate choice might be the Mat\'ern 5/2 kernel.
%The 5/2 version is 
%\[
%k(\x,\x^{\prime })=\alpha_{0} \left( 1+\sqrt{5r^{2}(\x,\x^{\prime };\ell)} +\frac{5}{3} r^{2}(\x,\x^{\prime };\ell)\right)  \exp \left( -\sqrt{5r^{2}(\x,\x^{\prime };\ell)} \right)  
%\]
%\noindent where $r^2(x,x'; \ell) = (x-x')^T\text{diag}(\ell)^{-2}(x-x')$ and hyperparameters $\alpha_0$ and $\ell$. This function is twice differentiable.

 
\begin{figure}[!tbp]
  \centering
  \subfloat[Samples from the prior distribution.]{\includegraphics[width=0.93\textwidth]{figures/lidija22}\label{fig:f1}}
  \hfill
  \subfloat[Samples from the posterior distribution.]{\includegraphics[width=0.93\textwidth]{figures/lidija0}\label{fig:f2}}
  \caption{Shaded region depicts our belief of where the unknown objective function may lie. It reflects our level of uncertainty. Black dots in (b) represent observations. The solid line is the mean function.}
\label{f1}
\end{figure}


It is not obvious any of the functions in Figure \ref{fig:f1} could accurately model some unknown objective function. Fortunately, GPs have a nice property: they are closed under Bayesian conditions. Let $\mathcal{D}=\left[(\x_1, y_1), \ldots, (\x_n, y_n)\right]$ be a data set where $\x_i$ is the $i$-th training point and $y_i= f(\x_i)$ the observed value. Closed under Bayesian conditions means if $f\sim\text{GP}$, then $f\mid \mathcal{D}\sim \text{GP}$. However, this does require incorporating $\mathcal{D}$ into both the mean and covariance functions. %(See \cite[pg.~16]{gpbook}). 

Let $X$ and $X_*$ denote the set of training points and testing points, respectively, with $Y$ and $Y_*$ the corresponding observation sets.   The joint distribution of training observations and test point outputs is 

\begin{equation}
\left[ \begin{matrix}Y\\ Y_{* }\end{matrix} \right]  \sim \mathcal{N}\left( \begin{matrix}\mathbf{0}\  ,&\left[ \begin{matrix}K(X,X)&K(X,X_*)\\ K(X_*,X)&K(X_*,X_*)\end{matrix} \right]  \end{matrix} \right) 
\label{eqn:1} 
\end{equation}

\noindent $K(X, X_*)$ represents a covariance matrix between the training sample points and the test sample points. The other covariance submatrices are defined similarly. Note that all functions are defined with $\mu(\cdot)=0$ for convenience. 


Equation (\ref{eqn:0}) represents a prior distribution because it does not contain any information about $\mathcal{D}$. Each function in Figure \ref{fig:f1} is $f\sim\text{GP}(0,k(x,x'))$. The posterior distribution should filter out all functions that are not consistent with observed data. In other words, any function that does not associate all of the $(\x_i, y_i)$ pairs in $\mathcal{D}$ should be rejected. We can get a correct posterior distribution by redefining the mean and covariance functions. That is, 
\begin{equation}
f(x)\mid \mathcal{D}\sim\text{GP}\left(\bar{m}(\x), \bar{k}(\x)\right)
\end{equation}

\noindent with
\begin{equation}
\begin{aligned}
%\begin{matrix}
\bar{m}(\x)&=K(X_*, X) K(X,X)^{-1}Y \\ 
\bar{k}(\x)&= k(x,x^{\prime })-K(X_*, X)K(X,X)^{-1}K(X,X_*) 
%\end{matrix} 
\end{aligned}
\end{equation}

The mean function $\bar{m}(\x)$ is used to make predictions whereas the variance function $\bar{k}(\x)$ measures the confidence of the prediction. Figure \ref{f1} shows both the prior and the posterior distribution of functions. The blue regions indicate our belief or uncertainty about where the unknown objective function might exist. In the prior distribution the blue area is large because there is no information available. Conversely, the posterior distribution has much smaller uncertainty because some data is available.


\begin{figure}[htbp]
   \centerline{\includegraphics[width=11.2cm,height=10cm]{figures/acq2}}
   \caption{Using an acquisition function to choose observation points of an unknown objective function $f$. The dashed line is the objective function whereas the solid line is the mean function $\bar{m}(\x)$. This figure appeared in~\cite{bayes}.}
   \label{fig:3}
\end{figure}

The more data acquired, the more accurate the probability proxy model.  However, taking observations of $f$ is costly so the number of observations should be limited. This restriction causes an inevitable conflict between exploration and exploitation in the search. This tradeoff is handled by an \textit{acquisition function} $\alpha(\x)$ which helps guide the search. Several different acquisition functions have been proposed such as probability of improvement~\cite{pi}, expected improvement~\cite{ei} and upper confidence bound~\cite{srin}.  For example, $\alpha_{\text{UCB}}(\x)=m(\x)+\kappa\sigma(\x)$ is the upper confidence bound acquisition function where $\kappa>0$ is user-selected to balance exploitation and exploration. But regardless of which one is selected, the basic idea is the same: at each iteration of the Bayesian optimization algorithm a global optimizer finds $\hat{\x}=\argmax_{\x\in \mathcal{X}} \alpha(\x)$. $\hat{\x}$ is then the next location of an observation. This process is shown in Figure \ref{fig:3}.

A global optimizer is needed to find the maximum of $\alpha(\cdot)$. This gives rise to a question: if a global optimizer is available, why not just use it find the minimum of $f(\cdot)$?  The answer is each observation of $f$ is costly and it will take many samples to find the minimum. On the other hand, it is relatively cheap to find the maximum of $\alpha(\cdot)$. 

In summary, Bayesian optimization creates a probabilistic model of an unknown objective function $f$. It does this by fitting a GP to observed data. The mean function $\bar{m}(\x)$ from this GP is a proxy function that approximates $f$. More samples produces a more accurate model but sampling $f$ is costly. It is therefore wise to choose sample points believed to be near a global optima of $f$. An acquisition function, which trades off exploration and exploitation, figures out where to take the next sample. Optimizing a proxy function is more efficient than optimizing an objective function.

Algorithm \ref{alg:1} shows the steps in a Bayesian optimization. $\bf{x}_t$ is a sample point at time $t$ and $\bf{y}_t$ is the corresponding observed value of $f(\cdot)$. Each iteration the acquisition function $\alpha(\cdot)$ tells where to take the next sample. 
The algorithm assumes you are trying to find the minimum of an objective function $f(\x)$. Only line 8 must be changed if the maximum of $f$ is required.

\begin{algorithm}
\caption{Bayesian Optimization}
\begin{algorithmic}[1]
\label{Algorithm 1}
\State {\bf{INIT:}} Data $\mathcal{D}_0 = \{\x_0, y_0\}$
\For {$t = 1, 2, \ldots$}
\State Update GP using $\mathcal{D}_{t-1}$
\State Find $\x_{t}=\arg\max_{\x} \:  \alpha (\x)$
\State Evaluate $y_t = f(\x_t)$
\State $\mathcal{D}_t = \mathcal{D}_{t-1}\cup \{(\x_t, y_t)\}$
\EndFor \\
\Return $\x^{*}=\arg\min_{\x} \:   y(\x)$
\end{algorithmic}
\label{alg:1}
\end{algorithm}

%\begin{wrapfigure}{L}{0.5\textwidth}
%  \begin{center}
%    \includegraphics[width=0.48\textwidth]{figures/EA}
%  \end{center}
%  \caption{Block diagram of an evolutionary algorithm.}
%\end{wrapfigure}

\subsubsection{Evolutionary Algorithms}


A global optimizer is needed to find the maximum of the acquisition function $\alpha(\x)$. This function is multimodal so gradient techniques (e.g., hill-climbing) are likely to fail. Evolutionary algorithms (EAs) are stochastic search techniques. Although they are not guaranteed to find the global optima, they are quite effective at locating near optimal points even on high dimensional, multimodal functions~\cite{EAbook}.

There are different varieties of EAs (e.g., genetic algorithms, evolution strategies, etc.) but all of them generally follow the steps shown in Figure \ref{fig:6}.  Individuals (solutions) are collected into a population which evolves over time. The initial population is randomly generated. Thereafter, each generation (iteration) individuals in the current population are evaluated for fitness (solution quality). Highly fit parents are selected for reproduction (random variation) to generate offspring (new candidate solutions). In some EAs the offspring replace the parents whereas in other EAs parents and offspring compete for survival.  Survivors become the parents in the next generation. This evolutionary process of reproduction,  fitness evaluation and survival continues until a termination criteria is met,  typically after a fixed number of generations have evolved. The best fit solution in the final population is the optimal solution.

We will use an evolution strategy to find the global maxima of $\alpha(\x)$.

\begin{figure}
\centerline{\includegraphics[width=7.9cm,height=11cm]{figures/EA}}
\caption{A basic evolutionary algorithm.}
\label{fig:6}
\end{figure}


\subsubsection{PID Control}

\begin{figure}[htbp]
   \centerline{\includegraphics[width=11.2cm,height=3cm]{figures/pid}}
   \caption{A closed loop system with a PID controller.}
   \label{fig:4}
\end{figure}

The \textit{proportional-integral-derivative} (PID) controller is simple, versatile and in wide-spread use. Referring to the block diagram in Figure \ref{fig:4}, the controller output $v(t)$ is the input to the plant. In the time domain
\[
v(t)=K_{p}e(t)+K_{I}\int e(\tau )d\tau +K_{d}\frac{de(t)}{dt} 
\]
\noindent where $e(t) = r(t)-y(t)$ is an error signal. The controller transfer function is

\[
C(s)=\frac{K_{d}s^{2}+K_{p}s+K_{I}}{s} 
\]


\begin{figure}[htbp]
   \centerline{\includegraphics[width=11.2cm,height=6.3cm]{figures/step}}
   \caption{Unit step response}
   \label{fig:5}
\end{figure}

Figure \ref{fig:5} shows a typical linear system output response to a unit step input.
The purpose of the PID controller is to make a plant's step response meet specific design criteria such as a maximum settling time or steady-state error. The plant behavior is set by adjusting the PID controller gains $K_p, K_I,$ and $K_d$. The table below shows how these gains effect a plant's step response.

% Please add the following required packages to your document preamble:
% \usepackage{booktabs}
\begin{table}[htbp]\centering
\begin{tabular}{@{}ccccc@{}}
\toprule
gain & rise  & overshoot  & settling  & steady-state \\
& time & & time & error  \\ \midrule
 $K_p$ &  $\downarrow$ & $\uparrow$ &  small change &  $\downarrow$ \\
 $K_I$ &  $\downarrow$ &  $\uparrow$ & $\uparrow$ & $\downarrow$ \\
 $K_d$ & small change & $\downarrow$  & $\downarrow$ & no change \\ \bottomrule
\end{tabular}
\end{table}

The derivative term reduces transient overshoots and is effective for plants with large dead-time. Nevertheless, it can amplify high frequency sensor noise. It is therefore not uncommon to just use a PI controller---i.e., $K_d=0$. 

It is not an easy task to determine the controller gains. Manual tuning can be done although such trial-and-error techniques can be tedious. More often tuning heuristics such as the Zeigler-Nichols rules or the Rivera et.~al~\cite{riv} methods are used. Some automatic tuning tools are also now available.




\section{Evaluation/Experimentation Plan}

The research will be split into two phases: a proof of concept phase and a validation phase.

\subsection{Proof of Concept Phase}

The plant of a CPS will be simulated by a 3rd-order, type 0 transfer function $G(s)$ with a PID controller. (This system is depicted in Figure \ref{fig:4}.) Faults in this system will be simulated by adding additional (possibly complex-conjugate) poles and/or zeros to $G(s)$.  It is assumed the number and location of these poles/zeros are unknown so the PID controller gains cannot be found using conventional tuning methods. Several different scenarios will be investigated. 

There is always some underlying objective function that relates PID gain values to step response behavior. For example, the steady-state error $\theta$ to a step input is a function of the integration gain $K_I$---i.e., $\theta=f(K_I)$.  This function is unknown because the exact nature of the fault is unknown. Bayesian optimization can create a probabilistic model of $f(K_i)$ and the proxy function minimized to find a new $K_I$ gain that minimizes the steady-state error in the faulty system. The other PID controller gain values can be optimized in a similar manner.

The Bayesian optimizer software will propose new candidate $K_p, K_I$ and $K_d$ gain values which will be implemented in the physical system PID controller. A step input will produce a new step response. The overshoot, rise time, settling time or steady-state error can be extracted and constitute data observations of some underlying objective function. The Bayesian optimizer uses this data to create more accurate probabilistic models of the underlying objective function from which a new gain value can be obtained. This new gain value is sent to the CPS for evaluation. This iterative process repeats until acceptable behavior in the physical system is restored.

Full fault recovery is desired, but not guaranteed. Therefore, acceptable behaviors will be defined and the fault recovery method will be considered successful if acceptable behavior is achieved. 

All simulation will be performed using MATLAB and SIMULINK. The GP and EA programming will be done in Python.


\subsection{Validation Phase}

A ball \& beam system will be used to validate our proposed fault recovery method. A ball is placed on a beam as shown in Figure \ref{fig:3a}. The ball rolls to the left or right depending on the beam angle $\alpha$. A servo motor rotates through an angle $\theta$ which moves a lever arm up or down thus changing $\alpha$. The objective is to get the ball to stop at some desired position along the beam length $L$.  Sensors determine the ball position on the beam.

Funding is requested to purchase a commercial quality, off-the-shelf ball \& beam system. A suitable candidate system is manufactured by Acrome and is shown in Figure \ref{fig:3b}. The Acrome ball \& beam system has a built-in PID controller to drive the servo motor. The PID controller gains are user defined and will be initially set for a slightly under-damped step response. This will be the fault-free system.  Data will be obtained by applying a unit step input and then recording the response (settling time, steady-state error, etc.) A faulty system will be created by attaching small weights at the end of the beam where the lever arm is attached. Bayesian optimization determines the new PID controller gains and the new behavior is compared against the design specifications. 

It is not expected that the new gains will completely restore the original behavior. Therefore the specifications will be slightly relaxed. Fault recovery is achieved if the system meets the relaxed specifications. For example, the original specification may require zero volt steady-state error (for a unit step input). Fault recovery may be defined as successful if the faulty system has a steady-state error of say at most $\pm 20$mV. Similar relaxed specifications will be formulated for rise time and settling time.

\begin{figure}[htbp]
   \centerline{\includegraphics[width=8.5cm,height=5.4cm]{figures/bb2}}
   \caption{A ball and beam system}
   \label{fig:3a}
\end{figure}

\begin{figure}[htbp]
   \centerline{\includegraphics[width=12.2cm,height=5.4cm]{figures/BallBeam}}
   \caption{Acrome Ball \& Beam~\cite{BandB}}
   \label{fig:3b}
\end{figure}

 
\section{Project Management and Collaboration Plan}

The PI intends to employ one upper-division, undergraduate student in this research project. (No graduate research assistant is required.) 
The undergraduate research assistant has the following tasks:

\begin{itemize}
\item Construct models using MATLAB+SIMULINK. Both fault-free and faulty system models will be required.
\item Determine PID controller gains for the fault-free system
\item Collect data from the models as needed for the Bayesian Optimization
\item Become familiar with the ball \& beam system hardware and software.
\item Conduct testing during the validation phase
\item Assist in writing a peer-reviewed paper describing this research project
\end{itemize}



\section{Broader Impacts}

A full description of this research project plus all MATLAB and Python source code will be uploaded to GitHub, making it available to the general public. 

The results of this research will be documented in a peer-reviewed conference paper. The IEEE Computational Intelligence Society annually runs the Symposium Series in Computational Intelligence. Either the control and automation symposium or the evolvable systems symposium would be suitable. Funding is requested for the PI and the student to attend the conference. Although not mandatory, the student will be encouraged (and coached) to orally present the paper at the conference.

The PI believes in the full participation of women, persons with disabilities and underrepresented minorities in STEM research. Selection priority for the undergraduate research assistant will go to an upper-division electrical or computer engineering student from a historically underrepresented group. This student must have appropriate programming language skills (MATLAB/simulink) and a strong control systems background. All engineering students at Portland State University learn MATLAB and use it in several courses. The control systems background requirement is met by completing ECE 317 (\textit{Signals \& Systems III}) or equivalent with a grade of \textbf{B} or better.

The PI annually teaches a graduate level course ECE 584 \textit{Foundations of Cyber-physical Systems}. One topic in that class is the design of fault-tolerant embedded systems. The results of this research project will be incorporated into those lectures.

\end{document}

\newpage

\bibliographystyle{unsrt}
\bibliography{gp}

\end{document}  


