\section{Introduction}
Quantum walks (QWs) are dynamical tools that are used to control the motion of a quantum particle in space and time. Due to their potential applications in quantum information and physical sciences, they have been at the focus of much research work since its first introduction \cite{aharonov}. Owing to their quantum mechanical resources, i.e., quantum superposition, quantum interference, and entanglement, QWs hold the promise to develop new algorithms for computations on quantum computers \cite{deutsch1992rapid,divincenzo1995quantum,farhi1998quantum, kempe,ambainis2003quantum,computation1,computation2, chandrashekar2010discrete,venegas2012quantum}. In Physics QWs provide a versatile platform to simulate various physical phenomena, e.g., topological phases \cite{Kitagawa2010, Kitagawa2012, Kitagawa2012a, Asboth2012, Asboth2013, Asboth2014, Tarasinski2014, Asboth2015, Obuse2015, Cedzich2016BECorspnds, Groh2016, Xiao2017, Zhan2017, Sajid2019}, Anderson localization \cite{disorder1, Ahlbrecht2011, Ahlbrecht2011a, disorder2, disorder3}, Bloch Oscillations \cite{Genske2013, Cedzich2013, Arnault2020}, molecular binding \cite{Ahlbrecht2012, Lahini2012, Krapivsky2015}, and Hofstadter spectrum \cite{SajidThesis, Cedzich2020}, to name just a few. Due to their broad spectrum of applications, QWs have been realized in experiments using different physical systems, e.g., neutral atoms trapped in optical lattices \cite{karski, robens}, trapped ions on a line \cite{exp2, Xue2009, zahringer}, photons in free space \cite{broome, schreiber}, correlated photons on continuously evanescently coupled waveguides \cite{exp1}, and integrated photonics \cite{sansoni}.

QWs exhibit different features from their classical counterparts owing to their quantum mechanical resources.
In a quantum walk (QW) a quantum particle, which can exist in a superposition of several quantum states, moves to explore several paths simultaneously.
Quantum interference takes place when different trajectories cross each other. As a result the probability distribution of a QW is strikingly different from a classical random walk. For example, the variance of a quantum walk grows quadratically faster with the number of steps of the walk compared to the linear growth of a classical random walk \cite{kempe,venegas2012quantum}.
QWs have shown speedups in comparison to their classical counterparts and, hence, are useful tools to design new fast algorithms for computation on a quantum computer \cite{ambainis2003quantum,search1}, and to simulate and control certain quantum computational tasks \cite{chandrashekar2010discrete}.
This has sparked a great interest to engineer different types of QWs and investigate their properties in different settings in the context of controlling and manipulating a desired quantum state for quantum computation and simulation.
In this vein, QWs have been investigated on a large scale which has resulted in the exploration of various types of walks. This includes QWs with decoherence \cite{decoherent,Romanelli2005,Kendon,Zhang2013,albertinDecoherence2014}, QWs with time-dependent coin \cite{Banuls2006,Xue2015,Cedzich2016,panahiyan,Katayama2020}, QWs with position dependent coin \cite{Wojcik2004,Uzma2020}, and QWs with different types of phase defects \cite{zhang2014one,Zhang2014TW0,Farooq2020}.

In this paper we engineer QW protocols (one standard QW protocol and one split-step protocol) which imprint a time and spin-dependent phase shift (TSDPS) to the wave function of a quantum particle undergoing a QW on a one-dimensional (1D) lattice. We get inspiration from the previous works \cite{Xue2015,Cedzich2016,panahiyan,Katayama2020,zhang2014one,Zhang2014TW0,Farooq2020,Cedzich2019} where a desired time evolution of a QW is achieved either by manipulating the coin parameter by making it time or position dependent, or by introducing different types of spatial phase defects into the wavefunction of the quantum particle. Our engineered QWs share common features, e.g., complete revivals and partial revivals in the probability distribution of QWs which are investigated in \cite{Xue2015,Cedzich2016,Katayama2020,Farooq2020}. We use the phase factor ($\phi$) as a control knob to control and manipulate the evolution of a QW. By numerically computing the probability distribution $P(x,n)$ and the standard deviation $\sigma(n)$ of the walk, and the probability of the particle to return to its initial position $P(x=x_{\text{i}},n)$ as a function of the number of steps of the walk, we show revivals in the evolution of QWs with TSDPS.
For rational values of the phase factor, i.e., $\phi/2\pi=p/q$, where $p$ and $q$ are mutually coprime integers, periodic revivals occur in the QW driven by the standard protocol where the period of revivals depends on the denominator $q$. In this case the quantum particle takes a finite excursion in the lattice and then comes back to its initial position.
For an irrational $\phi/2\pi$, the QW shows partial revivals with irregular periods, and the walker remains localized in a small region of the lattice, hence, showing no transport. The periods of partial revivals in this case are not regular due to the incommensurability of the $\phi/2\pi$ with the step size of the QW.
In the case of a QW walk driven by the split-step protocol, partial revivals occurs for rational values of $\phi/2\pi$ with periods of revivals double of the standard QW.

In experiments, imperfections are inevitable and, hence, imprinting an exact rational or irrational $\phi/2\pi$ is challenging. From this perspective we investigate the robustness of revivals against imperfections in the imprinted phase factor. We show with our numerical results that in the presence of a linear random noise (random fluctuations) in the imprinted phase factor, signatures of revivals persist for smaller values of the noise (fluctuations) parameter in the case of $\phi/2\pi=p/q$. Increasing the strength of the noise parameter results in suppression of revivals for both rational and irrational $\phi/2\pi$. In this case the quantum particle starts to spread out showing quantum transport.

The rest of the paper is organized as follows. In Sec.~\ref{QW_TSDPS} we introduce our system, the standard QW protocol, coin operator, and shift operator with TSDPS. We present our numerical results for revivals and partial revivals in the evolution of the QW driven by the standard protocol for rational and irrational values of $\phi/2\pi$, respectively. In Sec.~\ref{SPLIT_STEP} we introduce the split-step protocol with TSDPS, and investigate revivals in the probability distribution for rational values of $\phi/2\pi$. Section~\ref{QW_TSDPS_Noise} presents the effects of random fluctuations in the phase factor on revivals in the probability distribution of the QW driven by the standard protocol. We summarize our results and conclude with a brief outlook in Sec.~\ref{conc}.

\section{Quantum Walks with a Time and Spin-dependent Phase Shift} \label{QW_TSDPS}
We consider a single quantum particle (also called a walker) with two internal degrees of freedom undergoing a QW on a 1D lattice. The internal states of the particle (also called spin states due to its analogy with a spinor particle) are represented by basis vectors $\{\ket{s}:\ s\in \{\uparrow, \downarrow\}\}$ which span a two-dimensional Hilbert space $\mathcal{H}^s$. The position states of the walker are represented by the lattice coordinates $ x$ with basis vectors $\{ \ket{x}:  x \in \mathbb{Z} \}$ spanning a 1D Hilbert space $\mathcal{H}^x$. The quantum particle undergoing a quantum walk resides in a Hilbert space that is the tensor product of the two Hilbert spaces, i.e., $\mathcal{H}^s \otimes \mathcal{H}^x$. For simplicity, we use dimensionless units by assuming the lattice constant and the time duration of a single-step of the QW to be equal to 1.

The 1D standard protocol of the QW with TSDPS consists of a set of unitary operators. We call this protocol as the ``walk operator'' which is defined as,
\begin{align}
\label{eq:_1}
\hat{W}_{\phi}(n) = \hat{S}_x(\phi n) \, \hat{C},
\end{align}
here $\hat{C}$ is known as a coin operator
and $\hat{S}_x(\phi n)$ as a spin-dependent shift operator. The coin operator acts on the internal states of the walker and rotates its spin state in the two-dimensional Hilbert space $\mathcal{H}^s$. In this work we employ the so-called Hadamard coin which is defined as,
\begin{align}
\label{eq:_coin}
\hat{C} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & \ \ \ \ 1 \\ 1 & \ \ -1 \end{pmatrix} \otimes \ket{x}\bra{x}.
\end{align}
The spin-dependent shift operator $\hat{S}_x(\phi n)$ translates the walker by one lattice site to the right or to the left depending on its internal state, and imprints a step and spin-dependent phase shift $\phi$ to the walker's wavefunction. The shift operator is defined as,
\begin{align} \label{eq:shift-x}
      \hat{S}_x (\phi n) = & \sum_{x}
           \Big[ \exp\big[i \phi n\big] \ket{\uparrow}\bra{\uparrow} \otimes \ket{x+1}\bra{x} \nonumber \\ &
                + \exp\big[ - i \phi n\big] \ket{\downarrow}\bra{\downarrow} \otimes \ket{x-1}\bra{x} \Big],
\end{align}
which shifts the walker in the spin-up (spin-down) state to the right (left). By applying the walk operator ( given in eq. (\ref{eq:_1})) to an initial state $(\ket{\Psi_{\text{i}}})$ of the walker one time constitutes a single step of the QW.
The evolution of the walk results by applying the walk operator to the initial state of the walker repeatedly for a large number of times. After certain $n$ number of steps of the walk (where $n\in \mathbb{N}$) the final state of the walker can be written as,
\begin{equation} \label{eq:final-state}
\ket{\Psi_n(\phi)}=\hat{W}_{\phi}(n) \ket{\Psi_{\text{i}}}.
\end{equation}
For a fixed value of $\phi$ the factor $n$ ensures the time dependence of the phase factor that is imprinted to the wavefunction of the walker at each step of the walk. The spatial probability distribution $P(x,n)$ is obtained by tracing out the coin degrees of freedom, i.e.,
\begin{equation} \label{eq:prob-dist}
P(x,n)=\sum_{s=0,1} \braket{\Psi_n(\phi)|\Psi_n(\phi)}.
\end{equation}
For rational values of the phase factor, i.e., $\phi/2\pi=p/q$, the walk operator $\hat{W}_{\phi}(n)$ is periodic, i.e., $\hat{W}_{\phi}(n+q*r) = \hat{W}_{\phi}(n)$ for some $r\in\mathbb{N}$. As a result the evolution of the walk (and hence, $P(x,n)$) is also periodic \cite{Cedzich2016}, i.e., $P(x,n+q*r) =P(x,n)$. To investigate revivals in the evolution of the QW we compute the return probability $P(x=x_{\text{i}},n)$ of the walker to its initial position $x_{\text{i}}$, and the standard deviation $\sigma(n)$ of the walk. After certain $n$ number of steps of the QW, $P(x=x_{\text{i}},n)$ and $\sigma(n)$ are obtained using the following expressions,
\begin{equation} \label{eq:prob-init-state}
P(x=x_{\text{i}},n)=\sum_{s=0,1} \braket{\Psi_{\text{i}}|\Psi_n(\phi)},
\end{equation}
\begin{equation} \label{eq:standard-deviation}
\sigma(n)= \sqrt{<x^2> - <x>^2}.
\end{equation}
\begin{figure*}[t]
  \centering
      \includegraphics[width=170mm]{TDPD2PiBy200S_100_49ym.eps} \\
      \includegraphics[width=170mm]{TDPD2PiBy200S_100_49AntSym1.eps} \\
      \includegraphics[width=170mm]{Simga_RtrnProb_2Piby200_100_49B.eps}
\caption{(First two rows). Probability distribution of the QW with TSDPS in the position-time plane for rational values of $\phi/2\pi$ (red, white, and blue colors indicate the maximum ($1.0$), moderate ($\sim0.1$), and minimum ($\sim0.$) probabilities, respectively).
First row: The evolution of the walk with a quantum particle initially prepared in the initial state $\Psi_{i1}$. The values of $\phi/2\pi$ are indicated in the inset.
Second row: The evolution of the walk with a quantum particle initially prepared in the initial state $\Psi_{i2}$. The values of $\phi/2\pi$ remain same in a given column.
The difference in the probability distribution for the two initial states $\Psi_{i1}$ and $\Psi_{i2}$ can be noticed by a careful inspection.
Third row: The evolution of the return probability of the walker to its initial position (solid curve), and the standard deviation (dotted curve) of the walk upto $200$ steps of the walk.
The behaviors of $P(x=x_{\text{i}}, n)$ and $\sigma(n)$ are same for both $\Psi_{i1}$ and $\Psi_{i2}$, and hence, the results are shown only for the first choice of the initial state $\Psi_{i1}$.
(a) For $\phi/2\pi=1/200$ and (b) $\phi/2\pi=1/100$ complete revivals occur after $n=q=200$ and $n=q=100$ steps of the walk, respectively. The return probability is equal to $1$ as it can be clearly seen from the corresponding subfigures (g) and (h) in the third row. (c) For $\phi/2\pi=1/49$ there occurs a partial revival when $n$ is an odd integral multiple of $q=49$ (i.e., $P(x=x_{\text{i}}, n)<1$), and a complete revival occurs when $n$ is an even integral multiple of $q=49$. The occurrence of revivals, their periodicities, and the boundedness of the standard deviation are independent of the initial state of the walker.}
\label{fig:a}
\end{figure*}
Here $<x>$ represents the average or expected value of the walker position $x$. Like $P(x,n)$, for $\phi/2\pi=p/q$ the return probability and the standard deviation of the walk are also periodic, i.e., $P(x=x_{\text{i}},n+q*r)$=$P(x=x_{\text{i}},n)$ and $\sigma(n+q*r)=\sigma(n)$. If initially a walker spreads from its initial position $x=x_{\text{i}}$, and after certain $n$ number of steps of the walk the return probability $P(x=x_{\text{i}}, n)$ becomes equal to $1$ (correspondingly $\sigma(n)$ becomes zero) we call this a complete revival. Similarly, if the walker returns to its initial position $x=x_{\text{i}}$ but with $P(x=x_{\text{i}},n)<1$ (correspondingly $\sigma(n)>0$), we call it a partial revival.
In the following sections we will investigate revivals in the probability distribution of a QW with TSDPS driven by the standard protocol (given in Eq.~(\ref{eq:_1})), and in a QW with TSDPS driven by the split-step protocol (given in Eq.~(\ref{eq:_SS_W_operator})) for various values of $\phi/2\pi$.

\subsection{Standard protocol with TSDPS for rational phase factor} \label{QW_TSDPS_Rational}
In this section we investigate the evolution of a QW driven by the standard protocol Eq.~(\ref{eq:_1}) for rational values of the phase factors, i.e., $\phi/2\pi=p/q$, where $p$ and $q$ are mutually coprime integers.
We consider a walker initially residing at the origin of a 1D lattice, and consider two possibilities for the internal state of the walker, i.e., an equal superposition of the two spin states, and a spin up state only. The two initial states of the walker can be written as,
\begin{align} \label{eq:initial-state}
\begin{aligned}
\ket{\Psi_{i1}} &=\frac{1}{2} \Big( \ket{\uparrow} + i \ket{\downarrow} \Big) \otimes \ket{0},  \\
\ket{\Psi_{i2}} & = \ket{\uparrow} \otimes \ket{0}.
\end{aligned}
\end{align}
We evolve these initial states by periodically applying the walk operator given in Eq.~(\ref{eq:_1}) for a large number of steps. To demonstrate revivals in the walk we calculate the probability distribution ($P(x,n)$) of the walk in the position-time plane, the probability of the walker to return to its initial position ($P(x=x_{\text{i}}, n)$), and the standard deviation ($\sigma(n)$) of the walk using Eqs.~(\ref{eq:prob-dist}), (\ref{eq:prob-init-state}), and (\ref{eq:standard-deviation}), respectively.

In Fig.~\ref{fig:a}, we show our numerically simulated results of the evolution of $P(x,n)$ in the position-time plane for $n=200$ steps of the walk for $\phi/2\pi=1/200$ (Fig.~\ref{fig:a}(a)), $\phi/2\pi=1/100$ (Fig.~\ref{fig:a}(b)), and $\phi/2\pi=1/49$ (Fig.~\ref{fig:a}(c)). The evolution shown in Figs.~\ref{fig:a}(a)-(c) is computed with the initial state $\psi_{i1}$ as defined in Eq.~\ref{eq:initial-state}. Revivals in the probability distribution are apparent in all three cases with period of revival equal to $n=q$ for even values of $q$ (Fig.~\ref{fig:a}(a) and (b)) and $n=2q$ for odd values of $q$ ((Fig.~\ref{fig:a}(c)). The difference in the periods of the revivals for the even and odd values of $q$ is due to the fact that after odd number of steps of the walk the probability of the walker at its initial position is zero \cite{Cedzich2016}. In the case of odd $q$ the probability distribution shows partial revivals when $n$ is an odd multiple of $q$, and complete revivals when $n$ is an even multiple of $q$. In all three cases, i.e., Fig.~\ref{fig:a}(a)-(c), the walker initially takes an excursion in a finite region of the 1D lattice when $n<q/2$, then reverts its direction after $n=q/2$ to come back to its initial position. As a whole the walker remains bounded in a finite region of the 1D lattice during its evolution.

Figures~\ref{fig:a}(d)-(f) show the evolution of the walk for the same values of $\phi/2\pi$ as in Figs.~\ref{fig:a}(a)-(c) (the value of $\phi/2\pi$ remains same in a given column) but with the second choice of the initial state, i.e., $\psi_{i2}$ defined in Eq.~\ref{eq:initial-state}. From the evolution of the probability distribution it is apparent that the occurrence of revivals and their periodicity are same as in the case of the first initial state $\psi_{i1}$. However, the spatial profile of the probability distribution $P(x,n)$ for a given number of steps $n$ is different for the two initial states. In the case of $\psi_{i1}$ the probability remains higher in position states indexed by positive integers, i.e., $x>0$ (let us call it lower branch of the $P(x,n)$) throughout the evolution of the walk. In the case of $\psi_{i2}$ higher probabilities alternate between the upper ($x<0$) and the lower ($x>0$) branches of the probability distribution, i.e., $P(x,n)$ remains higher in a single branch for certain range of $n$ and then switches to the other branch.

In order to give a clear demonstration of complete and partial revivals in the QW with TSDPS, in Figs.~\ref{fig:a}(g)-(i) we show the evolution of the return probabilities $P(x=x_{\text{i}}, n)$ of the walker to its initial position, and the standard deviation $\sigma(n)$ of the walk. Both the initial states $\psi_{i1}$ and $\psi_{i2}$ show identical behavior of $P(x=x_{\text{i}}, n)$ and $\sigma(n)$, we therefore show these results only for the first choice of the initial state, i.e., $\psi_{i1}$. The values of the phase factor remain same in a given column. For $\phi/2\pi=1/200$ and $\phi/2\pi=1/100$ complete revivals are clearly visible, i.e., $P(x=x_{\text{i}}, n)=1$ when $n$ is an integral multiple of $q$. In the case of $\phi/2\pi=1/49$, there are partial revivals when $n$ is an odd integral multiple of $q$ as $P(x=x_{\text{i}}, n)<1$, and complete revivals occur when $n$ is an even integral multiple of $q$ as $P(x=x_{\text{i}}, n)=1$. The extent of the region in which the walker remains bounded during its evolution can be estimated from the standard deviation as these are directly related. For a smaller value of the phase factor $\phi/2\pi$ (e.g., $\phi/2\pi=1/200$) the walker has a large excursion, i.e., $\sigma(n) \sim 26 $, and a smaller excursion ($\sigma(n) \sim 7 $) for a larger phase factor (e.g., $\phi/2\pi=1/49$).
In general, the TSDPS inhibits the ballistic expansion of the QW and confine the walker to a finite region of the lattice.
The extension of the bounded region decreases with the increase in $\phi/2\pi$. This behaviour is apparent from both the probability distribution and the standard deviation of the walk.
\begin{figure*}[t]
    \centering
    \includegraphics[width=57mm]{TDPD_Golden_Ratio_Sym2.eps}
    \includegraphics[width=57mm]{RtrnProb_Golden_Ratio.eps}
    \includegraphics[width=57mm]{Sigma_Golden_Ratio.eps} \\
	\includegraphics[width=57mm]{TDPD_Golden_Ratio_AntSym1.eps}
    \caption{Evolution of the QW with TSDPS and with the irrational phase factor $\phi/2\pi = (\sqrt{5}-1)/2$ (red, white, and blue colors indicate maximum ($1.$), moderate ($\sim0.1$), and minimum ($\sim0.$) probabilities, respectively). (a) and (d) show the probability distribution $P(x, n)$ of the walk in the position-time plane for $n=200$ and with symmetric (asymmetric) initial state $\Psi_{i1}$ ($\Psi_{i2}$) given in Eq.~(\ref{eq:initial-state}). (b) and (c) respectively show the long time evolution of the return probability $P(x=x_{\text{i}}, n)$ and the standard deviation of the walker with the initial state $\Psi_{i1}$ for $n=1000$. From (a) and (d) it is clear that the probability distribution is no more periodic for the irrational value of $\phi/2\pi$. The walker is completely localized in a finite region of the 1D lattice, which is also apparent from (c) which shows that $\sigma(n)$ has a finite upper bound, i.e., $\sigma<3$. The return probability in (b) shows a number of partial revivals with unpredictable periods in the long time evolution of the walk. $P(x=x_{\text{i}}, n)$ of the walker is bounded away from 0, i.e., $P(x=x_{\text{i}}, n) \ge 0.32$.
    By carefully inspecting (a) and (d) one can notice that the probability distributions for the two choices of the initial states are not exactly similar. However, the return probability and the standard deviation show similar behavior for the two choices of the initial states. We have, therefore, shown $P(x=x_{\text{i}}, n)$ and $\sigma(n)$ for the first choice of the initial state, i.e., $\Psi_{i1}$. }
\label{fig:b}
\end{figure*}
\subsection{Standard protocol with TSDPS for an irrational phase factor} \label{QW_TSDPS_Irrational}
We now consider the case of an irrational phase factor $\phi/2\pi$. The well known approximate irrational number is the Golden ratio $(\sqrt{5} - 1)/2$. The walk operator (given in Eq.~(\ref{eq:_1}) with the phase factor equal to the Golden ratio is no more periodic as $(\sqrt{5} - 1)/2$ is incommensurate with the number of steps of the walk. However, the behaviour of the probability distribution of the walker still shows a number of partial revivals in its long time evolution. In Fig.~\ref{fig:b} we show the evolution of the probability distribution $P(x,n)$ for $n=200$ for both choices of the initial states ($\Psi_{i1}$ and $\Psi_{i2}$), the evolution of the return probability of the walker with initial state $\Psi_{i1}$ to its initial position $P(x=x_{\text{i}}, n)$ for $n=1000$, and the standard deviation $\sigma(n)$ of the walker with initial state $\Psi_{i1}$ for $n=1000$. Figure~\ref{fig:b}(a) shows that the walker remains completely localized in a small region of the lattice throughout its evolution. Figure~\ref{fig:b}(b) shows a number of revivals with unpredictable periods in the long time evolution of the walk. These revivals are not strictly complete as the return probability $P(x=x_{\text{i}}, n)$ is not exactly equal to unity for all $n>1$. Note that $P(x=x_{\text{i}}, n)$ is bounded away from $0$, i.e., $P(x=x_{\text{i}}, n) \ge 0.32$ for all $n$ showing that there is a significant probability of the walker to remain at its initial position during the evolution of the walk. The reason for this behavior is that the initial state of the walker has a significant overlap with a bound state \cite{Cedzich2013}. Figure~\ref{fig:b}(c) shows the evolution of the standard deviation of the walk which has an upper bound, i.e., $\sigma(n)<3$ for all values of $n$. This shows that the walker takes a very short excursion during the evolution of the walk and remains localized in a small region of the 1D lattice. Figure~\ref{fig:b}(d) shows the same evolution of the walk as in the Fig.~\ref{fig:b}(a) but with the second choice of the initial state, i.e., $\Psi_{i2}$. By a careful inspection of Fig.~\ref{fig:b}(a) and Fig.~\ref{fig:b}(d) one can notice that the probability distributions for the two initial states are not similar. However, the return probability, the standard deviation, and the localization behaviour of the walk are same for the two initial states. We, therefore, show the return probability and the standard deviation for the first initial state, i.e., $\Psi_{i1}$ only.
\section{Split-Step Protocol with TSDPS for a rational phase factor} \label{SPLIT_STEP}
We now investigate revivals in a QW with a TSDPS driven by a split-step protocol. We consider only rational values of the phase factor, i.e., $\phi/2\pi=p/q$. The walk operator for this protocol is defined as,
\begin{align}
\label{eq:_SS_W_operator}
\hat{W}_{ss, \phi}(n) = \hat{S}_{x}^{\downarrow}(\phi n) \, \hat{C}_2 \, \hat{S}_{x}^{\uparrow}(\phi n) \, \hat{C}_1,
\end{align}
\begin{figure*}[t]
    \centering
    \includegraphics[width=71mm]{TDPD2PiBy100_Split_Step_400Steps.eps}
    \includegraphics[width=75mm]{Simga_RtrnProb_2Piby100_Split_Step_400Steps.eps} \\
    \includegraphics[width=71mm]{TDPD2PiBy49_Split_Step_400Steps.eps}
    \includegraphics[width=75mm]{Simga_RtrnProb_2Piby49_Split_Step_400Steps.eps}
\caption{Evolution of the split-step QW with TSDPS with rational phase factor $\phi/2\pi$ (red, white, and blue colors indicate maximum ($1.$), moderate ($\sim 0.1$), and minimum ($\sim 0.$) probabilities, respectively). (a) and (c) show the Probability distribution $P(x, n)$ of the walk in the position-time plane for $\phi/2\pi=1/100$ and $\phi/2\pi=1/49$, respectively.  (b) and (d) show the evolution of the return probability $P(x=x_{\text{i}}, n)$ (indicated by the blue solid curve) and $\sigma(n)$ (indicated by the orange dotted curve) corresponding to (a) and (c), respectively.
It is clear from the probability distribution that the walker remains bounded in a finite region of the lattice, as it was the case in the QW with TSDPS driven by the standard protocol. The return probability remains smaller than unity in both (b) and (d), showing partial revivals. The periods of the partial revivals are double of the standard protocol case (compare with Fig.~(\ref{fig:a})). At each period of revivals the return probability decreases and the standard deviation increases with the number of steps of the walk. The standard deviation has an upper bound due to the fact that the walker remains bounded during the evolution of the walk.}
\label{fig:Split_Step}
\end{figure*}
which consists of two coin operators $\hat{C}_1$ and $\hat{C}_2$, and two shift operators $\hat{S}_{x}^{\uparrow}(\phi n)$ and $\hat{S}_{x}^{\downarrow}(\phi n)$. We consider both the coin operators to be the Hadamard coins as defined in Eq.~(\ref{eq:_coin}).
The shift operator $\hat{S}_{x}^{\uparrow} (\phi n)$ shifts only the spin-up state of the walker to the right by a unit length and imprints a TSDPS to its wavefunction leaving the spin-down state unchanged. Mathematically, it is defined as,
\begin{align} \label{eq:shift-x-up}
     \hat{S}_{x}^{\uparrow} (\phi n) =  & \sum_{x}
             \Big[ \exp\big[i \phi n\big] \ket{\uparrow}\bra{\uparrow} \otimes \ket{x+1}\bra{x} \nonumber \\ &
              + \ket{\downarrow}\bra{\downarrow} \otimes \ket{x}\bra{x}\Big].
\end{align}
Similarly, the operator $\hat{S}_{x}^{\downarrow}(\phi n)$ shifts only the spin-down state of the walker to the left by a unit length and imprints a TSDPS to its wavefunction leaving the spin-up state unchanged. It is defined as,
\begin{align} \label{eq:shift-x-up}
     \hat{S}_{x}^{\downarrow} (\phi n) =  & \sum_{x}
             \Big[ \ket{\uparrow}\bra{\uparrow} \otimes \ket{x}\bra{x} \nonumber \\ &
              + \exp\big[- i \phi n\big] \ket{\downarrow}\bra{\downarrow} \otimes \ket{x-1}\bra{x}\Big].
\end{align}
To investigate revivals in the QW driven by the split-step protocol with TSDPS, we consider a walker in an initial state $\Psi_{i1}$ defined in Eq.~(\ref{eq:initial-state}), and evolve it through the walk operator defined in Eq.~(\ref{eq:_SS_W_operator}) for a large number of steps. We numerically compute the probability distribution $P(x,n)$, the return probability $P(x=x_{\text{i}},n)$ to the initial position, and the standard deviation $\sigma(n)$ of the walker similar to the QW driven by the standard protocol.

In Fig.~(\ref{fig:Split_Step}) we show our numerical results for $P(x,n)$, $P(x=x_{\text{i}},n)$, and $\sigma(n)$ of the split-step QW with TSDPS. The evolution is carried out with the initial state $\Psi_{i1}$ for $n=400$ steps of the walk. Figures~\ref{fig:Split_Step}(a) and (c) show $P(x,n)$ in the position-time plane for $\phi/2\pi=1/100$ and $\phi/2\pi=1/49$, respectively. The probability distributions show periodic behavior. The walker spreads initially for $n<q$, then reverts its direction around $n=q$, and comes back at its initial position $x=x_{\text{i}}$ after $n=2q$. However, the return probability $P(x=x_{\text{i}},n)$ remains smaller than unity for both even and odd $q$. This is clearly visible in Figs.~\ref{fig:Split_Step}(b) and (d) which show partial revivals in the probability distribution. The periods of revivals are double ($n=2q$ for even $q$) of the QW driven by the standard protocol (compare these results with Fig.~(\ref{fig:a}) for the corresponding values of $\phi/2\pi$). For an odd value of $q$ the period of revival is $n=2q$ which is same as it was in the case of QW with TSDPS driven by the standard protocol.

Figures.~\ref{fig:Split_Step}(b) and (d) show that at a given period of a partial revival the value of $P(x=x_{\text{i}},n)$ decreases compared to its value at a previous revival. This shows that revivals will be faint after a few periods. However, the walker will still be bounded in a finite region of the lattice.
The standard deviation $\sigma(n)$ also shows the occurrence of partial revivals as it shows quasi-periodic behaviour. At a given period of a partial revival the value of $\sigma(n)$ is larger compared to its value at a previous revival. However, it remains bounded throughout the evolution of the walk.
In general, in the split-step QW with TSDPS the walker remains bounded in a finite region of the 1D lattice. The extent of the bounded region varies inversely to the phase factor $\phi/2\pi$ similar to the standard protocol case.
\begin{figure*}[t]
  \centering
  \includegraphics[width=75mm]{RtrnProb_2PiBy49_E1_to_E5By100.eps}
  \includegraphics[width=75mm]{RtrnProb_Golden_Ratio_E1_to_5By100.eps}
      \\
  \includegraphics[width=75mm]{Sigma_2PiBy49_E1_to_E5By100.eps}
  \includegraphics[width=75mm]{Sigma_Golden_Ratio_E1_to_E5By100.eps}
  \caption{Return probability $P_{\epsilon}(x=x_{\text{i}},n)$ and standard deviation $\sigma_{\epsilon}(n)$ of the QW driven by the standard protocol for rational and irrational values of $\phi$ with $\phi_{\epsilon}$ varying randomly from one step of the walk to another. (a) $P_{\epsilon}(x=x_{\text{i}},n)$, (b) $\sigma_{\epsilon}(n)$ for $\phi=2\pi/49$ and $\epsilon = 1/20$ (solid blue curve), $\epsilon = 1/5$ (orange dashed curve), and $\epsilon = 1$ (green dotted curve). Signatures of revivals are apparent for the first two values of $\epsilon$, and in the third case the walker starts diffusion and, hence, there are no revivals. (c) $P_{\epsilon}(x=x_{\text{i}},n)$, (d) $\sigma_{\epsilon}(n)$ for $\phi=2\pi(\sqrt{5} -1)/2$ and with the same three values of $\epsilon$ (with the same color scheme) as in (a) and (b). For $\epsilon = 1/20$ the return probability (represented by the solid blue curve) is bounded away from zero which is similar to the case discussed in \ref{fig:b}(b). By increasing $\epsilon$, i.e., $\epsilon=1$ the walker starts to spread and $P_{\epsilon}(x=x_{\text{i}},n)$ is no more bounded from below. $P_{\epsilon}(x=x_{\text{i}},n)$ approaches to zero for higher values of $n$. The results shown in this figure are average values of $20$ runs.} %
\label{fig:c}
\end{figure*}
\section{Effects of Noise in the Phase Factor} \label{QW_TSDPS_Noise}
In the QW driven by the standard protocol we have shown that exact revivals in $P(x=x_{\text{i}},n)$ occur for exact rational values of $\phi/2\pi$. However, imprinting an exact rational value of $\phi/2\pi$ in experiments will be challenging as imperfections are inevitable. This can possibly eradicate any sign of revivals in the probability distribution of the QW in experiments. It is therefore worth to investigate the effects of noise in the phase factor $\phi/2\pi$ on the revivals in the probability distribution. The purpose is to find how far the above studied revivals in $P(x_{\text{i}},n)$ are robust to fluctuations in the phase factor.

We consider a phase factor $\phi_{\epsilon}(n)$ that is fluctuating randomly. Following \cite{Cedzich2016}, we model the fluctuating phase factor as,
\begin{equation} \label{eq:fluctuation_moeld}
\phi_{\epsilon}(n)=\phi + \epsilon \mathcal{R}_n,
\end{equation}
here, $\epsilon \in [0, 1]$ is the fluctuation parameter and $\mathcal{R}_n \in [-1, 1]$ is chosen randomly at each step of the walk. We investigate the evolution of the QW driven by the standard protocol with TSDPS that is fluctuating randomly from one step of the walk to another. Since we are interested to investigate the effects of random fluctuations in the phase factor on revivals, we compute only the return probability $P_{\epsilon}(x=x_{\text{i}},n)$ and the standard deviation $\sigma_{\epsilon}(n)$ of the QW.

Figure.~\ref{fig:c} shows our numerical results of the evolution of the return probability $P_{\epsilon}(x=x_{\text{i}},n)$ and the standard deviation $\sigma_{\epsilon}(n)$ of the walker undergoing a QW with TSDPS subjected to random fluctuations in the phase factor. The results shown in this figure are averaged over $20$ runs. We consider both rational and irrational values of $\phi/2\pi$. Figures~\ref{fig:c}(a) and (b) show the effect of random fluctuations on revivals for $\phi=2\pi/49$. For $\epsilon=1/20$ (red solid line) and $\epsilon=1/5$ (blue solid line) the signatures of revivals persist, while for $\epsilon=1$ (black solid line) the walker starts to spread out and, hence, revivals are absent. This clearly indicates that by increasing the fluctuation parameter from certain limit, quantum transport takes place in the QW with TSDPS as the walker is no more localized.

Figure~\ref{fig:c}(c) shows the evolution of the return probability $P_{\epsilon}(x=x_{\text{i}},n)$ and Fig.~\ref{fig:c}(d) shows the standard deviation $\sigma_{\epsilon}(n)$ of the walk for $\phi=(\sqrt{5}-1)/2$ and $\epsilon=1/20$ (red solid line), $\epsilon=1/5$ (blue solid line), and $\epsilon=1$ (black solid line). For $\epsilon=1/20$ the return probability is bounded away from zero and the walker is completely localized. This is similar to the behavior shown in Fig.~\ref{fig:b}(b). However, for $\epsilon=1/5$ and $\epsilon=1$ the return probabilities are no more bounded from below as it decreases with $n$ and approaches to zero for larger $n$. This is in stark contrast to the exact irrational $\phi$ as the walker is no more localized and spreads out from its initial position.
\section{Conclusion}\label{conc}
We have studied the evolution of 1D QWs (standard protocol and split-step protocol) with time and spin-dependent phase shifts. For the QW with standard protocol we have investigated the cases of both rational and irrational phase factors $\phi/2\pi$. For rational  values, i.e., $\phi/2\pi=p/q$ the walk operator is periodic in time and our results show revivals in the evolution of the QW. The period of revivals is $n=q$ for even values of $q$ and $n=2q$ for odd values. For the irrational case the walk operator, and hence, the return probability is no more periodic. However, we found a number of partial revivals with unpredictable periods in the evolution of the walk. The walker is completely localized in a finite region of the lattice.
For the QW with split-step protocol we considered only rational values of $\phi/2\pi$. We found partial revivals in the evolution of the walk with periods of partial revivals double (for even $q$) of the QW driven by the standard protocol.
We have shown that the results obtained for revivals, i.e., return probabilities and standard deviation do not depend on the choice of the initial state of the walker. However, the probability distributions are different for different initial states.

From the point of view of an experimental realization of revivals in QWs with TSDPS driven by the standard protocol, we further investigated the effects of random fluctuations on revivals. For both rational and irrational cases of the phase factor, we found that signatures of revivals persist for smaller values of the noise parameter. By increasing the value of the fluctuation parameter from certain limit, signatures of revivals vanish and quantum transport takes place. This so-called noise-induced transport is demonstrated for both the rational and irrational values of the phase factor. Our study is important from the point of view that under which conditions one should expect revivals in QWs with TSDPS, how far it is robust against noise in the phase factor, and when quantum transport takes over revivals. A natural extension of this work is to investigate the role of TSDPS in QWs in higher dimensions.
\section*{Acknowledgements}
Q. A., H. Q, and T. T acknowledge support from physics department, Kohat University of Science and Technology, Pakistan.

\twocolumngrid


\section{Hybrid Modeling}

\end{document}

\section{Introduction}
\label{sec:intro}

Safety is of significant importance in many modern control applications, with complex systems demanding the rigorous encoding of safety properties in the controller design process. One particularly challenging example is the safe locomotion of bipedal robots in difficult terrain; a concrete example of this is the ``stepping-stone problem'' where the feet must be placed in precise locations. In practice, the models used to design safety-critical controllers are imperfect---especially in the context of highly dynamic systems like walking robots---with model uncertainty arising due to parametric error and unmodeled dynamics. This uncertainty can lead to unsafe or dangerous behavior, thus requiring that safety-critical control be understood in the presence of model uncertainty.

In this work we consider a machine learning based approach for attenuating the impact of model uncertainty on the safe behavior of a system. The application of learning to safety-critical control of systems with uncertain models has recently exhibited significant success \citep{taylor2019learning2,cheng2019accelerating, choi2020reinforcement, lee2020learning, castaneda2020gaussian}. We look to achieve safety defined in terms of \textit{set invariance} \citep{blanchini1999set}, which is an area of active research at the intersection of data-driven methods and nonlinear control theory \citep{berkenkamp2016safe2, fisac2018general, wang2018safe, fan2019bayesian, dean2020guaranteeing}.

We will leverage control barrier functions (CBFs) for synthesizing safety controllers \citep{ames2014control, ames2017control} that achieve set invariance. CBFs have become a popular method for ensuring safety \citep{nguyen2016exponential, wang2018safe}, but they require an accurate model of the system dynamics and model uncertainty can result in a loss of safety guarantees \citep{kolathaya2018input}. Efforts have been made to make the safety guarantees of CBFs robust to the effects of model uncertainty, but in doing so these methods can cause the system's behavior to be overly conservative \citep{gurriet2018towards, xu2015robustness, taylor2019adaptive}.

One challenge in incorporating learning methods with nonlinear control is the need for diverse data capturing system behavior, and in particular the input-to-state relationship. While in simulated environments it is possible to collect high coverage data sets that accurately determine how inputs affect the evolution of the system, collecting such data on real physical systems may be prohibitively costly or damaging to the system. The lack of this data can lead to challenges with under-determination in supervised learning problems that seek to preserve the underlying structure of dynamic systems \citep{taylor2019episodic, taylor2019learning2}. Often this leads to model fits with low training loss that result in poor control performance when integrated in closed-loop controllers.

One setting for which collecting diverse data is infeasible is bipedal locomotion, as consistent walking requires well defined input profiles to ensure stability \citep{grizzle2001asymptotically}. The stepping-stone task is a historical benchmark for evaluating the safety-critical control of biped platforms \citep{nguyen20163d, nguyen2018dynamic, nguyen2020dynamic}. As compared to the statically stable motion which has been successfully employed to accomplish this task on a variety of legged robotic platforms \citep{griffin2019footstep, fankhauser2018robust}, the underactuated nature of the bipedal robot under consideration in this work requires dynamic gaits to successfully traverse a set of stepping-stones. During such dynamic motion, satisfaction of safety constraints is predicated on having accurate model information, a requirement that cannot be met on hardware platforms, thus suggesting the need for model uncertainty reduction via learning.

\textbf{Our Contribution:} The main contributions of this work are two-fold. First, we propose an episodic learning framework for iteratively reducing the impact of disturbances on the safety-critical behavior of a system. Compared to previous work, our method resolves challenges with input data diversity by recursively defining disturbance estimators and controllers that are agnostic to the underlying control-affine structure. We successfully deploy this method in the context of the classical bipedal robotic stepping-stone problem. We validate that this method is able to effectively reduce model uncertainty and improve the safe behavior of a complex system both in simulation and experiment. The second contribution is the first experimental demonstration of CBFs for safety-critical control on a bipedal robot. This challenge has been the focus of significant prior work and we demonstrate that it can be achieved through the introduction of learning.

\section{Background}
\label{sec:background}
In this section we provide a review of CBFs and projection-to-state safety (PSSf), which will be utilized in Section \ref{sec:learning} to formulate a learning approach for mitigating model uncertainty.

\subsection{Safety via Control Barrier Functions}
Consider the nonlinear control affine system given by:
\begin{equation}
    \label{eqn:eom}
    \dot{\mb{x}} = \mb{f}(\mb{x})+\mb{g}(\mb{x})\mb{u},
\end{equation}
where $\mb{x}\in\mathbb{R}^n$, $\mb{u}\in\mathbb{R}^m$, and $\mb{f}:\mathbb{R}^n\to\mathbb{R}^n$ and $\mb{g}:\mathbb{R}^n\to\mathbb{R}^{n\times m}$ are locally Lipschitz continuous on $\mathbb{R}^n$. Given a Lipschitz continuous state-feedback controller $\mb{k}:\mathbb{R}^n\to\mathbb{R}^m$, the closed-loop system dynamics are:
\begin{equation}
    \label{eqn:cloop}
    \dot{\mb{x}} = \mb{f}_{\textrm{cl}}(\mb{x}) \triangleq  \mb{f}(\mb{x})+\mb{g}(\mb{x})\mb{k}(\mb{x}).
\end{equation}
The assumption of local Lipschitz continuity of $\mb{f}$, $\mb{g}$, and $\mb{k}$ implies that $\mb{f}_\textrm{cl}$ is locally Lipschitz continuous. Thus, for any initial condition $\mb{x}_0 := \mb{x}(0) \in \mathbb{R}^n$ there exists a maximum time interval $I(\mb{x}_0) = [0, t_{\textrm{max}})$ such that $\mb{x}(t)$ is the unique solution to \eqref{eqn:cloop} on $I(\mb{x}_0)$ \citep{perko2013differential}. In the case that $\mb{f}_{\textrm{cl}}$ is forward complete, $t_{\textrm{max}}=\infty$.

The notion of safety that we consider in this paper is formalized by specifying a \textit{safe set} in the state space that the system must remain in to be considered safe. In particular, consider a set $\mathcal{C}\subset \mathbb{R}^n$ defined as the 0-superlevel set of a continuously differentiable function $h:\mathbb{R}^n\to \mathbb{R}$, yielding:
\begin{align}
    \mathcal{C} &\triangleq \left\{\mb{x} \in \mathbb{R}^n : h(\mb{x}) \geq 0\right\}, \label{eqn:safeset}
\end{align}
with $\partial\mathcal{C} \triangleq \{\mb{x} \in \mathbb{R}^n : h(\mb{x}) = 0\}$ and
    $\textrm{Int}(\mathcal{C}) \triangleq \{\mb{x} \in \mathbb{R}^n : h(\mb{x}) > 0\}$. We assume that $\mathcal{C}$ is nonempty and has no isolated points, that is, $\textrm{Int}(\mathcal{C}) \not = \emptyset \textrm{ and }\overline{\textrm{Int}(\mathcal{C})} = \mathcal{C}$. We refer to $\mathcal{C}$ as the \textit{safe set}. This construction motivates the following definitions of forward invariance and safety:

\begin{definition}[\textit{Forward Invariance \& Safety}]
A set $\mathcal{C}\subset\mathbb{R}^n$ is \textit{forward invariant} if for every $\mb{x}_0\in\mathcal{C}$, the solution $\mb{x}(t)$ to \eqref{eqn:cloop} satisfies $\mb{x}(t) \in \mathcal{C}$ for all $t \in I(\mb{x}_0)$. The system \eqref{eqn:cloop} is \textit{safe} on the set $\mathcal{C}$ if the set $\mathcal{C}$ is forward invariant.
\end{definition}

Before defining CBFs, we denote a continuous function $\alpha:[0,a)\to\mathbb{R}_+$, with $a>0$, as \textit{class $\cal{K}$} ($\alpha\in\cal{K}$) if $\alpha(0)=0$ and $\alpha$ is strictly monotonically increasing. If $a=\infty$ and $\lim_{r\to\infty}\alpha(r)=\infty$, then $\alpha$ is \textit{class $\cal{K}_\infty$} ($\alpha\in\cal{K}_\infty$). A continuous function $\alpha:(-b,a)\to\mathbb{R}$, with $a,b>0$, is said to belong to \textit{extended class $\cal{K}$} ($\alpha\in\cal{K}_e$) if $\alpha(0)=0$ and $\alpha$ is strictly monotonically increasing. If $a,b=\infty$, $\lim_{r\to\infty}\alpha(r)=\infty$, and $\lim_{r\to-\infty}\alpha(r)=-\infty$, then $\alpha$ is said to belong to \textit{extended class $\cal{K}_\infty$} ($\alpha\in\cal{K}_{\infty,e}$).

Certifying the safety of the closed-loop system \eqref{eqn:cloop} with respect to a set $\mathcal{C}$ may be impossible if the controller $\mb{k}$ was not chosen to enforce the safety of $\mathcal{C}$. CBFs can serve as a synthesis tool for attaining the forward invariance, and thus the safety of a set:
\begin{definition}[\textit{Control Barrier Function (CBF)} \citep{ames2017control}]
Let $\mathcal{C}\subset\mathbb{R}^n$ be the 0-superlevel set of a continuously differentiable function $h:\mathbb{R}^n\to\mathbb{R}$ with $0$ a regular value. The function $h$ is a \textit{control barrier function} (CBF) for \eqref{eqn:eom} on $\mathcal{C}$ if there exists $\alpha\in\mathcal{K}_{\infty,e}$ such that for all $\mb{x}\in\mathbb{R}^n$:
\begin{align}
\label{eqn:cbf}
     \sup_{\mb{u}\in\mathbb{R}^m} \dot{h}(\mb{x},\mb{u}) \triangleq \nabla h(\mb{x})\left(\mb{f}(\mb{x})+\mb{g}(\mb{x})\mb{u}\right) = L_{\mb{f}}h(\mb{x}) + L_{\mb{g}}h(\mb{x}) \u \geq-\alpha(h(\mb{x})),
\end{align}
where $L_{\mb{f}}h(\mb{x})$ and $L_{\mb{g}}h(\mb{x})$ are Lie derivatives.
\end{definition}
\noindent Given a CBF $h$ for \eqref{eqn:eom} and a corresponding $\alpha\in\cal{K}_{\infty,e}$, we define the point-wise set of all control values that satisfy \eqref{eqn:cbf}:
\begin{equation*}
    K_{\textrm{cbf}}(\mb{x}) \triangleq \left\{\mb{u}\in\mathbb{R}^m ~\left|~ \dot{h}(\mb{x},\mb{u})\geq-\alpha(h(\mb{x})) \right. \right\}.
\end{equation*}
We have the following result relating controllers taking values in  $K_{\textrm{cbf}}(\mb{x})$ to the safety of \eqref{eqn:eom} on $\mathcal{C}$:
\begin{theorem}[\cite{ames2014control}]
Given a set $\mathcal{C}\subset\mathbb{R}^n$ defined as the 0-superlevel set of a continuously differentiable function $h:\mathbb{R}^n\to\mathbb{R}$, if $h$ is a CBF for \eqref{eqn:eom} on $\mathcal{C}$, then any Lipschitz continuous controller $\mb{k}:\mathbb{R}^n\to\mathbb{R}^m$, such that $\mb{k}(\mb{x})\in K_{\textrm{cbf}}(\mb{x})$ for all $\mb{x}\in\mathbb{R}^n$, renders the system \eqref{eqn:eom} safe with respect to the set $\mathcal{C}$.
\end{theorem}
Given a nominal (but not necessarily safe) locally Lipschitz continuous controller $\mb{k}_d:\mathbb{R}^n\to\mathbb{R}^m$, a controller taking values in the set $K_\textrm{cbf}(\mb{x})$ is the safety-critical CBF-QP:
\begin{align}
\label{eqn:CBF-QP}
\tag{CBF-QP}
\mb{k}(\mb{x}) =  \,\,\underset{\mb{u} \in \mathbb{R}^m}{\argmin}  &  \quad \frac{1}{2} \| \mb{u} -\mb{k}_d(\mb{x})\|_2^2  \\
\mathrm{s.t.} \quad & \quad \dot{h}(\mb{x},\mb{u})
\geq -\alpha(h(\mb{x})). \nonumber
\end{align}

\subsection{Model Uncertainty and Projection-to-State Safety}
\label{sec:uncertainty}

In practice, the system dynamics \eqref{eqn:eom} are not known during  control design  due to parametric error and unmodeled dynamics. Instead, a nominal model of the system is utilized:
\begin{equation}
\label{eqn:eomnominal}
    \widehat{\dot{\mb{x}}} = \widehat{\mb{f}}(\mb{x}) + \widehat{\mb{g}}(\mb{x})\mb{u},
\end{equation}
where $\widehat{\mb{f}}:\mathbb{R}^n\to\mathbb{R}^n$ and $\widehat{\mb{g}}:\mathbb{R}^n\to\mathbb{R}^{n\times m}$ are assumed to be Lipschitz continuous on $\mathbb{R}^n$. By adding and subtracting the right hand side of \eqref{eqn:eomnominal} to \eqref{eqn:eom}, the dynamics of the system are:
\begin{equation}
    \dot{\mb{x}} = \widehat{\mb{f}}(\mb{x}) + \widehat{\mb{g}}(\mb{x})\mb{u} + \overbrace{\underbrace{\mb{f}(\mb{x})-\widehat{\mb{f}}(\mb{x})}_{\mb{b}(\mb{x})}+\underbrace{\left(\mb{g}(\mb{x})-\widehat{\mb{g}}(\mb{x})\right)}_{\mb{A}(\mb{x})}\mb{u}}^{\mb{d}(\mb{x},\u)},
\end{equation}
where the unknown disturbance $\mb{d}(\mb{x},\u)=\mb{b}(\mb{x})+\mb{A}(\mb{x})\mb{u}$ is assumed to be time invariant, but depends on the state and input to the system. If the function $h:\mathbb{R}^n\to\mathbb{R}$ is a CBF for the nominal model \eqref{eqn:eomnominal} on $\mathcal{C}$, the uncertainty in the dynamics directly manifests in the time derivative of $h$:
\begin{align}
\label{eqn:hdotdecomp}
    \dot{h}(\mb{x},\mb{u})  =  \underbrace{\nabla h(\mb{x})(\widehat{\mb{f}}(\mb{x})+\widehat{\mb{g}}(\mb{x})\mb{u})}_{\widehat{\dot{h}}(\mb{x},\mb{u})}
    + \underbrace{\nabla h(\mb{x})\mb{b}(\mb{x})}_{b(\mb{x})} + \underbrace{\nabla h(\mb{x})\mb{A}(\mb{x})}_{\mb{a}(\mb{x})^\top}\mb{u}.
\end{align}
Given that $h$ is a CBF for \eqref{eqn:eomnominal} on $\mathcal{C}$, let $\mb{k}:\mathbb{R}^n\to\mathbb{R}^m$ be a Lipschitz continuous state-feedback controller such that $\widehat{\dot{h}}(\mb{x},\mb{k}(\mb{x})) \geq -\alpha(h(\mb{x}))$. Defining the \textit{projected disturbance} as:
\begin{equation}
\label{eqn:projdist}
\delta(\mb{x}) \triangleq \dot{h}(\mb{x},\mb{k}(\mb{x})) - \widehat{\dot{h}}(\mb{x},\mb{k}(\mb{x})) = b(\mb{x})+\mb{a}(\mb{x})^\top\mb{k}(\mb{x}),
\end{equation}
yields:
\begin{equation}
\label{eqn:hdottrueineq}
    \dot{h}(\mb{x},\mb{k}(\mb{x})) \geq -\alpha(h(\mb{x})) - \delta(\mb{x})
\end{equation}
Assuming that $\delta$ is essentially bounded in time (there exists $M\in\mathbb{R}$, $M>0$, such that $\vert\delta\vert_\infty  \triangleq\esssup_{t\geq 0}\vert\delta(\mb{x}(t))\vert< M$), we may make use of the following definition:
\begin{definition}[\textit{Projection-to-State Safety (PSSf)} (\cite{taylor2020control})]
Given a feedback controller $\mb{k}$, the closed-loop system \eqref{eqn:cloop}, $\dot{\mb{x}} = \mb{f}_{\textrm{cl}}(\mb{x}) = \widehat{\mb{f}}(\mb{x}) + \widehat{\mb{g}}(\mb{x})\mb{k}(\mb{x}) + \mb{d}(\mb{x})$ with $\mb{d}(\mb{x})=\mb{b}(\mb{x})+\mb{A}(\mb{x})\mb{k}(\mb{x})$,  is \emph{projection-to-state safe} (PSSf) on $\mathcal{C}$ with respect to the function $h$ and projected disturbances $\delta:\mathbb{R}^n\to\mathbb{R}$ if there exists $\overline{\delta}>0$ and $\gamma\in\mathcal{K}_\infty$ such that the set $\mathcal{C}_{\delta}\supset\mathcal{C}$,
\begin{align}\label{eq:delta_safe_set}
    \mathcal{C}_{\delta} &\triangleq \left\{\mb{x} \in \mathbb{R}^n : h(\mb{x})+\gamma(\vert\delta\vert_\infty) \geq 0\right\},
\end{align}
is forward invariant for all $\delta$ satisfying $\vert\delta\vert_\infty\leq\overline{\delta}$.
\end{definition}

PSSf captures the fact that in the presence of model uncertainty, satisfying the CBF condition \eqref{eqn:cbf} for the estimated time derivative $\widehat{\dot{h}}$ is not sufficient for safety, as the projected disturbance $\delta$ appears in the lower bound on true time derivative of $h$ as in \eqref{eqn:hdottrueineq}. This results in a larger forward invariant set, given by $\mathcal{C}_\delta$, that grows with the magnitude of the projected disturbance.

\begin{figure}
    \centering
    \includegraphics{Figures/PssfCBF.pdf}
    \caption{Geometric visualization of Projection-to-State Safety. The system is able to leave the safe set $\mathcal{C}$, but remains within the larger set $\mathcal{C}_\delta$. Maximum possible deviation from the safe set grows larger with $\vert\delta\vert_\infty$.}
    \label{fig:pssf}
    \vspace{-3mm}
\end{figure}

\section{Learning Projected Disturbances}
\label{sec:learning}
In this section we explore how an estimate of the projected disturbance $\delta$ can be learned episodically from data and incorporated into control synthesis to improve PSSf behavior.

As seen in Section \ref{sec:background}, the projected disturbance $\delta$ appears in the time derivative of the barrier function $\dot{h}$, and potentially leads to unsafe behavior since it compromises the CBF condition as in \eqref{eqn:hdottrueineq}. If an upper bound $\overline{\delta}$ on $\vert\delta\vert_\infty$ is known (or determined heuristically), it could be directly incorporated into the inequality enforced in the controller:
\begin{align}
    \label{eqn:dumbreject}
    \tag{$\overline{\delta}$-CBF-QP}
    \mb{k}(\mb{x}) =  \,\,\underset{\mb{u} \in \mathbb{R}^m}{\argmin}  &  \quad \frac{1}{2} \| \mb{u}-\mb{k}_d(\mb{x}) \|_2^2  \\

    \mathrm{s.t.} \quad & \quad \widehat{\dot{h}}(\mb{x}, \u) - \overline{\delta}
    \geq -\alpha(h(\mb{x})). \nonumber
\end{align}
While this will enforce safety of the original set $\mathcal{C}$, it can be exceedingly conservative if $\overline{\delta}$ is larger than the actual projected disturbance. Furthermore, as the projected disturbance is a function of the state, its magnitude (and possibly sign) may change along a trajectory, leading to additional conservativeness in this approach.

Instead, we consider a learning approach to resolve the impact of $\delta$. To motivate such an approach, consider the following setting: in an experiment, the system is allowed to evolve forward in time from a particular initial condition and under a given state-feedback controller. During this experiment, data is collected which provides a discrete-time history of the CBF, $h$. This time history is smoothed and numerically differentiated to compute an approximate time history of the true value of the time derivative of the CBF, $\dot{h}$. This yields a collection of input-output pairs:
\begin{equation}
    {D}_i = ((\mb{x}_i, \mb{k}(\mb{x}_i)), \dot{h}_i) \in (\mathbb{R}^n \times \mathbb{R}^m) \times \mathbb{R}
\end{equation}
whereby a dataset $\mathfrak{D}=\{{D}_i\}_{i=1}^N$ can be constructed. Given a nonlinear function class $\mathcal{H}: \mathbb{R}^n \to \mathbb{R}$ and a loss function $\mathcal{L}: \mathbb{R} \times \mathbb{R} \to \mathbb{R} $, a learning problem can be specified as finding a function $\widehat{\delta} \in \mathcal{H}$ to estimate $\delta$ via empirical risk minimization:
\begin{equation}
    \label{eqn:erm}
    \tag{ERM}
    \inf_{\widehat{\delta} \in \mathcal{H}}\frac{1}{N}\sum_{i=1}^N\cal{L}\left(\widehat{\dot{h}}(\mb{x}_i,\mb{k}(\mb{x}_i))+ \widehat{\delta}(\mb{x}_i),\dot{h}_i\right).
\end{equation}
A controller can be synthesized which incorporates $\widehat{\delta}$ as follows:
\begin{align}
    \label{eqn:delta-CBF-QP}
    \tag{$\widehat{\delta}$-CBF-QP}
    \mb{k}(\mb{x}) =  \,\,\underset{\mb{u} \in \mathbb{R}^m}{\argmin}  &  \quad \frac{1}{2} \| \mb{u}-\mb{k}_d(\mb{x}) \|_2^2  \\

    \mathrm{s.t.} \quad & \quad \widehat{\dot{h}}(\mb{x}, \u) + \widehat{\delta}(\mb{x})
    \geq -\alpha(h(\mb{x})). \nonumber
\end{align}
Note that compared with \ref{eqn:CBF-QP}, the extended safe set with \ref{eqn:delta-CBF-QP} shrinks from \eqref{eq:delta_safe_set} to  $\mathcal{C}_{\delta} = \left\{\mb{x} \in \mathbb{R}^n : h(\mb{x})+\gamma(\vert\delta-\hat{\delta}\vert_\infty) \geq 0\right\}$.

We directly build upon the episodic learning framework from \citep{taylor2019episodic, taylor2019learning2} by seeking to learn $\delta$. Our approach is outlined in Algorithm 1.  In each episode, the  algorithm runs the current controller to collect data, learns a new $\widehat{\delta}$ using the newly collected data, and synthesizes a new controller. In this prior work, which was applied to less complex dynamical systems, the collected data was rich enough to determine a control affine structure. In many contexts, such as bipedal robots, such a degree of diversity is infeasible without damaging the system. We instead directly learn $\delta$ as a function of the previous controller $\mb{k}$ via a recursive relationship, as updating the estimator leads to the definition of a new projected disturbance $\delta' = b(\mb{x})+\mb{a}(\mb{x})^\top\mb{k}'(\mb{x})$. This yields a projected disturbance $\delta$  learned iteratively by modifying $\widehat{\delta}$ over the course of multiple episodes. This episodic approach to safety-critical control is captured in Algorithm \ref{alg:episodic_learning}.

\begin{algorithm}[H]
    \SetAlgoLined
    \caption{Projected Disturbance Learning (PDL)}
    \SetKwFunction{Sample}{sample}\SetKwFunction{Experiment}{experiment} \SetKwFunction{ERM}{ERM}\SetKwFunction{Augment}{augment}
    \SetKwFunction{LCBFQP}{$\widehat{\delta}$-CBF-QP}
    \SetKwInOut{Input}{input} \SetKwInOut{Output}{output}
    \label{alg:episodic_learning}
    \Input{CBF $h$, CBF derivative estimate $\widehat{\dot{h}}$, model class $\mathcal{H}$, loss function $\mathcal{L}$, nominal state-feedback controller $\mb{k}_0$, number of episodes $T$, initial condition $\mb{x}_0$}
    \Output{Augmented Controller $\mb{k}_T$}
    \BlankLine
    \For{$j = 1, \ldots, T $}{

        $\mathfrak{D}_j \leftarrow$\Experiment{$\mb{x}_0,\mb{k}_{j-1}$}  \hspace{2.438cm}\tcp*[h]{Execute experiment}\;\\
        $\widehat{\delta}\leftarrow$\ERM{$\cal{H}, \cal{L},\mathfrak{D}_j,\widehat{\dot{h}}_0$}  \hspace{3.7cm}\tcp*[h]{Fit estimator}\;\\
        $\widehat{\dot{h}}_j\leftarrow \widehat{\dot{h}}_0+\widehat{\delta}$  \hspace{5.49cm}\tcp*[h]{Update derivative estimator}\;\\
        $\mb{k}_j\leftarrow$ \LCBFQP{$\widehat{\dot{h}}_j$} \hspace{3.9cm}\tcp*[h]{Synthesize new controller}\;
    }
\end{algorithm}

\section{Bipedal Robotics}
\label{sec:bipeds}
In this section we specify the notion of learning projected disturbances to the setting of bipedal locomotion. We briefly introduce the theory of bipedal locomotion and then describe the barrier function formulations which allow us to achieve safe bipedal locomotion across stepping-stones. A deeper exploration of this material may be found in \citep{westervelt2003hybrid}.

\subsection{Bipedal Robotic Dynamics}
\label{sec:bipedsdynamics}
The bipedal robotic system we consider is the AMBER-3M robotic platform seen in Figure \ref{fig:amber_and_barrier}, modeled as an underactuated, planar five-link robot with point feet \citep{ambrose2017toward} whose physical parameters are reported in \cite[Table 1]{ma2014human}. The configuration coordinates $\mb{q}\in \mathcal{Q}\subset \mathbb{R}^5$ are given by
$

    \mb{q} = [q_{sf},\ q_{sk},\ q_{sh},\ q_{nsh},\ q_{nsk}]^\top,
$
with stance foot angle $q_{sf}$, stance knee angle $q_{sk}$, stance hip angle $q_{sh}$, non-stance hip angle $q_{nsh}$ and non-stance knee angle $q_{nsk}$. The continuous-time equations of motion, derived from the Euler-Lagrange equations, are given by:
\begin{align}
    \mb{D}(\mb{q})\ddot{\mb{q}}+\mb{C}(\mb{q},\dot{\mb{q}})\dot{\mb{q}}+\mb{G}(\mb{q})&=\mb{Bu}, \label{roboDyn}
\end{align}
where $\mb{D}(q)\in{\mathbb{S}}_{++}^{5}$ is the positive definite mass-inertia matrix,
$\mb{C}(q,\dot{q}) \in \mathbb{R}^{5\times 5}$ contains the centrifugal and Coriolis forces, $\mb{G}(q)\in \mathbb{R}^{5}$ contains the gravitational forces, $\mb{B}\in \mathbb{R}^{5\times 4}$ is the actuation matrix, and $\mb{u}\in\mathcal{U}\subset\mathbb{R}^4$ is the input. Note that this is the pinned model of the robot dynamics; for the unpinned model, refer to \citep{hereid2018dynamic}. For AMBER-3M, the number of inputs is one fewer than the degrees of freedom, meaning the system has one degree of underactuation.

Taking $p^v:\mathcal{Q}\to \mathbb{R}$ to represent the vertical position (height) of the swing foot, the admissible states are given by the \textit{domain} $\mathcal{D} = \{(q,\mb{\dot{q}})\in T\mathcal{Q}\ |\ p^v(\mb{q}) \ge 0\}$. The switching surface on which the impact events occur, also known as the \textit{guard}, is defined by:
\begin{align}
    \mathcal{S} = \{(\mb{q},\dot{\mb{q}})\in T\mathcal{Q} ~|~ p^v(q) = 0, \dot{p}^v(q,\dot{q}) < 0\}\subset \mathcal{D}.
\end{align}
With the full system state given by $\mb{x}=(\mb{q},\dot{\mb{q}})\in T\mathcal{Q}$, the impact dynamics \citep{westervelt2007feedback} are defined by a \textit{reset map} $\bs{\Delta}:\mathcal{S}\to \mathcal{D}$ relating pre-impact states $\mb{x}^-(t) \triangleq \lim_{\tau \nearrow t} \mb{x}(\tau)$ and post-impact states  $\mb{x}^+(t) = \lim_{\tau \searrow t} \mb{x}(\tau)$ via $\mb{x}^+(t) = \bs{\Delta}(\mb{x}^-(t))$. Combining these concepts and rearranging (\ref{roboDyn}) into control affine form yields the following \textit{hybrid control system}:
\begin{equation}\label{eq:hybrid_model_1}
\mathcal{HC} = \begin{cases} \dot{\mb{x}} = \mb{f}(\mb{x}) + \mb{g}(\mb{x}) \mb{u} \quad &\mb{x}^- \in \mathcal{D} \setminus \mathcal{S},\\
\mb{x}^+= \bs{\Delta}(\mb{x}^-) &\mb{x}^- \in \mathcal{S}.\end{cases}
\end{equation}
\subsection{Bipedal Robotic Control}
\label{sec:bipedcontrol}
Control of bipedal robotic systems centers around a \textit{phasing variable}, $\tau:\mathcal{Q}\to [0,1]$, given by:
\begin{equation} \label{eq:phase_variable}
    \tau(\mb{q}) = \frac{\delta_{hip}(\mb{q}) - \delta_{hip}^+}{\delta_{hip}^- -\delta_{hip}^+},
\end{equation}
where $\delta_{hip}:\mathbb{R}^5\to\mathbb{R}$ defined as $\delta_{hip}(\mb{q}) = [-l_t-l_f,\ -l_f,\ 0,\ 0,\ 0]\mb{q}$ is the linearized hip position with $l_t$ and $l_f$ the length of the robots tibia and femur, respectively. The constants \(\delta_{hip}^+\) and \(\delta_{hip}^-\) are the linearized hip positions at the beginning and the end of a step, ensuring that $\tau(\mb{q})$ increases monotonically in time within a step. Desired trajectories resulting in walking gaits for the robot can be rapidly synthesized via a hybrid zero dynamics framework (\cite{westervelt2003hybrid,ames2014rapidly}).
We are now well equipped to define the relative degree 2 (\cite{sastry1999nonlinear}) outputs $\mb{y}:\mathcal{Q}\to \mathbb{R}^4$ as the difference between the actual output $\mb{y}_a$ and the desired output trajectory $\mb{y}_d$:
\begin{equation}
    \mb{y}(\mb{q},\bs{\alpha}) \triangleq \mb{y}_a(\mb{q}) - \mb{y}_d(\tau(\mb{q}),\bs{\alpha}),
\end{equation}
with $\bs{\alpha}$ being the coefficients of a Bézier polynomial coming from the trajectory generation step. The actual output is given by the actuated coordinates:
$
    \mb{y}_a(\mb{q}) = \begin{bmatrix} \mb{0}_{4\times 1} & \mb{I}_{4\times 4} \end{bmatrix}\mb{q}.
$
The nominal controller for this system is then given by the proportional-derivative controller $\mb{k}_d(\mb{x}) = \mb{k}_{PD}(\mb{x}) \triangleq -\mb{K}_P \mb{y}(q) - \mb{K}_D \dot{\mb{y}}(q)$ with proportional gain $\mb{K}_P\in\mathbb{S}_{++}^4$ and derivative gain $\mb{K}_D\in\mathbb{S}_{++}^4$.

\begin{figure}[ht!]
\centering
\includegraphics[width=\textwidth]{Figures/confBarrier.pdf}
	\caption{\textbf{(Left):} Schematic diagram of the AMBER-3M robot with position coordinates. \textbf{(Center):} Schematic of the foot placement in the stepping-stone problem. The boundaries of virtual stepping-stones are captured via the blue and orange vertical lines. \textbf{(Right):} Virtual stepping stone width as function of the phase variable $\tau(\mb{q})$.}
	\label{fig:amber_and_barrier}
\end{figure}

\subsection{Control Barrier Functions for Stepping-Stones}

\noindent The stepping-stone problem is captured through the use of \textit{virtual} stepping-stones, which shrink over the course of a step to confine foot placement to a safe region defined on a targeted stone \citep{grandia2020multi}. The CBFs used to specify these foot position constraints are given by:
\begin{align}\label{eqn:ss-CBF1}
    h_1(q) &= R(\tau(q))-(O_x-F_x(q)),\\
    h_2(q) &= R(\tau(q))+(O_x-F_x(q)),\label{eqn:ss-CBF2}
\end{align}
where $F_x(q)$ is the horizontal position of the swing foot and $O_x>0$ is the horizontal position of the center of stepping-stone. The virtual stone width is given by the function $R:\mathbb{R}\to \mathbb{R}$:
\begin{equation}
    R(\tau(q)) = \frac{ar-1}{1  +ar(e^{-m(\tau(q)-1)}-1)}+1+r
\end{equation}
where  $m>0$ determines the decay rate of the barrier function, $(1+a)r$ is half of the targeted stone width, and $1+r$ defines the half the width of the virtual stepping-stone when $\tau = 0$. These functions are visualized in Figure \ref{fig:amber_and_barrier}.
The safety constraints can be interpreted as keeping the swing foot horizontal position in an interval centered at the middle of the stepping-stone, where the interval shrinks as $\tau$ increases.
As this formulation of CBFs is position-based and therefore relative degree two, we employ the exponential control barrier function (ECBF) extension technique \citep{nguyen2016exponential} to both CBFs to attain the relative degree 1 CBFs: $h_{e,i}(\mb{x}) \triangleq L_{\mb{f}}h_i(\mb{x})  + \alpha_{e}h_i(q).$

Combining the results in Section \ref{sec:learning} and \ref{sec:bipeds}, with $\hat{\mb{f}}$ and $\hat{\mb{g}}$ the nominal model in \eqref{eq:hybrid_model_1}, the final Stepping Stone QP \eqref{SS-QP} controller combines the robustifying term of \eqref{eqn:delta-CBF-QP} with the stepping-stone ECBF extensions of \eqref{eqn:ss-CBF1} and \eqref{eqn:ss-CBF2}:
\begin{align}
    \mb{k}(\mb{x}) =  \,\,\underset{\mb{u} \in \mathbb{R}^m}{\argmin}  &  \quad \frac{1}{2} \| \mb{u}-\mb{k}_{PD}(\mb{x}) \|_2^2 \tag{SS-QP} \label{SS-QP}\\

    \mathrm{s.t.} \quad & \quad L_{\hat{\mb{f}}}^2h_1(\mb{x})+L_{\hat{\mb{g}}}L_{\hat{\mb{f}}}h_1(\mb{x})\u + \alpha_e L_{\hat{\mb{f}}}h_1(\mb{x}) + \widehat{\delta}_1(\mb{x})
    \geq -\alpha(h_{e,1}(\mb{x})) \label{finalCBF1}\\
    & \quad L_{\hat{\mb{f}}}^2h_2(\mb{x})+L_{\hat{\mb{g}}}L_{\hat{\mb{f}}}h_2(\mb{x})\u + \alpha_e L_{\hat{\mb{f}}}h_2(\mb{x}) + \widehat{\delta}_2(\mb{x})
    \geq -\alpha(h_{e,2}(\mb{x})). \label{finalCBF2}
\end{align}

\noindent We assume that this \eqref{SS-QP} is feasible and we encountered no infeasibilities in simulation or experimentation.
\section{Simulation and Experimental Validation}
\label{sec:exp}
In this section we apply our episodic learning framework (Algorithm \ref{alg:episodic_learning}) to the AMBER-3M platform in both simulation with injected model uncertainty and on hardware with the model error inherent to real-world systems. In each instance the estimator $\widehat{\delta}$ was implemented as a neural network with two hidden layers of 50 hidden units using the ReLU activation function. The network was trained minimizing mean absolute error using mini-batch gradient descent. Mean absolute error was chosen over other loss functions for its robustness to outliers. The same controller (\ref{SS-QP}) was deployed in the RaiSim \citep{raisim} simulation environment and on the AMBER-3M hardware platform, as seen in the supplementary video (\cite{videoVideo}).
The complete learning code used in simulation and experiments can be found at \citep{git}.
\begin{figure}[b!]
\centering
\includegraphics[width=\textwidth]{Figures/SimHardHComb.pdf}
\caption{Simulation (S) and Hardware (H) data where model mismatch causes violations. \textbf{(Far-Left)}: Simulation where the barrier functions $h_1$ (solid blue) and $h_2$ (solid orange) are enforced via a \ref{eqn:CBF-QP}. The \ref{eqn:dumbreject} is also shown for $\overline{\delta}_1$ (dashed blue) and $\overline{\delta}_2$ (dashed orange), which results in more conservative behavior over many steps. \textbf{(Mid-Left)}: After three episodes of learning the \ref{SS-QP} in simulation, the maximum barrier violation decreases from 2.0 to 0.3 cm. \textbf{(Mid-Right)}: Hardware where the barrier functions $h_1$ (blue) and $h_2$ (orange) enforced via a \ref{eqn:CBF-QP}. \textbf{(Far-Right)}: After two episodes of learning on hardware, the maximum barrier violation decreases from 9.2 to 1.9 cm via the \ref{SS-QP}.}
\label{fig:eps}
\end{figure}

\subsection{Simulation}
The controllers and learning algorithm were first validated in simulation. Model error was introduced by increasing the inertia of all limbs on the true model by a factor of ten while maintaining constant mass. Due to the underactuated nature of the robot and the relationship between step length and zero dynamics stability, not every set of stepping stones is navigable, even if safety is perfectly enforced with respect to the CBFs. Therefore, a feasible stepping stone configuration was first generated for the robot to traverse with stones of 4 cm in width. Without knowledge of the modified model ($\widehat{\delta}_1(\mb{x}) = \widehat{\delta}_2(\mb{x}) = 0$), the controller did not satisfy the CBF constraints (\ref{finalCBF1}-\ref{finalCBF2}), resulting in a maximum violation at foot placement of 2.0 cm, causing the robot to miss the stepping stone and fall over. Three episodes of the PDL algorithm were run, after which the maximum violation was reduced to be 0.3 cm, only 15\% of the original violation. Additionally, the \ref{eqn:dumbreject} controller was implemented, which ensured safety but resulted in extremely conservative behavior, resulting in poor qualitative walking, i.e. harsh foot strikes and an over-bending torso. A comparison of the barrier functions $h_1$ and $h_2$ over the steps with these controllers can be seen in Figure \ref{fig:eps}.

\subsection{Hardware}
The same nominal model for the robot was used in the hardware experiments as in simulation, with model uncertainty presenting itself as significant friction in the joints, as well as imperfect mass and inertia measurements. The PDL algorithm was implemented on the AMBER-3M robot across a sequence of two episodes. The controllers ran on an off-board i7-6700HQ CPU @ 2.6GHz with 16 GB RAM, which computed desired torques and communicated them with the ELMO motor drivers on the 137 cm tall, 22 kg robot. The motor driver communication ran at 2kHz, and the \ref{SS-QP} ran at 1kHz.
The stepping stone configuration was specified to the controller with stones of 8 cm in width.
As with simulation, the CBF-QP resulted in a maximum violation of the barriers of 9.2 cm due to model error. After running the PDL algorithm for two episodes, the maximum violation of the barriers was 1.9 cm, only 21\% of the original violation, as depicted in Figure \ref{fig:eps}. Although learning improves our estimate of the safe set size, as discussed in Section \ref{sec:learning}, there is still uncertainty in the stone size. This was accommodated by utilizing physical stones which were 10 cm in width.
The 7.3 cm reduction in stone size mismatch captures the change from $|\delta|_\infty$ to $|\delta-\widehat{\delta}|_\infty$. Gait tiles for this improved traversal of the stepping stones are shown in Figure \ref{fig:gaitTiles}.

\begin{figure}[t!]
    \centering
    \includegraphics[width=\textwidth]{Figures/Hardware/GaitTiles.png}
    \caption{Gait tiles for Episode 2 of learning showing the AMBER-3M robot safely traversing a set of stepping stones. Notice the change in step width and added lean of the torso induced by the barrier functions.}
    \label{fig:gaitTiles}
    \vspace{-7mm}
\end{figure}

\section{Conclusion}
\label{sec:conclude}
In this paper we presented an episodic learning approach for reducing the impact of model uncertainty on safety-critical control using control barrier functions. Our method is able to learn a projected disturbance and incorporate the learned model information into an optimization-based controller, as demonstrated in both a high-fidelity simulation and on the AMBER-3M robot hardware platform. Future work includes theoretical analysis of the convergence properties of Algorithm \ref{alg:episodic_learning}, modification of the controller to ensure trajectories remain close to previously-seen data points to reduce the effects of generalization error, and an extension to three dimensional bipedal walking with vision sensing for stepping-stone identification.

\section{Hybrid motion planning}

\begin{figure}[!t]
	\begin{center}
	\vspace{3mm}
		\includegraphics[width=0.43\columnwidth]{Figures/Vision/V35_2.png}
		\includegraphics[width=0.46\columnwidth]{Figures/Vision/model.png}
		\caption{On the left is the rendering of the robot design, and on the right is the illustration of the configuration coordinates for the robot. The leg indices $\mathrm{l}_*$ are shown on the vertices of the \textit{body} link. Each leg has three actuated joints and equipped with a point foot.}
		\label{fig:config}
	\end{center}
	\vspace{-5mm}
\end{figure}

We consider the nonlinear model of quadrupedal locomotion as a hybrid dynamical system, which is an alternating sequence of continuous- and discrete-time dynamics. The order of the sequence is dictated by contact events. In comparison with bipedal walking, the increased number of contact points of quadrupeds increased the complexity of the hybrid model substantially. In this section, we introduce a unified model for quadrupedal behaviors including walking, ambling and trotting, based on which we design full model based optimal controllers as well as simulation validation.

\subsection{The robot}
The robot we study in this paper is the Vision 60-v3.2 (see \figref{fig:robot}), a quadrupedal robot built by Ghost Robotics, with the total weight of $26$ kg and maximum standing height of $0.5$ m. As shown in \figref{fig:config}, we model the quadruped as a $13$-link system: a \textit{body} link and four legs, each of which has three children links ---the \textit{hip}, \textit{upper} and \textit{lower} links. The configuration variables of Vision 60 are denoted by $q\in\mathcal{Q}\subset\mathbb{R}^n$ where $n=18$ is the total number of degrees of freedom (DOF) without considering any contact constraints. Utilizing the floating base convention \cite{Grizzle2014Models}, we can have $q^T =( q_b^T, q_l^T )$, in which $q_b\in\mathbb{R}^3\times\mathrm{SO}(3)$ represents the global Cartesian position and orientation of a frame attached to the \textit{body} linkage, and the local coordinates $q_l\in\mathbb{R}^{12}$ denote the $12$ joint variables: hip roll, hip pitch and knee angles. These angles are denoted by $\theta_{hr_j}, \theta_{hp_j}, \theta_{k_j}$ for the $j$-th leg, all of which are actuated by Brushless (BL) DC motors. This yields the system's total DOF to be $18$ and control inputs $u\in\mathbb{R}^{12}$.
With different scenarios of foot contacts with the ground, we have a mixture of overactuated, fully-actuated, and underactuated domains (i.e., phases) for the dynamics. Futher, we can define the state space $\mathcal{X}=T\mathcal{Q}\subseteq \mathbb{R}^{2n}$ with the state vector $x^T = (q^T,\dot q^T)$, where $T\mathcal{Q}$ is the tangent bundle of the configuration space $\mathcal{Q}$.

We now define the different modes of the quadrupedal dynamics, separated by discrete events. This fusion of continuous- and discrete-time dynamics yields a hybrid automaton model \cite{derIntro}. For more detailed definitions of hybrid modeling for bipedal robots, we refer the readers to \cite{Ames2013Humana}.

\subsection{The continuous-time domain: constrained dynamics}
Given the floating base coordinates, we can derive the unconstrained dynamics, i.e. without any contact constraints, by the Euler-Lagrange equations as:
\begin{align*}
    D(q)\, \ddot q + H(q,\dot q) &= B u
\end{align*}
where $D(q)\in\mathbb{R}^{n\times n}$ is the mass-inertia matrix, $ H(q,\dot q)\in\mathbb{R}^n$ includes the Coriolis, centrifugal and gravity terms and $B(q)\in\mathbb{R}^{n\times m}$ is the actuation matrix and $u\in\mathbb{R}^m$ is the torque applied at joints with $m=12$.

Now consider $k\in\{1,2,3,4\}$ feet standing on the ground, which means $\bar k=4-k$ feet are swinging in the air. This creates a variety of contact scenarios, hence generates different types of quadrupedal behaviors (see \figref{fig:graph}). For each contact situation, we associate a \textit{continuous domain}:

\vspace{-3mm}
\small
\begin{align}
	\label{eq:domain}
	\mathcal{D}_v := \{ (x,u) | h_v(q)=\dot h_v(q,\dot q)=0, N_v \succeq 0, h_{\mathrm{s},v}(q)\succeq 0\}
\end{align}
\normalsize
where $v\in\{1,2,3...\}$ is the domain index and $h_v(q)\in\mathbb{R}^k$ is the height of all standing feet with ground reaction force $N_v(x,u)\in\mathbb{R}^k$. The height of the other feet, referred as the swing feet, is shown by $h_{\mathrm{s}, v}(q)\in\mathbb{R}^{\bar k}$.
We then have the constrained dynamics for $\mathcal{D}_v$ as:
\begin{align}\begin{cases}
D(q)\, \ddot q + H(q,\dot q) = B_v\, u + J_{v}^T(q)\, \lambda_v \\
J_v(q)\, \ddot q  + \dot J_v(q,\dot q)\, \dot q =0
\end{cases}
\label{eq:eom-general}
\end{align}
where $J_v(q):={\partial p_v(q)} / {\partial q}$ represents the Jacobian matrix of the Cartesian position of the standing feet $p_v(q)\in\mathbb{R}^{3k}$, with the corresponding constraint wrench $\lambda_v\in\mathbb{R}^{3k}$.
Note that the actuation matrix $B_v$ is domain dependent. This is because the double and triple support phases (\figref{fig:graph}) create closed-chain structures that induce redundancy in control and constraints.
This is an underdetermined problem which often appears in multi-contact locomotion \cite{zhaomulti} that yields nonunique controllers for $u$. Therefore, we manually turn off the rear standing leg's hip pitch motor for double support phase and turn off the diagonal standing legs' hip pitch motors and the other standing leg's hip roll motors for triple supporting phase. This implementation in return yields underactauted dynamics for the full system. Now we can convert the constrained EOM \eqref{eq:eom-general} into a controlled ODE:
\begin{align}
\label{eq:ode}
	\dot x = f_v(x) + g_v(x) u.
\end{align}
To track a given set of time-based trajectories $\mathcal{B}_v(t)$, which will be detailed in a later section, we deployed an \textit{input output feedback linearization} controller:
\begin{align}
\label{eq:io}
	u_{\mathrm{io}}(x,t) = \mathcal{A}(x) \big(  \mathcal{L}(x,t) - 2\varepsilon y(q,t) - \varepsilon^2 \dot y(q,\dot q, t)  \big).
\end{align}
with the outputs $y(q,t) = y_a(q)-\mathcal{B}_v(t)$ and $\varepsilon>0$. In this formulation, we chose the \textit{actual outputs} $y_a(q)$ as all of the actuated joints. The notations follow directly from Eq. (28) of \cite{Ames2013Humana}. The time-based controller has been justified for improved robustness over state-based methods both in theory \cite{DURUS_timedependent} and experiments \cite{ma2017bipedal}. As a result, \eqref{eq:io} forces the system to converge to a desired gait exponentially, that is, $y_a(q)\rightarrow \mathcal{B}_v(t)$.
The output dynamics become
\begin{align}
    \label{eq:ydyn}
	\ddot y = -2\varepsilon \dot y -\varepsilon^2 y
\end{align}
for which the origin $(y,\dot{y})=(0,0)$ is exponentially stable.

\begin{figure}[!t]
\vspace{2mm}
	\begin{center}
		\includegraphics[width=0.9\columnwidth]{Figures/hybrid/directedGraph.png}
		\caption{The cyclic directed graph for the multi-domain hybrid dynamics of walking, ambling and trotting gaits.
		The dashed lines represent a \textit{relabeling} map \cite{Ames2013Humana} that flips the left and right legs' contact attributes.

		}
		\label{fig:graph}
	\end{center}
	\vspace{-5mm}
\end{figure}

Prior to introducing the edges and the discrete dynamics, some assumptions are necessary to construct a feasible model both for computation and experiment:
\begin{itemize}
	\item There is no ground slippage. This is partially guaranteed by enforcing a friction cone condition. However, slipping locomotion has been observed on quadrupedal animals for energy efficiency \cite{gamus2013analysis, Maslipery}.

	\item The ethology studies \cite{Muybridgeanimal} have observed a pattern of $4\times2$-domain\footnote{The term $m \times 2$ denotes a gait with $m$-domain. See \figref{fig:graph} for the directed graph of the gaits with $m$ domains. The second half of the motion is directly a left-right mirror (i.e., symmetry) of the gait. See \figref{fig:gaits} for the full motion with $m \times 2$ domains.}
	locomotion on quadrupedal animals, and the authors have investigated this modeling method in \cite{Hamed_Ma_Ames_Vision60}. But in this paper, we assume the \textit{stance leg transition domains ---one leg strikes while another leg lifts} are instantaneous and passive for walking and ambling. Hence $\mathbf{s}_1$ and $\mathbf{s}_2$ in \figref{fig:graph} become edges and we can have a $2\times 2$-domain behavior for \textit{walk} and \textit{amble}.
	But a $4\times 2$-domain model is still used for \textit{trot}.

\end{itemize}

\subsection{The discrete-time domain: impact and lift-off}
On the edge of $\mathcal{D}_v$ in \eqref{eq:domain}, one of the conditions reaches its bound. Thus we have two switching mechanisms:
\begin{itemize}

	\item \textit{Lift off}: a standing foot of leg ${l_*}$ lifts off from the ground, meaning
	$N_v^{l_*}(x,u) = 0$.

	\item \textit{Impact}: a swing foot of leg $l_*$ impacts the ground, meaning $h_{s,v}^{l_*}(q) = 0$, $\dot h_{s,v}^{l_*}(q,\dot q) < 0$.
\end{itemize}
For \textit{lift off}, an identity map $x^+_{v+1} = x^-_{v}$ is sufficient to represent the transition from current to the next domain, where $x^-_{v}$ is the state at the end of the domain $\mathcal{D}_v$ and  $x^+_{v+1}$ is the state at the beginning of $\mathcal{D}_{v+1}$. However, the \textit{impact} shall cause a jump in the velocity terms. To describe this, we make use of a plastic impact model \cite{Grizzle2014Models},
\begin{align}
\label{eq:impact}
	\begin{bmatrix} D & -J_{v+1}^T\\ J_{v+1} &0 \end{bmatrix}
	\begin{bmatrix} \dot q^+_{v+1}\\ \Lambda \end{bmatrix}
	=
	\begin{bmatrix} D \dot q^-_v\\ 0  \end{bmatrix}
\end{align}
by using the conservation of momentum while satisfying the next domain's holonomic constraints.

\input{sections/opt}

\subsection{Gait Synthesis}
\label{sec:gaitsyn}
The desired trajectory $\mb{y}_d(\tau(\mb{q}),\alpha)$ is given by a degree $N$ Bézier polynomial
\[\mb{y}_d^i(\tau) = \sum_{k=0}^{N}\alpha_k^i \frac{N!}{k!(M-k)!}\tau^k (1-\tau)^{M-k}, \quad i\in \{1,...,m\} \]
With the relative degree 2 output defined for the hybrid control system $\Sigma$, the zero dynamic surface is defined as:
\begin{equation*} \label{eq:zdsurface}
    \mathcal{Z}(\delta_{hip}^-,\delta_{hip}^+,\bs{\alpha}) = \{ \mb{x} \in \mathcal{X}  \big| \mb{y}(\mb{q}) = 0 \mbox{ and } L_\mb{f} \mb{y}(\mb{q},\dot{\mb{q}}) = 0 \}
\end{equation*}

For underactuated systems, \cite[Thm. 2]{ames2014rapidly} proved that with sufficiently rapidly exponentially stabilizing outputs, an exponentially stable periodic orbit in zero dynamics surface is an exponentially stable periodic orbit for the full-order dynamics.  The goal of the gait optimization is to find a periodic gait in the zero dynamics surface.

The gait optimization is then performed with a feedback linearizing control law
\begin{equation}
    \u(\mb{x}) = -\left(L_\mb{g} L_\mb{f} \mb{y}(\mb{x})\right)^{-1} \left(\frac{1}{\epsilon}K_D L_\mb{f} \mb{y}(\mb{x}) +\frac{1}{\epsilon^2}K_P \mb{y}(\mb{x}) + L^2_\mb{f}\mb{y}(\mb{x})\right)
\end{equation}
which results in rapidly exponentially stabilizing output dynamics as $\Ddot{\mb{y}} = -\frac{1}{\epsilon}K_D \dot{\mb{y}} -\frac{1}{\epsilon^2}K_P \mb{y}$.

The optimization is performed to find a periodic orbit on the two-dimensional zero dynamics submanifold of $\mathcal{X}$, i.e. to find the set of parameters $\delta_{hip}^-,\delta_{hip}^+,\bs{\alpha}$ for desired outputs such that
\[ \Delta(\mb{x}^-) \in \mathcal{Z}(\delta_{hip}^-,\delta_{hip}^+,\bs{\alpha}) ,\quad  \forall \mb{x}^- \in \mathcal{S} \cap \mathcal{Z}(\delta_{hip}^-,\delta_{hip}^+,\bs{\alpha}) \tag{HZD}\]
under a list of physical constraints including foot clearance, torso moving range, friction cone constraints and actuation limits.

The gait optimization is a nonlinear program (NLP) performed under a direct collocation framework \cite{hereid2017frost}:
\begin{align*}
    \mathcal{X}^*  = \operatorname*{argmin}_{\chi} \quad & \Phi(\mathcal{X})\\
    \textrm{s.t.   }& \text{HZD}\\
    &\text{physical constraints}\\
\end{align*}
where $\mathcal{X}$ is a set of parameters in the NLP including $\delta_{hip}^-,\delta_{hip}^+$, and $\bs{\alpha}$.