\section{Introduction}

Robots that act in the real world often require the use of contacts.
Contacts are important to locomote by walking \cite{Escande2013, Orthey2013, Bouyarmane2018} or climbing \cite{B06}, to manipulate the environment \cite{Toussaint2020} or to grasp objects \cite{CHR19}. Algorithms that enable robots to move using contacts are studied in the field of contact-based motion planning \cite{Escande2013, TPM18, Kingston2019}.
The objective of contact-based motion planning is to complete a manipulation or locomotion task by reaching a goal configuration while incorporating contact constraints where chosen points on the robot have to be in contact with chosen surfaces in the environment.

One approach to deal with contact constraints is projection-based constraint planning \cite{Berenson2011, Kingston2019}. In projection-based constraint planning, we implicitly define contact constraints, which we can sample using dedicated projection methods \cite{Berenson2011}. Such a framework has been used to plan contact motions for the Robonaut 2 robot \cite{Kingston2019}. While projection-based planning works well for fixed contact planning \cite{Kingston2019}, there does not yet exist an extension to incorporate sliding constraints or transitions between sliding contacts. Sliding contacts however, are a fundamental requirement if we want to explore the shape of an object \cite{DET17}, adjust an object grasp \cite{Shi2017}, keep the robot in balance \cite{Samadi2020}, minimize state uncertainty \cite{Pall2018}, plan more efficiently \cite{Roussel2019} or animate computer characters \cite{Hu2019}.

To overcome those limitations of constraint planning, we develop an extension of constraint planning for sliding contacts. This extension consists of the following items:

\begin{enumerate}
    \item Definition of sliding contact constraint and transition constraint by extending the constraint-based planning framework \cite{Kingston2019}.
    \item Sampling method for the transition constraint, which uses constraint modes for efficient sampling.
    \item Evaluations with various sampling-based planning algorithms for one-step contact cycles using manipulator arms with two, three and seven internal degrees of freedom (dof).
\end{enumerate}

Note that we concentrate in this paper exclusively on one-step sliding contact cycles, because solving such cycles robustly is of fundamental importance for planning long-horizon contact sequences \cite{Toussaint2020}. Furthermore, while our approach makes simplifying assumptions on the control of sliding contacts, our motivation is to eventually have a concise framework in which we can guarantee completeness and optimality. As an example, Figure \ref{fig:example} shows a one-step sliding contact and transition cycle for a robot with three internal dof and two contact points, where we exhaustively sample the transition constraints of sliding contacts to find a solution path.

\begin{figure}[t]
\centering
	\includegraphics[width=\linewidth]{graphics/pullFigure.png}
	\caption{We present a method to plan constraints for sliding contact and transition motions on arbitrary mesh surfaces, where predefined contact points of a robot (yellow) are required to stay in constant contact with a surface or to transition.}
	\label{fig:example}
\end{figure}
\section{Related Work}

Planning contact sequences through an environment can be addressed using constraint-based motion planning \cite{Berenson2011, Kingston2019}. This approach is based on building probabilistic road maps or trees by randomly sampling configurations, which are projected onto the contact constraints to make them feasible \cite{Berenson2011}. The transition between contact constraints is often captured using a constraint graph \cite{Mirabel2017}, through which we can identify sequences of contacts to reach a goal \cite{Escande2013}. Such sequences are often used to plan more complicated contact or manipulation tasks \cite{Garrett2020, Toussaint2020}. Because contact constraints can often not be described analytically, we have to either learn them \cite{Fernandez2020} or use methods to project
onto them \cite{Berenson2011}. To plan longer contact sequences, it often becomes necessary to exploit the structure of the state space. Different approaches exist, for example by factoring the state space for more efficient sampling \cite{Vega2020}, by computing guide paths from a simplified robot geometry \cite{TPM18} or by precomputing convex regions in workspace which we can exploit using mixed-integer programming \cite{Deits2014}. Our approach is complementary, in that we add sliding transition constraints to the constraint-based framework \cite{Kingston2019}, which makes it possible to plan sliding transition constraints using sampling-based planner \cite{Lav06} or asymptotically optimal planner \cite{Karaman2011}.

While constraint-based planning is often used to plan non-sliding contacts, there exists substantial research on the control aspect of sliding motions \cite{Trinkle1990, Cruciani2020}. To execute and control sliding motions, we could use physical models which make use of friction cone models \cite{Samadi2020} to directly optimize sliding motions to push objects along planar surfaces \cite{Toussaint2020}. An alternative approach would be to learn the parameters of a controller for specific object categories \cite{Li2019}. Recent results using learned sliding controllers can be used for dexterous in-hand manipulation, where objects are either rotated in-hand without making contact with the environment \cite{Kumar2016, Andrychowicz2020} or features in the environment are used as helper tools \cite{Cruciani2020}. Our approach is complementary, in that we make simplifying assumptions on the control of sliding contacts (deterministic, no uncertainty) and concentrate on the computational challenges of planning sliding contacts \cite{Lee2015}, which could be used to provide guarantees on completeness and optimality.
\section{Background}
This section provides some background knowledge.
Large parts of this section are based on Part II of \cite{Lav06} and \cite{KMK18}.

A \textbf{configuration} $q$ defines the independent variables of a given robot needed to uniquely specify the robot's position relative to a reference frame.
Considering a manipulator arm with $m$ number of rigid links connected by $n$ joints that is fixed at one end, one possible configuration would define $n$ internal variables plus external variables to describe a root frame fixed to the robot.
The \textbf{configuration space} $Q$ of a robot is the set of all such possible configurations $q$.
To avoid collisions with any obstacles or the robot itself, typically a free space $\Q_{\text{free}} \subseteq Q$ is defined that only contains collision-free configurations.

A basic motion planning problem is defined as finding a continuous path in $\Q_{\text{free}}$ that connects from a given start configuration $\q_{\text{start}}  \in \Q_{\text{free}}$ to a goal configuration $\q_{\text{goal}} \in \Q_{\text{free}}$.
Explicitly computing $\Q_{\text{free}}$ is a complex problem and the complexity grows with an increase of a robot's degrees of freedom \cite{Canny1988}.
We use the concept of sampling-based motion planning to avoid this by working with an implicit representation of $\Q_{\text{free}}$ probed through sampling strategies.

\textbf{Sampling-based motion planning} \cite{Lav06} relies on sampling collision-free configurations $q$ in $\Q_{\text{free}}$ which are connected to a tree or a graph.
With a longer planning time and increasing numbers of samples, more space can be mapped out.
If the problem is solvable, the probability of finding a path then converges to one if time goes to infinity.
The algorithms provide an efficient solution for finding feasible paths for high-dimensional problems.

\textbf{Constraint-based planning} can be used whenever a robot performs tasks that limit its possible motions to subspaces of the configuration space \cite{KMK18}.
When we use a contact-based method with point contacts, sliding contacts or contact transitions, we put such limits on a robot's configuration space $Q$.
In order to express these limits, we formulate a contact as a single task constraint and a transition as a set of successive task constraints.
We extend the motion planning objective of finding collision-free configurations to simultaneously satisfy given constraints.
These constraints are specified by a constraint function $$f(q): Q \rightarrow \mathbb{R}^{n}$$ that is satisfied when the real-value vector  $f(q) = \textbf{0}$ for a given configuration $q$.
This function is then used to construct an implicitly constrained configuration space. $$X = \{q \in Q \mid f(q) = \textbf{0}\}$$
It contains all configurations that satisfy the defined constraint.
We can now define a configuration space that includes all collision-free configurations from $\Q_{\text{free}}$ that also satisfy the constraints.
$$\X_{\text{free}} = X \cap \Q_{\text{free}}$$
A basic constrained motion planning problem is thus defined as finding a continuous path in $\X_{\text{free}}$ that connects from a given start configuration $\q_{\text{start}}  \in  \X_{\text{free}}$ to a goal configuration $\q_{\text{goal}} \in  \X_{\text{free}}$.

We assume that a robot has $k$ \textbf{contact points}, which are predefined points on the robots geometry.
To each contact point, we associate a mapping $f_k$ which maps a joint configuration and a surface to the distance of the $k$-th contact point to the nearest point on the surface.
For a surface $s$ in a workspace $W$, we say that the contact point fulfills a contact constraint if $f_k(q,s)=0$.

A continuous motion for a contact point which fulfills the contact constraint over a surface is a \textbf{sliding contact}.
A \textbf{contact transition} is the process of breaking a contact on an initial surface $s_I$, moving freely and creating a contact on a goal surface $s_G$.
In configuration space, we thereby transition from the constraint surface of all configurations $q$ in $Q$ for which $f_k(q,s_I)=0$ to being constraint-free to being in the set of configurations for which $f_k(q,s_G)=0$.
To keep track of sliding constraints and their transitions, we make use of the concept of a constraint graph \cite{Mirabel2017, Kingston2019}, which is a graph with vertices being active constraints and edges being possible transitions between active constraints. In the next section we elaborate on such a constraint graph.
\section{Sliding Contacts and Transition Constraints}

The goal of this work is realizing a new method for planning a one-step contact cycle for sliding contacts and the transition between two given surfaces. For this purpose, we develop a sliding contact constraint, a transition constraint and a sampler to efficiently exploit the transition constraint.

The sliding contact constraint keeps a contact point in constant touch with the surface of an object, while the transition constraint encompasses a contact break from a given start surface, the contact-free motion toward a goal surface and a contact creation at the goal.
The steps of a transition motion are divided into three separate transition modes.
We combine the concept of our contact constraint with a mode sampler that samples these three modes uniformly.

For a robot with $k$ contact points we predefine each to be either a constant sliding contact point or one that performs a transition and we specify $k$ constraints accordingly.
Only these $k$ designated contact points can be in contact, the rest of the body is planned to avoid any collision.

\subsection{Explanatory Example}

In Figure~\ref{fig:constraintGraph} we show a constraint graph example for a robot with $k = 2$ contact points and two contact surfaces. This graph corresponds to the situation of the robot in Figure~\ref{fig:2dof}, which is initially in contact with surface one (bottom box) and has to transition one contact point to surface two (top box) while the other contact point can slide along surface one.
The constraint graph depicts the three different modes a robot's contact point can be in and which modes can be accessed afterwards.
The $k$-th entry of a tuple corresponds to the assignment of the $k$-th contact point to the according surface or to be free. In our example, each contact point can thus be in one of three modes.

\begin{itemize}
	\item \textbf{Mode 0} contact point is unconstrained.
	\item \textbf{Mode 1} contact point is in contact with surface one.
	\item \textbf{Mode 2} contact point is in contact with surface two.
\end{itemize}

State (0,0) describes a free floating robot without any contacts.
From there either contact point can make contact with object surface one or two.
This mode change is depicted by the edges between the different states.
State (1,1) corresponds to a robot with both contacts on the same initial surface $1$ (see the green stance in Figure \ref{fig:2dof}) and state (2,2) corresponds to a double contact with surface $2$.
A step cycle describes the process of a joint breaking contact with one surface and creating a new contact on another surface.
At least one of the two joints is in contact with one of the surfaces during the whole process.
A step cycle is completed by changing modes along the depicted edges starting for example at state (1,1) and ending in state (2,1).
From there, the cycle can be repeated to form a step sequence of contacts and transitions.

Highlighted in green are the two particular one-step cycle examples. Starting from state (1,1), we consider the task of moving along the edges over (1,0) or (0,1) to state (1,2) or state (2,1), respectively.
This corresponds to breaking the first respectively second contact point from the first surface, freely moving it and finally making contact with the second surface while the second respectively first contact point stays in sliding contact with the first surface.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{graphics/constraintGraph.png}
    \caption{Constraint Graph. This graph depicts the possible modes and mode changes of a robot with $k = 2$ contact joints and two contact surfaces. Relevant for this work are the parts highlighted in green.}
    \label{fig:constraintGraph}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=0.48\linewidth]{graphics/2dofBottomUp.png}
	\includegraphics[width=0.48\linewidth]{graphics/2dofBottomUpPath.png}
	\caption[Two dof manipulator rectangle scenario]{Left: Two dof rectangle scenario. $q_{start}$ is depicted in green and $q_{goal}$ in red. Right: Path with intermediate step.}
	\label{fig:2dof}
\end{figure}

\subsection{Sliding Constraint\label{chap:contactconstraint}}

First we present the sliding contact constraint that is put on the $k$-th contact point of the robot specified to be in sliding contact.
We compute the distance between this contact point and the closest point on a given surface $S$ at the current configuration $q$ and assign this value to the constraint function $f_k(q,S)$.

Algorithm \ref{alg:contact} shows the pseudo-code for computing sliding constraints.
The inputs are the robot $r$, the contact point on the robot $p$ and the given contact surface $S$.
First the contact point's position in world coordinates $p_w$ is determined.
The distance $d$ is first declared infinite.
We iterate over partitions of the whole $S$ (meaning the triangles of a given surface mesh).
During each iteration we save the coordinate on one such partition $s$ that is closest to $p_w$ into $c$ and the distance between these two coordinates into the variable $d_s$.
We then compare the two distance variables and if $d_s$ is smaller than $d$, we assign the value of $d_s$ to $d$.
After this loop, we have determined the coordinate of the surface partition closest to the current position of the robot contact joint and return this smallest distance $d$.

\begin{algorithm}
 \KwInput{Robot $r$, Contact point $p$, Surface $S$}
 \KwOutput{Distance d}
 $p_w \leftarrow getWorldPosition(r, p)$\;
 $d = +\infty$\;
 \For{$s\text{ in }S$}{
  $c \leftarrow closestPosition(s, p_w)$\;
  $d_s \leftarrow getDistance(c, p_w)$\;
  \uIf{$d_s < d$}{
    $d \leftarrow d_s$\;
   }
  \textbf{return} $d$\;
 }
 \caption{Sliding Constraint}
 \label{alg:contact}
\end{algorithm}

\subsection{Transition Constraint}
This section presents the implementation of our transition constraint.
It combines the sliding contact constraint with a constraint-free state into three transition modes.

Algorithm \ref{alg:transition} shows the pseudo-code for the transition constraints.
It requires input on the robot $r$, the contact point on the robot $p$, the transition mode $mode$, the desired contact start surface $S_I$ and the goal surface $S_G$.
First the contact point's position in world coordinates $p_w$ is determined.
The distance $d$ is first declared infinite.
The transition $mode$ of the current iteration is set by our sampling method explained in the next section.
We enforce a different constraint depending on the mode.
In the case of $mode = 0$ the returned distance $d$ is 0.
This corresponds to a contact-free state and by definition the constraint is satisfied.
If $mode = 1$ then we compute a sliding contact constraint using surface $S_I$.
The last case of $mode = 2$ computes a sliding constraint on surface $S_G$.

\begin{algorithm}
 \KwInput{Robot $r$, Contact point $p$, Mode $mode$, Start surface $S_I$, Goal surface $S_G$}
 \KwOutput{distance d}
  $p_w \leftarrow getWorldPosition(r, p)$\;
  \uIf{$mode = 0$}{
    \textbf{return} $d = 0$
   }\uElseIf{$mode = 1$}{
    $d \leftarrow slidingConstraint(r, p, S_I)$\;
    \textbf{return} $d$\;
  }\uElseIf{$mode = 2$}{
    $d \leftarrow slidingConstraint(r, p, S_G)$\;
    \textbf{return} $d$\;
 }
 \caption{Transition Constraint}
 \label{alg:transition}
\end{algorithm}

\subsection{Sampling Method}

In order to use the transition constraint for path planning, we implement a method to sample the three different constraint modes.
The mode is set before calling the transition constraint algorithm and the transition constraint is enforced accordingly.

In Algorithm \ref{alg:sampler} we show the pseudo-code for sampling transition modes.
The input is the configuration space $Q$ and a list of the $k$ specified constraints $C$. Each constraint could be either a sliding or a transition constraint. We first sample a configuration $q$ uniformly in $Q$.
We then iterate over the constraints and check each entry if it is a transition constraint.
For each transition constraint we then assign a uniformly sampled integer from $[0, 1, 2]$ to the transition $mode$. Once the modes are set, we call the projection method which minimizes the constraint distance. Details on the project method can be found in \cite{Kingston2019}. To verify that we successfully projected, we then compute the distance and return the projected configuration.

\begin{algorithm}
 \KwInput{Configuration Space Q, Constraints C}
 \KwOutput{Configuration q, Distance d}
 \SetAlgoLined
 $q \leftarrow sampleUniform(Q)$\;
 \For{$c\text{ in }C$}{
  \uIf{$c$  \textbf{is}  $TransitionConstraint$}{
   $c.mode \leftarrow uniformSampledInt(0, 2)$\;
   }
 }
 $q \leftarrow C.project(q)$\;
 $d \leftarrow C.distance(q)$\;
 \textbf{return} $q, d$\;
 \caption{Transition Mode Sampler}
 \label{alg:sampler}
\end{algorithm}

\section{Evaluation and Results}

We test our transition constraints with manipulator arms of two, three and seven internal dof and different obstacle scenarios in two-dimensional (2D) and three-dimensional (3D) space.
The scenarios are run ten times with each planning algorithm and with a specified computation time.
We compare five sampling-based planning algorithms, which performed well in initial experiments. In detail, we use Rapidly Exploring Random Tree (RRT) \cite{Kuffner2000}, Sparse Stable RRT (SST) \cite{Li2016}, Search Tree with Resolution Independent Density Estimation (STRIDE) \cite{Gipson2013}, Probabilistic Roadmap Method (PRM) \cite{kavraki_1996} and Sparse Roadmap Spanners (SPARS) \cite{dobson_2014}. We evaluate the performance of our method by comparing the average computation times of the different algorithms.

\subsection{Programming Setup}

The sliding contact and transition constraints and sampler are implemented as extension of the constraint-based planning framework \cite{Kingston2019} of the Open Motion Planning Library \cite{OMPL, OMPLBenchmark}. To simulate and visualize the output, we use Kris' Locomotion and Manipulation Planning Toolbox \cite{Klampt}. The code is written in the programming language C++ and is freely available\footnote{\url{https://github.com/mtkhoury/MotionPlanningExplorerGUI/tree/contact_feature}}.

\subsection{Two dof Manipulator Arm}

Our first scenario is a 2D scenario with two rectangular objects and a manipulator arm with two rigid links, two joints and two contact points (yellow) to make contacts with (Figure \ref{fig:2dof}).
The robot starts in full contact with the lower rectangle and is required to slide along the surface with the first contact point (left), while making a sliding transition with the second contact point (right). A computed path is depicted in purple. In its goal position the robot keeps one contact with the initial surface and makes one contact with the upper rectangle.

Figure \ref{fig:time2dof} shows the average computation time of ten runs of this scenario with a specified maximum planning time of 0.5 seconds.
We can see that all planning algorithms successfully find a path in the allocated time.

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=\linewidth]{graphics/benchmark_02D_bottomUp_firstLinkFixed_0_5seconds.jpg}
	\caption[Computation Time Graph of Figure \ref{fig:2dof}]{Computation Time. Two dof Scenario. Maximum planning time of 0.5 seconds.}
	\label{fig:time2dof}
\end{figure}

\subsection{Seven dof Manipulator Arm}

The next scenario contains three rectangular objects and a manipulator arm with seven joints and two contact points (Figure \ref{fig:7dofCube}).
The robot starts with one contact point on the lower rectangle and the other contact point on the left-hand rectangle.
The task is to transition the upper end joint toward the right-hand rectangle.
To its goal position the robot has to slide the lower contact point along the lower rectangle's surface and has to make contact with the upper contact point on the upper right rectangle.

Figure \ref{fig:time7dofCube} shows the ten run average computation time graph of this scenario with a specified maximum planning time of 5 seconds.
We can observe that all planning algorithms find a path successfully in five seconds or less (on average). Both SPARS and SST require more time in comparison to the other three planners but they still solve the problem quickly.

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=0.48\linewidth]{graphics/7dofPlatforms.png}
	\includegraphics[width=0.48\linewidth]{graphics/7dofPlatformsPath.png}

	\caption[Seven dof manipulator rectangle scenario]{Seven dof manipulator rectangle scenario. $q_{start}$ is depicted in green and $q_{goal}$ in red.}
	\label{fig:7dofCube}
\end{figure}

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=\linewidth]{graphics/benchmark_07D_platforms_firstLinkFixed_5seconds.png}
	\caption[Computation Time Graph of Figure \ref{fig:7dofCube}]{Computation Time. Seven dof Scenario. Maximum planning time of 5 seconds.}
	\label{fig:time7dofCube}
\end{figure}

\subsection{Three dof Manipulator Arm on Sphere}

The last scenario contains a sphere, a cuboid object and a manipulator arm with three dof (Figure \ref{fig:3dofSphere}).
The robot starts in full contact with the sphere and has to transition up along the sphere to make contact with the upper right cuboid. During the transition, we restrict one contact point to remain in sliding contact with the sphere. A solution path is shown on the right.

Figure \ref{fig:time3dofSphere} shows the average computation time of ten runs of this scenario with a specified maximum planning time of 30 seconds.
We can see that not all planning algorithms find a path in the maximum planning time.
The algorithms STRIDE and PRM solve the planning problem the fastest in under 5 seconds and SPARS exceeds the time limit.

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=0.48\linewidth]{graphics/3dofSphere.png}
	\includegraphics[width=0.48\linewidth]{graphics/3dofSpherePath.png}
	\caption[Three dof manipulator sphere scenario]{Left: Three dof cuboid scenario. $q_{start}$ is depicted in green and $q_{goal}$ in red. Right: Path with intermediate state.}
	\label{fig:3dofSphere}
\end{figure}

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=\linewidth]{graphics/benchmark_03D_sphere_firstLinkFixed_30seconds.jpg}
	\caption[Computation Time Graph of Figure \ref{fig:3dofSphere}]{Computation Time. Three dof sphere scenario. Maximum planning time of 30 seconds.}
	\label{fig:time3dofSphere}
\end{figure}

\subsection{Limitations}

During the setup of our testing scenarios we came across a limitation due to the different ways to set the projection.
One setup, shown in Figure \ref{fig:firstFixed}, uses the first contact point (end of thinner link) to stay in contact with the initial surface while the last contact point (end of thicker link) transitions.
The other setup is vice-versa, the first contact point transitions and the last contact point keeps contact (Figure \ref{fig:lastFixed}).
Figure \ref{fig:timeFirstFixed} and Figure \ref{fig:timeLastFixed} show the computation time graphs of the aforementioned setups with a maximum planning time of 1 second (for better visibility) and 30 seconds respectively.
We can see a significant difference in the planning duration until a path was found.
The setup where the first contact point is in sliding contact was solved in under two seconds by all five algorithms.
In the other setup the path was found after a much longer time by three of the algorithms and the other two exceeded the maximum planning time.

This observed limitation may stem from the way the configurations are projected. In this work, we only consider setups in which the root of the robot is positioned on the robot's first contact joint.
When we place the root on the robot's end joint that underlies a transition constraint, the configurations are projected differently, which affects the sampling of the transition constraint. Further research is needed to determine why a discrepancy between projections exist and which projection would be most useful for a specific robot and task.

\begin{figure}[hpbt]
	\centering
	\begin{minipage}{0.485\linewidth}
		\centering
		\includegraphics[width=\textwidth]{graphics/3dofBottomUpFirstFixed.png}
		\caption[Three dof manipulator First Joint Sliding]{Three dof with first joint in sliding contact. $\q_{\text{start}}$ is depicted in green and $\q_{\text{goal}}$ in red.}
		\label{fig:firstFixed}
	\end{minipage}\hfill
	\begin{minipage}{0.485\linewidth}
		\centering
		\includegraphics[width=\textwidth]{graphics/3dofBottomUpLastFixed.png}
		\caption[Three dof manipulator Last Joint Sliding]{Three dof with last joint in sliding contact. $\q_{\text{start}}$ is depicted in green and $\q_{\text{goal}}$ in red.}
		\label{fig:lastFixed}
	\end{minipage}
\end{figure}

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=\linewidth]{graphics/benchmark_03D_bottomUp_firstLinkFixed_1second.jpg}
	\caption[Computation Time Graph Figure of \ref{fig:firstFixed}]{Computation times for the first joint in sliding contact. Maximum planning time of 1 second.}
	\label{fig:timeFirstFixed}
\end{figure}

\begin{figure}[hpbt]
	\centering
	\includegraphics[width=\linewidth]{graphics/benchmark_03D_bottomUp_lastLinkFixed_30seconds.png}
	\caption[Computation Time Graph Figure of \ref{fig:lastFixed}]{Computation times for the last joint in sliding contact. Maximum planning time of 30 seconds.}
	\label{fig:timeLastFixed}
\end{figure}

Finally, we can also see that the required computation time differs depending on the scene setup.
We considered three kinds of objects in our scenarios: rectangles in the 2D scenarios, and cuboids and a sphere in the 3D scenarios.
The sphere in the 3D scene leads to a longer computation time than the cuboid and rectangular objects.
This can be attributed to a much higher number of surface polygons that has to be iterated over during the planning process. Note that analytical expressions exist for the sphere and using mesh polygons is not strictly necessary. However, our goal is to have a method which works for any mesh-based object representation. It is clear that a 3D scenario leads to more potential contact surfaces and thus also results in more polygons. To address this problem in future work, we believe it to be advantageous to use additional projection methods \cite{Orthey2020IJRR}, which could create equivalence classes of polygons.
\section{Conclusion}

In this work we contributed a geometric approach to planning sliding contacts using constrained planning and developed a sampler method for sliding transition constraints.
We formulated two constraints for a robot's contact joints: a sliding contact constraint and a transition constraint.
For the transition constraint we factored the problem into three constraint modes and implemented a sampling method that samples them uniformly.
We addressed the planning of a single contact joint breaking contact, moving it toward a goal and creating a new contact while the other joints remain in fixed or sliding contact.

Our tests and evaluations using sampling-based planning algorithms showed that the presented concept works on manipulator arms in 2D and 3D scenarios.
We observed that the planning time significantly depends on the projection used for the transition motion. To overcome this limitation, we believe it fruitful to further study the structure of the projection itself \cite{Berenson2011, Orthey2020IJRR} or to learn explicit representations of transition constraints \cite{Fernandez2020}.

Despite these limitations, we can successfully combine sliding contacts with transition constraints in a constraint-based sampling framework.
Using this approach we were able to solve planning problems for a robot with sliding contact points, which is fundamental and a first step towards probabilistically complete and optimal planning of long-horizon contact sequences which involve sliding contacts.


\section{Introduction}
Entanglement, a very powerful and efficient quantum resource, is a cornerstone of many quantum communication and quantum computation protocols.
A lot of quantum schemes depend heavily on the property of entanglement \cite{horodecki2009entanglement}, such as quantum key distribution \cite{ekert1991qkd}, quantum teleportation \cite{bennett1993teleporting}, quantum metrology \cite{giovannetti2006metrology}, and quantum sensing \cite{degen2017quantumsensing}.

There is no doubt that the preparation of entangled quantum states is of great concern in the past decades.
Many schemes have been proposed theoretically and achieved experimentally \cite{dicarlo2010preparation, neeley2010generation, delteil2016generation, Hu2020entangledqudits, Erhard2020review}.
However, the preparation of high-dimensional entangled states is so difficult that it has been implemented almost exclusively in photonic systems \cite{krenn2014generation}.
Entanglement swapping \cite{zukowski1993entanglementswapping}, also known as teleportation of entanglement, can entangle two particles that are not related at first by performing a joint Bell-state measurement on them, which plays an important role in long distance quantum communication, such as quantum repeater \cite{briegel1998repeater, santra2019repeater}. It is also the heart for quantum network and entanglement distribution.
However, Ref. \cite{John2002generalizedmeasurment} has pointed out that it is impossible to perform perfect Bell state measurement for qudits by using only linear elements and particle detectors.
It is of interest to know if we can sidestep this thorny issue but still entangle designate parts.
And in addition to the simple case involving two-qubit entangled quantum state, GHZ entanglement swapping by using three pairs EPR entangled state has also been discussed \cite{Lu2009GHZ}. Then, Ref. \cite{Su2016multipartite} demonstrated entanglement swapping between two multipartite entangled state, such as GHZ state and EPR state, and two GHZ states.
However, a relatively general framework for generating entangled state via entanglement swapping is missing, let alone the experimental platform chosen to test.

Motivated by this, we develop a new scheme to generate multipartite entanglement in high dimensions by using quantum walks instead of direct Bell state measurement in the framework of entanglement swapping and thus avoid the difficulty of the realization of Bell state measurement.
Our scheme is mainly based on the model of quantum walks with multiple coins \cite{brun2003quantum} which can bring out entanglement between position space and coin space and has been proved to be useful in quantum communication protocols\cite{wang2017generalized, shang2019pst}.
Here we mainly consider quantum walks with multiple coins on 2-line, 2-circle and 2-complete ($d$-complete) graph, which are all the building blocks in quantum network.
In this way, we can generate two-qubit entangled  state, two-qudit entangled state, three-qubit GHZ state and even three-qudit GHZ state between several designate parties.

In terms of experimental realization of entangled state, it is very difficult to find an appropriate experimental preparation device especially for the high-dimensional entangled state, let alone universal method and experimental device.
However, the physical implementations of the quantum walk have been realized in many different physical systems, such as trapped atom \cite{karski2009trappedatom}, trapped ions \cite{tamura2020trappedions}, photonic system \cite{tang2018photonic} and so on. Since quantum walk is a universal quantum computing model \cite{childs2009universal, childs2013universal, lovett2010universal, underwood2010universal}, these schemes may provide a universal platform for entanglement generation.

Recently, IBM Quantum Experience, as a superconducting flux qubit based open platform, attracts the attention of many researchers \cite{IBM_Q}. Many theoretical schemes in the field of quantum computation and quantum information have been detected and performed on this platform \cite{shang2019experiment, alsina2016experimental, devitt2016experiment, behera2017experimental, Wang2018entangledqubits}.
Here, we carry out experiments on this platform to generate Bell state and three-qubit GHZ state between several designate parties.
Also, the accuracy of our scheme can be verified well by performing quantum state tomography \cite{james2001measurement}.

Furthermore, we give its application in quantum communication. In detail, we provide an improved multiparty quantum secret sharing protocol based on Ref. \cite{zhang2005MQSS} by using our scheme.
Compared with the original protocols, our new scheme are high code capacity and high security. And it is easy to be implemented in experiment.

\section{Results}
\begin{figure}[htb]
 \centering
 \subfigure[]{\label{schematic_diagram1}
 \includegraphics[width=5cm]{schematic_diagram1.pdf}}\ \ \
 \subfigure[]{\label{schematic_diagram2}
 \includegraphics[width=5.5cm]{schematic_diagram2.pdf}}
 \caption{The schematic diagram of our scheme. The small boxes indicate particles and the dotted lines indicate entanglement. (a) represents the generation of bipartite entangled state. (b) represents the generation of GHZ state.
 Due to the symmetry of the particles, any two or three unrelated particles can be entangled by our scheme.
But for the sake of narration, let us take entanglement generation of particles $1, 4$ and particles $1, 4, 5$ as examples.}
 \label{diagram_scheme}
 \end{figure}

For the quantum walks with $k$ coins on an arbitrary graph, the $k$ coins were used cyclically.
At $i$-th step of quantum walks, the coin state can be flipped with the $i$-th coin operator $C_{i}$ and then the walker moves to another position according to the conditional shift operator $S_{i}$, which can be described as
\begin{equation}
U_{i}=(S_{i}\otimes I)\cdot(I\otimes C_{i}).
\end{equation}
After $t$ steps, the initial state $|\phi(0)\rangle$ evolves into
\begin{equation}
|\phi(t)\rangle=(U_{k}\cdots U_{1})^{t/k}|\phi(0)\rangle.
\end{equation}
We denote the line, circle and complete graph with $N$ vertices as $N$-line, $N$-circle and $N$-complete graph respectively. The three corresponding conditional shift operators acting on the position space and the coin space that is active at the step can be written as:
\begin{equation} \label{Ncomplete}
S^{N-complete}=\sum_{x,j=0}^{N-1}|(x+j)\bmod\,N\rangle\langle x|\otimes|j\rangle\langle j|,
\end{equation}
\begin{equation} \label{Ncycle}
S^{N-circle}=\sum_{x=0}^{N-1}(|(x+1)\bmod\,N\rangle\langle x|\otimes|0\rangle\langle0|+|(x-1)\bmod\,N\rangle\langle x|\otimes|1\rangle\langle1|),
\end{equation}
\begin{equation} \label{Nline}
\begin{split}
S^{N-line}&=\sum_{x=0}^{N-2}|(x+1)\bmod\,N\rangle\langle x|\otimes|0\rangle\langle0|+\sum_{x=1}^{N-1}|(x-1)\bmod\,N\rangle\langle x|\otimes|1\rangle\langle1| \\
&+|N-1\rangle\langle N-1|\otimes|1\rangle\langle0|+|0\rangle\langle0|\otimes|0\rangle\langle1|.
\end{split}
\end{equation}
In particular, when the number of vertices is $2$, the conditional shift operators are:
\begin{equation} \label{2complete}
S^{2-complete}=(|0\rangle\langle0|+|1\rangle\langle1|)\otimes|0\rangle\langle0|
              +(|1\rangle\langle0|+|0\rangle\langle1|)\otimes|1\rangle\langle1|,
\end{equation}
\begin{equation} \label{2cycle}
S^{2-circle}=(|1\rangle\langle0|+|0\rangle\langle1|)\otimes|0\rangle\langle0|
           +(|1\rangle\langle0|+|0\rangle\langle1|)\otimes|1\rangle\langle1|,
\end{equation}
\begin{equation} \label{2line}
S^{2-line}=|1\rangle\langle0|\otimes|0\rangle\langle0|+|0\rangle\langle1|\otimes|1\rangle\langle1|
          +|0\rangle\langle0|\otimes|0\rangle\langle1|+|1\rangle\langle1|\otimes|1\rangle\langle0|.
\end{equation}

\subsection{Generation of two-qubit entangled state}
We start by discussing the generation of two-qubit entangled state. And the experimental setup is illustrated in Figure \ref{schematic_diagram1}.
Without loss of generality, we can assume that the initial state of the four particles is
\begin{equation} \label{qubit_initial}
|\psi(0)\rangle=(a|01\rangle+b|10\rangle)_{1,2}(a|01\rangle+b|10\rangle)_{3,4},
\end{equation}
where $|a|^{2}+|b|^{2}=1$.
Now we try to generate entanglement between particle $1$ and $4$.
Here, we can view particle $1$ and $2, 3, 4$ as walker and three coins respectively.
And thus we can consider quantum walks with three coins on $2$-line, $2$-circle and $2$-complete graph.

For the case of quantum walks with three coins on $2$-line, we can perform two-step quantum walks where $C_{1}=H$ and $C_{2}=I$. According to conditional shift operator and initial state shown in equation (\ref{2line}) and (\ref{qubit_initial}), the evolution process is written as follows:
\begin{equation}
\begin{split}
|\psi(1)\rangle=&\frac{1}{\sqrt{2}}(a^{2}|1001\rangle-a^{2}|0001\rangle+ab|1010\rangle-ab|0010\rangle\\
&+ab|1101\rangle+ab|0101\rangle+b^{2}|1110\rangle+b^{2}|0110\rangle)_{1,2,3,4},
\end{split}
\end{equation}
\begin{equation}
|\psi(2)\rangle=\frac{1}{\sqrt{2}}(-a|00\rangle+a|01\rangle+b|10\rangle+b|11\rangle)_{2,3}(a|11\rangle+b|00\rangle)_{1,4}.
\end{equation}
Therefore, we can definitely obtain entangled state between particle $1$ and $4$ by performing two local measurements on the other two particles. In particular, the state will be Bell state, also known as maximal entangled state, when $a=b=\frac{1}{\sqrt{2}}$.

In the case of quantum walks with three coins on $2$-circle, it is impossible to get entangled between particle $1$ and $4$ no matter how many steps we take. In this case, the conditional shift operator happens to be the second-order Pauli matrix $X_{2}$.
Therefore, the general form of evolved quantum state under three-step quantum walks is $(I\otimes C_{1}\otimes C_{2}\otimes C_{3})(a^{2}|1101\rangle+ab|1110\rangle+ab|0001\rangle+b^{2}|0010\rangle)$.
And it has a definite explicit form by giving three unitary coin operators with parameters.
Then, it is not hard for us to find that there is no satisfied coin operator by analyzing the final state with parameters.

As regards quantum walks with three coins on $2$-complete graph, we can take three-step quantum walks where $C_{1}=I$, $C_{2}=H$ and $C_{3}=X_{2}$. Thus
the whole evolutionary process can be described as follows:
\begin{equation}
|\psi(1)\rangle=(a^{2}|1101\rangle+ab|1110\rangle+ab|1001\rangle+b^{2}|1010\rangle)_{1,2,3,4},
\end{equation}
\begin{equation}
\begin{split}
|\psi(2)\rangle=&\frac{1}{\sqrt{2}}(a^{2}|1101\rangle+a^{2}|0111\rangle+ab|1100\rangle-ab|0110\rangle\\
&+ab|1001\rangle+ab|0011\rangle+b^{2}|1000\rangle-b^{2}|0010\rangle)_{1,2,3,4},
\end{split}
\end{equation}
\begin{equation} \label{final_state_of_qubit_on_complete}
|\psi(3)\rangle=\frac{1}{\sqrt{2}}[(b|00\rangle+a|10\rangle)_{2,3}(a|10\rangle+b|01\rangle)_{1,4}+(b|01\rangle+a|11\rangle)_{2,3}(a|00\rangle-b|11\rangle)_{1,4}].
\end{equation}
So we can obtain entanglement between particle $1$ and $4$ by performing three-step quantum walks and two local measurements.
Of course, the state obtained happens to be Bell state when $a=b=\frac{1}{\sqrt{2}}$.
In addition, there are other ways to select coin operators to generate new entanglement. For example, $C_{1}=I$, $C_{2}=X_{2}$, $C_{3}=H$ or $C_{1}=X_{2}$, $C_{2}=X_{2}$, $C_{3}=H$, or $C_{1}=X_{2}$, $C_{2}=H$, $C_{3}=X_{2}$ can also be used to generate entanglement between particle $1$ and $4$.

\subsection{Generation of two-qudit entangled state}
High-dimensional entangled states are of paramount interest both from a theoretical and practical perspective.
Generalized Bell states in Hilbert space $C^{d}\otimes C^{d}$, a basis of maximally entangled high-dimensional bipartite states, can be described as \cite{bennett1993teleporting, sych2009generalizedBellstates}
\begin{equation}
|\psi_{k,l}\rangle=\frac{1}{\sqrt{d}}\sum_{m=0}^{d-1}\exp(\frac{2\pi \mathrm{i}}{d}mk)|m\rangle|m-l\rangle,
\end{equation}
where $m-l$ means $(m-l)\ \mathrm{mod}\ d$.
The general entangled state can be written as $\sum_{i=0}^{d-1}a_{i}|i\rangle|i\rangle$, where $\sum_{i=0}^{d-1}|a_{i}|^{2}=1$.
For the sake of simplicity, we denote $d$-order generalized Pauli operator as $X_{d}=\sum_{i=0}^{d-1}|i+1\rangle\langle i|$, which is also known as discrete Weyl operator.
Next, we propose two novel schemes to generate maximal and non-maximal entangled high-dimensional states (two-qudit entangled state) by performing quantum walks with three coins on $d$-complete graph.  And its schematic diagram is shown in Figure \ref{schematic_diagram1}.

For the case of maximal entangled state, the initial state is $|\psi(0)\rangle=|\psi_{k,l}\rangle_{1,2}|\psi_{k,l}\rangle_{3,4}=\frac{1}{d}\sum_{m,n=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}(m+n)k)|m,m-l,n,n-l\rangle$.
Now we can perform three-step quantum walks where $C_{1}=I$, $C_{2}=F$ (Fourier transformaiton) and $C_{3}=X_{d}$.
Thus, we can obtain the final state:
\begin{equation}
|\psi(3)\rangle=\frac{1}{d}\sum_{m,p=0}^{d-1}[\exp(\frac{2\pi\mathrm{i}}{d}mk)|m-l\rangle_{2}|p\rangle_{3}
(\frac{1}{\sqrt{d}}\sum_{n=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}(k+p)n)|2m+p+n-2l+1\rangle_{1}|n-l+1\rangle_{4})].
\end{equation}
Then we take two local measurements on particle $2$ and $3$ and assume that the result are $\exp(\frac{2\pi\mathrm{i}}{d}m_{0}k)|m_{0}-l\rangle_{2}|p_{0}\rangle_{3}$. So the reduced state would be $\exp(-\frac{2\pi\mathrm{i}}{d}(2m_{0}+p_{0}-2l+1)(k+p_{0}))|\psi_{k+p_{0},2m_{0}+p_{0}-l}\rangle_{1,4}$.

More generally, in the case of non-maximal entangled state, the initial state can be written as $|\psi(0)\rangle=(\sum_{i=0}^{d-1}a_{i}|i\rangle|i\rangle)_{1,2}\otimes(\sum_{j=0}^{d-1}b_{j}|j\rangle|j\rangle)_{3,4}$, where $\sum_{i=0}^{d-1}|a_{i}|^{2}=\sum_{j=0}^{d-1}|b_{j}|^{2}=1$.
After the same three-step quantum walks ($C_{1}=I$, $C_{2}=F$ and $C_{3}=X_{d}$), the final state will be
$|\psi(3)\rangle=\frac{1}{\sqrt{d}}\sum_{i,k=0}^{d-1}[a_{i}|i\rangle_{2}|k\rangle_{3}(\sum_{j=0}^{d-1}b_{j}\exp(\frac{2\pi\mathrm{i}}{d}jk)|2i+j+k+1\rangle_{1}|j+1\rangle_{4})]$.
Similarly, particle $1$ is entangled with $4$ after two local measurements.

\subsection{Generation of three-qubit GHZ state}
The three-qubit GHZ states  maximize entanglement monotones \cite{normalform2003, threequbit2000} and have extremely non-classical properties.
The setup for this problem is shown in the Figure \ref{schematic_diagram2}.
The initial state of the five particles can be written as
\begin{equation} \label{GHZ_initial}
|\psi(0)\rangle=(a|01\rangle+b|10\rangle)_{1,2}(\frac{|000\rangle+|111\rangle)}{\sqrt{2}})_{3,4,5},
\end{equation}
where $|a|^{2}+|b|^{2}=1$.
In order to generate entanglement between particle $1$, $4$ and $5$, we can view these five particles as the walker and four coins respectively.
So we can discuss the quantum walks with four coins on $2$-line, $2$-circle and $2$-complete graph.

Let us start by talking about the case of quantum walks with four coins on $2$-line. And in order to do that, we perform two-step quantum walks where $C_{1}=H$ and $C_{2}=X$. Thus the initial quantum state presented in equation (\ref{GHZ_initial}) will evolves over time into
\begin{equation}
|\psi(1)\rangle=\frac{1}{2}(-a|00\rangle+b|01\rangle+a|10\rangle+b|11\rangle)_{1,2}(|000\rangle+|111\rangle)_{3,4,5},
\end{equation}
and
\begin{equation}
|\psi(2)\rangle=\frac{1}{2}(-a|00\rangle+a|01\rangle+b|10\rangle+b|11\rangle)_{2,3}(|000\rangle+|111\rangle)_{1,4,5}.
\end{equation}
Then we just have to make local measurements of particle $2$ and particle $3$, and no matter what the measurements are, the remaining three particles always collapses to GHZ state.

For quantum walks with four coins on $2$-circle, particle $1$, $4$ and $5$ will never get entangled no matter how many steps we take.
Given an initial state in equation (\ref{GHZ_initial}), after four steps the state evolves to $\frac{1}{\sqrt{2}}(I\otimes C_{1}\otimes C_{2}\otimes C_{3}\otimes C_{4})(a|01000\rangle+a|01111\rangle+b|10000\rangle+b|10111\rangle)$, where these four operators can be any unitary operator of the second order.
By calculating the undetermined coefficients of these four coin operators, we find that particle $1$, $4$ and $5$ are never completely entangled.

In regard to the quantum walks with four coins on $2$-complete graph, we can generate entanglement between particle 1, 4 and 5 by performing four-step quantum walks where $C_{1}=C_{2}=I$ and $C_{3}=C_{4}=H$.
After the first two steps of quantum walks, the quantum state evolves into
\begin{equation}
|\psi(2)\rangle=\frac{1}{\sqrt{2}}(a|11000\rangle+a|01111\rangle+b|10000\rangle+b|00111\rangle)_{1,2,3,4,5}.
\end{equation}
The resulting state will be
\begin{equation}
\begin{split}
|\psi(3)\rangle=&\frac{1}{2}(a|11000\rangle+a|01010\rangle+a|01101\rangle-a|11111\rangle\\
&+b|10000\rangle+b|00010\rangle+b|00101\rangle-b|10111\rangle)_{1,2,3,4,5}
\end{split}
\end{equation}
after one step with coin operator $C_{3}=H$.
And then the final quantum state can be described as follows:
\begin{equation}
\begin{split}
|\psi(4)\rangle=\frac{1}{2\sqrt{2}}[&(b|00\rangle+a|10\rangle)_{2,3}(|100\rangle+|001\rangle+|010\rangle+|111\rangle)_{1,4,5}\\
&+(b|01\rangle+a|11\rangle)_{2,3}(|000\rangle-|101\rangle-|110\rangle+|011\rangle)_{1,4,5}].
\end{split}
\end{equation}
So we can get the GHZ-like state \cite{yang2009GHZlike} of particle $1$, $4$ and $5$, which is the entangled quantum state and belongs to the class of GHZ state, by taking two local measurements of particle $2$ and $3$.

\subsection{Generation of three-qudit GHZ state}
Three-qudit GHZ state can be written as $\frac{1}{\sqrt{d}}\sum_{n=0}^{d-1}|n\rangle^{\otimes^{3}}$.
And its generation process is presented in Figure \ref{schematic_diagram2}.
For the initial state
\begin{equation}
|\psi(0)\rangle=
|\psi_{k,l}\rangle_{1,2}(\frac{1}{\sqrt{d}}\sum_{n=0}^{d-1}|n\rangle^{\otimes^{3}})_{3,4,5}
=\frac{1}{d}\sum_{m,n=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}mk)|m,m-l,n,n,n\rangle_{1,2,3,4,5},
\end{equation}
we can perform four-step quantum walks with four coins on $d$-complete graph by choosing $C_{1}=C_{2}=I$ and $C_{3}=C_{4}=F$ to establish the entanglement between particle 1, 4 and 5 that are unrelated at first.
The evolutionary process can be described as follows:
\begin{equation}
|\psi(2)\rangle=\frac{1}{d}\sum_{m,n=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}mk)|2m-l+n,m-l,n,n,n\rangle_{1,2,3,4,5},
\end{equation}
\begin{equation}
\begin{split}
|\psi(4)\rangle&=\frac{1}{d^{2}}\sum_{m,n,p,q=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}(mk+np+nq))|2m-l+n+p+q,m-l,n,p,q\rangle_{1,2,3,4,5}\\
&=\frac{1}{d}\sum_{m,n=0}^{d-1}[\exp(\frac{2\pi\mathrm{i}}{d}mk)|m-l\rangle_{2}|n\rangle_{3}
(\frac{1}{d}\sum_{p,q=0}^{d-1}\exp(\frac{2\pi\mathrm{i}}{d}(np+nq))|2m-l+n+p+q,p,q\rangle_{1,4,5})].
\end{split}
\end{equation}
Then we can obtain an entangled state of particle $1$, $4$ and $5$ by performing two local measurements on particle $2$ and $3$.
Note that any state like $\sum_{p,q=0}^{d-1}|x_{0}+p+q,p,q\rangle$ is completely entangled, where $x_{0}\in\{0,1,2,\cdots,d-1\}$. It can be seen as the high-dimensional case of GHZ-like state.

\subsection{Experimental protocol}
We provide the experimental realization of our scheme using IBM quantum computer $ibmq_{\_}5_{\_}yorktown-ibmqx2$ (it is called $ibmqx2$ for short) and $simulator$ with 8192 number shots (see Supplementary Note 1 for the introduction). Here, we mainly discuss the generation of standard Bell state and three-qubit GHZ state based on quantum walks.

\subsubsection{Experimental generation of Bell state on $2$-complete graph}
\begin{figure}[htb]
 \centering
 \includegraphics[width=9cm]{ES_Bell_circuit.pdf}
 \caption{Quantum circuit for generating Bell state between two parties that are not related at first. The measurements on $q[0]$ and $q[3]$ are in X-basis. The measurements on $q[1]$ and $q[2]$ are in Z-basis.}
 \label{ES_Bell_circuit}
 \end{figure}
\begin{figure}[htb]
 \centering
 \subfigure[]{\label{ES_Bell_result_simulator}
 \includegraphics[width=6cm]{ES_Bell_simulator.pdf}}\ \
 \subfigure[]{\label{ES_Bell_result_ibmqx2}
 \includegraphics[width=6cm]{ES_Bell_ibmqx2.pdf}}
 \caption{The probability results of quantum circuit after four measurements in Z-basis. (a) is on $simulator$ and (b) is on $ibmqx2$.}
 \label{ES_Bell_result}
 \end{figure}
Based on the theoretical scheme, we provide the corresponding IBM quantum circuit that can be seen in Figure \ref{ES_Bell_circuit} to generate Bell state by performing three-step quantum walks on $2$-complete graph.
And the results, measured in Z-basis, are shown in Figure \ref{ES_Bell_result}.
According to the Figure \ref{ES_Bell_result_simulator}, the simulation seems very close to the theoretical result (there are eight outcomes appeared in Figure \ref{ES_Bell_result_simulator} and the probability of each is $12.5\%$). And the results on $ibmqx2$ shown in Figure \ref{ES_Bell_result_ibmqx2} are not as good as the simulation results. In fact, there are some noises on the quantum platform, such as decoherence, depolarizing, general noises and so on, which result in the discrepancy between theoretical values and the run results.
According to the quantum state tomography \cite{james2001measurement} and the detailed experimental data, we find that the fidelity between theoretical and simulation density matrix $F(\rho_{ij}^{T},\rho_{ij}^{E\_simulator})$ where $i,j\in\{0,1\}$ are approximately equal to 1 (see Supplementary Note 2).
And the fidelity between theoretical and experimental density matrix $F(\rho_{ij}^{T},\rho_{ij}^{E\_ibmqx2})$ are 0.8535, 0.8793, 0.7909, 0.8549 (see Supplementary Note 2).

\subsubsection{Experimental generation of three-qubit GHZ state on $2$-line}

\begin{figure}[htb]
 \centering
 \subfigure[]{\label{ES_GHZ_circuit}
 \includegraphics[width=7cm]{ES_GHZ_circuit.pdf}}\ \
 \subfigure[]{\label{ES_GHZ_visualization}
 \includegraphics[width=7cm]{ES_GHZ_visualization.pdf}}\ \
 \subfigure[]{\label{ES_GHZ_ibmqx2}
 \includegraphics[width=11cm]{ES_GHZ_ibmqx2.pdf}}
 \caption{(a): Quantum circuit for generating three-qubit GHZ state between three parties that are not related at first. (b): The statevector amplitude results of quantum circuit after five measurements in Z-basis.(c): The probability results of quantum circuit on \textit{ibmqx2}.}
 \label{ES_GHZ}
 \end{figure}

As for the generation of three-qubit GHZ state, we provide the corresponding IBM quantum circuit diagram that can be seen in Figure \ref{ES_GHZ_circuit}, which is realized by performing two-step quantum walks on $2$-line on the basis of theoretical scheme.
The whole design idea is similar to the experimental realization of Bell state.
The only thing that needs to be paid attention to is the design of shift operator on $2$-line.
Here we give the amplitude results and probability results on \textit{ibmqx2} in Figure \ref{ES_GHZ_visualization} and \ref{ES_GHZ_ibmqx2} respectively.
The analysis of the experimental results can be similarly obtained, so we will not repeat it.

\section{Application}
In order to share secret among some agents and reconstruct it only when they collaborate together, the scheme of quantum secret sharing (QSS) was first proposed in 1999 \cite{mark1999}.
Zhang and Man proposed a multiparty quantum secret sharing (MQSS) protocol by using entanglement swapping \cite{zhang2005MQSS}.
Then Lin and Zhang revised and improved the protocol respectively \cite{lin2007MQSS,zhang2007reply}.

Here, we improve and perfect the performance of MQSS protocol proposed by Zhang and Man \cite{zhang2005MQSS} by using our scheme.
According to our scheme, for the initial state
\begin{equation}
|\phi(0)\rangle=|\psi_{x,y}\rangle_{1,2}|\psi_{k,l}\rangle_{3,4}=\frac{1}{d}\sum_{m,n=0}^{d-1}\exp(\frac{2\pi \mathrm{i}}{d}mx+nk)|m,m-y,n,n-l\rangle_{1,2,3,4},
\end{equation}
after three-step quantum walks (the coin operators are $I$, $F$, $X_{d}$ in sequence), we could obtain the final state
\begin{equation} \label{application}
\begin{split}
|\phi(3)\rangle&=\frac{1}{d\sqrt{d}}\sum_{m,n,p=0}^{d-1}\exp(\frac{2\pi \mathrm{i}}{d}mx+nk+np)|2m+p+n-y-l+1,m-y,p,n-l+1\rangle_{1,2,3,4} \\
&=\frac{1}{d}\sum_{m,p=0}^{d-1}\exp(\frac{2\pi \mathrm{i}}{d}(m x-(2m+p-y-l+1)(k+p)))|m-y\rangle_{2}|p\rangle_{3}|\psi_{k+p,2m+p-y}\rangle_{1,4}.
\end{split}
\end{equation}
Assume that $d$ is odd, we can set up a bijection between the pairs $(m-y, p)$ and $(k+p, 2m+p-y)$, where $y$ and $k$ are known fixed numbers.

Let $U_{i}=X_{d} ^{i}=\sum_{j=0}^{d-1}|j+i\rangle\langle j|$ be $d$ local unitary transformations ($i=0,1,\cdots,d-1$) that correspond to $\lceil\log d\rceil$ classical bits respectively. To be specific, $U_{0}$ corresponds to $``0\cdots000"$, $U_{1}$ corresponds to $``0\cdots001"$, $U_{2}$ corresponds to $``0\cdots010"$ and so on.
Based on the property $U_{i}$ holds,
\begin{equation}
(I\otimes U_{i})|\psi_{k,l}\rangle=|\psi_{k,l-i}\rangle,
\end{equation}
we can establish correspondence between the unitary operators $\{U_{i}: i=0,1,\cdots,d-1\}$ and generalized Bell states $\{|\psi_{k,l}\rangle: l=0,1,\cdots,d-1\}$ for an arbitrary fixed $k$.
Next, we will show the improved MQSS based on \cite{zhang2005MQSS} in which $d$ is an odd number. And the schematic flow chart of the whole process is shown in Figure \ref{diagram_application}.
\begin{figure}[htb]
 \centering
 \includegraphics[width=10cm]{schematic_diagram3.pdf}
 \caption{The schematic diagram of the improved multiparty quantum secret sharing where the small boxes indicate particles and the dotted lines indicate entanglement. The particles in yellow belong to Alice. The particles in green and red are in Bob's and Charlie's possession. The purple particle has been performed by a unitary operator.}
 \label{diagram_application}
 \end{figure}

\begin{enumerate}
  \item Alice prepares three identical generalized Bell state pairs $|\psi_{k,l}\rangle_{1,2}$, $|\psi_{k,l}\rangle_{3,4}$ and $|\psi_{k,l}\rangle_{5,6}$. She sends qudit 2 and 3 to Bob, 4 and 5 to Charlie via their quantum channels respectively. Then, Alice enters step 2 with probability $q$, or enters step 3 with probability $1-q$.
  \item Alice chooses randomly from two mutually unbiased measurement basis, $\mathcal{M}_{1}=\{|i\rangle: i=0,1,\cdots,d-1\}$ and $\mathcal{M}_{2}=\{|\tilde{i}\rangle: i=0,1,\cdots,d-1\}$, where $|\tilde{i}\rangle=\frac{1}{\sqrt{d}}\sum_{j=0}^{d-1}\exp(\frac{2\pi \mathrm{i}}{d}ij)|j\rangle$. Alice tells Bob the basis she has chosen. They measure qudit 1 and 2 using this chosen measurement basis respectively. Then Bob tells Alice his outcome of qudit 2. Alice compares these two outcomes to judge the security of Alice-Bob quantum channel (For example, if the generalized Bell state pairs that Alice prepares is $|\psi_{0,0}\rangle$, these two measurement outcomes must be identical when the channel is safe). The security of Alice-Charlie quantum channel can also be checked in the same way. If the two channel are both secure, returns to step 1. Otherwise, this process will be terminated.
  \item Alice applies a local unitary operator chosen randomly from $\{U_{i}: i=0,1,\cdots,d-1\}$ on one of her qudits 1 and 6 ( say, on qudit 6). Firstly, they perform three-step quantum walks on qudit 5, 6, 1 and 2. Alice performs two local measurements on qudit 6 and 1 and announces her outcome. Next, they perform three-step quantum walks on qudit 5, 2, 3 and 4. Bob performs two local measurements on qudit 2 and 3 and writes down the outcomes. At this times, Charlie records his own state. Bob and Charlie can deduce the unitary transformation Alice performs and the classical bits Alice wants to share.
\end{enumerate}

Now let us explain how we derived $U_{i}$ in step 3 based on the three measurement results they announced and evolutionary properties.
Without generality, we can denote the outcome Alice announces as $|a,b\rangle_{6,1}$.
Suppose the quantum state of qudit 5 and 2 is $|\psi_{\bar{k}, \bar{l}}\rangle$ after the first quantum walks.
According to equation (\ref{application}), we can denote the outcome Bob writes down and the state Charlie possesses as $|m-\bar{l},p\rangle_{2,3}$ and $|\psi_{k+p,2m+p-\bar{l}}\rangle_{5,4}$ respectively.
Since $k$ is known, $p$, $m$ and $\bar{l}$ are also known.
Next we can infer that $\bar{k}=k+p$ and the state of qudit 5 and 6 after the local unitary transformation performed by Alice is $|\psi_{k,\bar{l}-b-2a}\rangle_{5,6}$ based on the first quantum walks.
Thus, the unitary transformation Alice performed is $U_{l-(\bar{l}-b-2a)}$. Note that the above mentioned addition and subtraction methods are all of modulo $d$. So it is not hard to infer the classical bits Alice wants to share.

In addition, our protocol can also be generalized to a multiparty case. Without loss of generality, we suppose that there are $N$ parties in total.
First, the sender prepares $N$ identical generalized Bell state pairs and distributes two qudits to the agents one by one and does not tell their the ordering of the two qudits. So the $i$-th agent and the $i+1$-th agent share a generalized Bell state pair. And the $N$-th agent will share a generalized Bell state pair with the sender.
Then, they check the security of quantum channel (step 2) with probability $q$ or transmit message (step 3) with probability $1-q$. The order of quantum walk and measurements, which is crucial for the whole process, is the same as the order of entanglement swapping in the previous protocol \cite{lin2007MQSS, zhang2007reply}.
The security (even in a noisy channel) can be guaranteed by the discussions in Ref. \cite{zhang2005MQSS, lin2007MQSS, zhang2007reply}.

Now, let us make some comparisons between our scheme and the other schemes \cite{zhang2005MQSS, zhang2007qutrit}.
Compared with Ref. \cite{zhang2005MQSS}, the advantage of our scheme is that the capacity of classical messages it can distribute is much larger.
Specifically, it is a task to distribute $N=\lceil\log d\rceil$ classical bits instead of just two bits.
In addition, Ref. \cite{zhang2005MQSS} needs $\frac{3N}{2(1-q)}$ Bell states in step 3. Our scheme just needs $\frac{3}{1-q}$ generalized Bell states. Thus, as for the resource, our scheme is much better.
In Ref. \cite{zhang2007qutrit}, Zhang et.al designed a MQSS protocol by swapping qutrit-state entanglement and then generalized to the qudit case, in which the so-called two-qutrit (two-qudit) entangled-state measurement is crucial but hard to implement in experiment.
Here, we mainly use the model of quantum walks with multiple coins which is practical and novel.
Therefore, our scheme avoids the difficulty in the experimental realization of Bell state measurement, especially for the measurement of high dimensional Bell state \cite{John2002generalizedmeasurment}.
So our protocol is more feasible than the protocol in Ref. \cite{zhang2007qutrit}.
Also, it may inspire more potential applications in quantum communication and quantum information.

\section{Discussion}
The generation of entanglement, especially in the high-dimensional case, between several designated parties is vital and essential in scalable quantum network.
Inspired by entanglement swapping, we have developed a scheme to generate entangled state by using the model of quantum walks with multiple coins, including two-qubit entangled state, two-qudit entangled state, three-qubit GHZ state and three-qudit GHZ state.
In comparison with original entanglement swapping, the benefit of using our scheme is expected to avoid performing joint Bell state measurement that is still an unsolved problem to date, and provide a relatively general framework for generating entangled state via entanglement swapping.

Also, we provide the experimental realization of our scheme to generate Bell state and three-qubit GHZ state on IBM platform.
And our scheme exhibits high fidelity by performing quantum state tomography.
Moreover, we present an improved multiparty quantum secret sharing by using our our scheme.
After our improvement, it can convey more classical information and require less quantum resources (Bell states).

It remains an interesting open question to design further the GHZ state and $W$ state for a system of $n$ qubit, where $n$ is an arbitrary integer.
It is also highly desired to extend our method to more practical applications.

\section{Data availability}
The experimental data that supports our finding are available from the corresponding author on reasonable request.

\section{Acknowledgements}
We thank the support of National Key Research and Development Program of China under grant 2016YFB1000902, National Natural Science Foundation of China (Grant No.61472412, 61872352), and Program for Creative Research Group of National Natural Science Foundation of China (Grant No. 61621003).

\section{References}


\section{System Model} \label{system_model}
	We are considering the device-tier of 5G D2D architecture mentioned in \cite{tehrani2014device}, where devices or UEs can communicate among themselves with or without the help from BS. The service region is discretized into various \textit{zones} or \textit{grids} as shown in figure \ref{fig_grid_time}(a) with one BS. Each zone may have many UEs and is assumed to have at least one D2D device which is ready to take part in D2D communication as a relay or source/destination node.
	A zone $i$ containing the source UE may form connection to a UE of another zone $j\in \mathbb{U}^i$ , where $\mathbb{U}^i$ is the \textit{viable} relay zones of the zone $i$ which is given by the BS. A viable relaying zone of zone $i$ is one which is nearer to the zone containing the destination UE and is in the communication range of the zone $i$. When the source UE in the fixed zone $i$ forms a connection with another UE of zone $j\in \mathbb{U}^i$, then the link formed between zones $i$ and $j$ is termed as link $j$.	Link is formed between UEs of two zones when they are in communication range of each other and the received signal strength is sufficient for the required data rate.
	Each UE can communicate with one another on mmWave channels using directional antennas. The received signal strength ($Q_{ij}$) on zone $j$ from zone $i$  is modeled as \cite{p14_7450161}:
	\begin{equation}
	Q_{ij}=\mu \cdot P \cdot G_t\cdot G_r \cdot PL_{ij}
	\end{equation}
	where, $\mu$ is the shadowing random variable, $P$ is the transmit power of the UE, $G_t$ \& $G_r$ are transmit and receive beam-forming gains respectively. $PL_{ij}$ is the distance dependent path loss function between zone $i$ and zone $j$.

	Time is discretized  as $(nN+l)\delta$ as shown in figure \ref{fig_grid_time}(b), where $n$ belongs to set of nonnegative integers, $l$ takes integer values in $[0,N-1]$, $\delta$ is the smaller discretized time slot when the UEs transmit packets locally. It is assumed that $\delta$ (for each $l\in[0,N-1]$) is large enough to send one packet of size $L$ bytes. Here, $N$ is the number of time slots (of $\delta$ duration) between two consecutive global decisions by the BS. Global decision by BS is made at time when $nN+l$ is divisible by $N$. At this time instant BS takes the channel state information from all UEs in the service region and gives the decision of best relaying UE of a given zone for a given source UE. Hence, in between two consecutive time instants when BS can make global decision, a UE can send at-most $N$ packets of size $L$ to another UE. Note that at time $l=0$, the UE chooses the relay link suggested by the BS and at time  $l\in \{1,2,\cdots,N-1\}$, UEs do not get channel state information from the BS.  At global time instants, BS sends two types of information to UEs, i) the best relaying UE for a given source UE and ii) viable relaying zones $\mathbb{U}^i$ for given source zone $i$, hence the zone $i$ may choose  an appropriate zone for relaying data from the set $\mathbb{U}^i$ by undergoing exploration.
	\begin{figure}[h!]
		\centering
		\includegraphics[width=0.44\textwidth]{grid_time.eps}
		\caption{(a) Service region divided into zones along with dynamic obstacles. (b) Discretized time slots with exploration time unit.}
		\label{fig_grid_time}
	\end{figure}

	\textbf{Exploration}: When the current link quality is not good enough then source UE locally explores for an appropriate relaying UE from given viable set. One of the basic procedure in the exploration phase involves sequentially searching space in all directions to align transmitter and receiver beams (beam alignment). The quality of link is piggy-back to the sending UE. This overall process of beam alignment is termed as the exploration phase. Any of the state of art approaches mentioned in \cite{8458146} can be applied for the exploration phase. We have considered an abstraction of this phase and thus it is considered as a black box. It is assumed that the beams are perfectly aligned after the exploration phase and thus we focus on the effects of dynamic obstacle's blockage on a given D2D link. This exploration by source UE in zone $i$ is done for UEs belonging only in the set $\mathbb{U}^i$ to find out the best relaying zone for that time instant. Note that the UE is using directional mmWave antennas for exploring the neighbors and hence explorations cause some significant delay. Each exploration time duration is denoted by $\epsilon < \delta$ as shown in figure \ref{fig_grid_time}(b). In a given duration $\delta$, we can perform a number of explorations. Each of the time instants during exploration time of duration $\epsilon$  is denoted by $m\in\{0,1,\cdots,V-1\}$, where $V=\floor*{\delta/\epsilon}$.
	The overall exploration time is bounded by the maximum value $M$ which can be more than $V$. It is assumed that once exploration is complete,  switching takes negligible time.

	There are static and dynamic obstacles in the service region. There is \textit{no} external facility like radars or vision cameras available at BS to track them. The behavior of dynamic obstacles  are not known a priori and need to be learned from the received ACKs of probe packets during exploration phase in an online fashion. Since mmWaves are highly susceptible to obstacles and suffer from severe penetration losses, we assume that even a single moving or static obstacle may break an already established D2D link and can cause packet loss.

	\section{Problem Formulation as POMDP} \label{problem_formulation}
	Zone $i$ containing source UE forms a link with an UE (relay or destination) of some other zone $j\in \mathbb{U}^i$. Global decision for the best relay is given by the BS at the time instant $nN\delta$ to relay data packet till $(n+1)N\delta$ time instant. The quality of link given by the BS may deteriorate due to presence of dynamic obstacles resulting in link outage and a new relay link needs to be explored from the viable set $\mathbb{U}^i$ locally. A newly found relay link might provide very good data rate initially during exploration phase, but it may gradually deteriorate in future data transmission time because of dynamic obstacles. This problem is subtle due to the fact that exploration time is very small (ranging from few microseconds to $10~ms$ \cite{7744807}) in contrast to the data transmission time which is in order of seconds (usually 100-1000 times more than search time). The main idea is to reduce the E2E delay on the cost of sending additional probe packets till decision of whether selecting or not selecting a given relay for data transmission considering dynamic obstacles can be learned. This process is limited upto  maximum $M$ exploration time units where each exploration time duration is of $\epsilon$ time. If the given relay is not selected within this time then the exploration process again starts for another possible relaying UE in a zone belonging to $\mathbb{U}^i$.

	We model this problem of relay selection during  exploration as that of POMDP considering uncertainty in relay link quality. We will define the POMDP as follows: For all the possible links $j\in \mathbb{U}^i$, the state is written as $y_{m}^j \in\{0,1\}$ for exploration time instant $m$. Values $1$ and $0$ signify if relay link $j$ is in  good or bad state respectively. Good and bad notion denotes that whether link will be formed successfully or not respectively. State represented by $G$ and  $\overline{G}$ respectively denote the good and bad state as shown in figure \ref{fig_prob_2}. The action set is defined as \{do not select relay link,  select relay link, decision cannot be made thus continue sending probe packets\}.  The first action signifies that the current link is bad, hence stop exploring it and start a new exploration. The second action signifies that the current link being explored is of good quality, hence choose it for data transmission. The third action signifies that the decision of relay selection cannot be made and exploration should continue for another $\epsilon$ unit of time to learn the relay link quality. The action is denoted by $a_m^j$ for link $j$ at time instant $m$. The link quality is observed by the ACKs of the received probe packet signal. The  ACK test for link $j$ at time $m$ is denoted as $w_{m}^j\in\{0,1\}$ for a given state $y_m^j$ and action $a_m^j$, may also be uncertain due to the unpredictable channel condition.  In figure \ref{fig_prob_2}, $A$ and $\overline{A}$  denotes whether ACKs are received successfully or not respectively.

	Figure \ref{fig_prob_2} represents the probabilistic structure of the problem as shown below.
	\[P(w_m^j=1|y_{m}^j=1)=k; P(w_m^j=0|y_{m}^j=1)=1-k\]
	\[P(w_m^j=1|y_{m}^j=0)=g; P(w_m^j=0|y_{m}^j=0)=1-g\]
	However, note that the ACK piggyback in current time instant will show quality of the channel for probe packet sent in previous time instant. Hence a transition probability of $g$ is introduced above.
	The probabilistic structure assumed for the system state transition is given below. We have assumed that $q>s$ and $k>g$ which signify respectively that the probability of a link becoming good from previous good state is higher than that of previous bad state and probability of successful ACK from good state is higher than that of bad state.
	\[P(y_{m+1}^j=1|y_{m}^j=1)=q; 	P(y_{m+1}^j=0|y_{m}^j=1)=1-q\]
	\[P(y_{m+1}^j=1|y_{m}^j=0)=s; 	P(y_{m+1}^j=0|y_{m}^j=0)=1-s\]
	\begin{figure}[h!]
		\centering
		\includegraphics[width=0.25\textwidth]{prob_structure.eps}
		\caption{Probabilistic structure of the exploration problem at a UE locally.}
		\label{fig_prob_2}
	\end{figure}
	For a given relaying zone $j$, let us define the information vector available locally to the zone $i$ till time instant $m$ as $H_m^j=(w_0^j,w_1^j,\cdots,w_m^j)$. The sufficient statistics or belief \cite{dp_book} (chapter 5) for this problem is defined locally for relaying zone $j$ as:
	\begin{equation}
	r_m^j=P(y_m^j=1|H_m^j)
	\end{equation}
	This equation gives the probability that a relaying link/zone is in good state given the previous history of information. The estimator function for the local system can be defined as:
	\begin{equation}
	r_{m+1}^j=\Phi'(r_m^j,w_{m+1}^j).
	\end{equation}
	Using Baye's rule we get,
	\begin{equation} \label{update_rule_2}
	r_{m+1}^j = \begin{cases}
	\frac{(qr_m^j+(1-r_m^j)s)k}{H}, & \text{if } w_{m+1}^j=1\\
	\frac{(qr_m^j+(1-r_m^j)s)(1-k)}{J} , & \text{if }w_{m+1}^j=0
	\end{cases}
	\end{equation}
	where,
	 \[H=(qr_m^j+(1-r_m^j)s)k+(r_m^j(1-q)+(1-r_m^j)(1-s))g\]
	  \[J=(qr_m^j+(1-r_m^j)s)(1-k)+(r_m^j(1-q)+(1-r_m^j)(1-s))(1-g)\]

	\textbf{Cost structure}: The cost for stopping the exploration for both cases when link cannot be and can be formed is $0$. However, if the link which seems to be good and selected, may become bad in upcoming time instants and thus causes packet loss.  So it will incur $D_1$ cost to compensate for the packet loss. Similarly, the link which was bad and not selected could have been a good relay link in upcoming time instants. Then it will incur some cost $D_2$ in order to compensate for further exploration.  If decision of relay selection cannot be made in current time instant, then probe packets are sent for another $\epsilon$ time unit to continue the exploration.  The cost incurred here is $c_\epsilon$. This will go on till time $M$ which is the upper bound on the exploration time for a single link. The objective is to derive a decision criterion, whether to select or not select the given relay link, or if this decision cannot be made then continue learning the current relay link quality. The expected cost is formulated as a dynamic program. At the end of the $(M-1)th$ period, the expected cost is:
	\begin{equation}
	K_{M-1}^j(r)= \min\{rD_1,(1-r)D_2\} \label{eq_j11}
	\end{equation}
	where $r$ is a variable denoting belief. For the time instant $m=M-2$, we have,
	\begin{equation}
	K_{M-2}^j(r)=\min\{ rD_1,(1-r)D_2,c_\epsilon+\mathbb{E}[K_{M-1}^j(r)] \label{eq_j12}
	\end{equation}
	where $\mathbb{E}[\cdot]$ is the mathematical expectation over the observations. Here the first term in minimization expression denotes the expected cost ($rD_1+(1-r)0$) incurred when the learning of the current link being explored  is stopped since it was in the bad state. The second term denotes the expected cost for selecting the current relay link being explored. The third cost is the expected cost for continuing exploring the current relay link for another round of $\epsilon$ unit of time. We can write a general expression as:

	\scriptsize
	\begin{equation}
	K_{m}^j(r)=\min\{rD_1,(1-r) D_2, c_\epsilon+\mathbb{E}[K_{m+1}^j(\Phi'(r,w))]\} \label{eq_j13}
	\end{equation}\normalsize
	where $w\in \{0,1\}$ is a variable denoting ACKs failure/success. Note that the notation $w_m^j$ defined earlier denotes the ACK failure/success of link $j$ at time instant $m$, whereas $w$ is used to denote a general variable for ACK failure/success. We will derive the optimal policy for this problem in section \ref{derivation_policy_2}.

	\section{Derivation of the Optimal Policy in  Exploration Phase} \label{derivation_policy_2}

	\subsection{Properties of $K_m^j(r)$}
	The general expression for time instant $m$ as mentioned in equation \eqref{eq_j13}, can be written equivalently as:
	\begin{equation}
	K_{m}^j(r)=\min\{rD_1,(1-r) D_2,E_m^j(r)\} \label{eq_dp_2}
	\end{equation}
	where,
	\begin{multline} \label{eq_j14}
	E_m^j(r)=c_\epsilon+P(w=1|r) K_{m+1}^j(\Phi'(r,1)) +\\ P(w=0|r)  K_{m+1}^j(\Phi'(r,0))
	\end{multline}
	For notation simplicity we will now remove the superscript $j$ from each of the respective notations, e.g., we will write $E_m^j()$ as $E_m()$. Hence $E_m^j(r)$ can now be denoted as $E_m(r)$.
	Now equation \eqref{eq_j14} can be further reduced as:

	\scriptsize
	\begin{multline} \label{eq_j15}
	E_m(r)=c_\epsilon +\\ \big((qr+(1-r)s)k+(r(1-q)+(1-r)(1-s))g\big)K_{m+1}(\Phi'(r,1)) +\\ \big(1-\big((qr+(1-r)s)k+(r(1-q)+(1-r)(1-s))g\big)\big)K_{m+1}(\Phi'(r,0))
	\end{multline}
	\normalsize
	The expected cost for the base case at the end of  $(M-1)th$ period is mentioned in equation \eqref{eq_j11}.
	We now show some of the properties of $E_m(r)$.
	\begin{prop} \label{prop_3}
		$E_m(r)$ is piece-wise linear and concave in $r$.
	\end{prop}
	\begin{proof}
		Proof given in Appendix \ref{A6}.
	\end{proof}

	\begin{prop}\label{prop_4}
		$\forall r \in [0,1]$, 	\[E_{m-1}(r) \le E_m(r) \le E_{m+1}(r)\]
	\end{prop}

	\begin{proof}
		Proof given in Appendix \ref{A7}.
	\end{proof}

	\subsection{Policy Structure}
	The structure of an optimal policy for our POMDP problem is provided in the following theorem.
	\begin{thm}\label{thm_2}
		The optimal policy for exploration in POMDP problem is a threshold policy.  At any time instant $m \in \{0,1,\cdots,M-1\}$, the optimal action is to stop exploration on current relay link and start exploring other relay links from $\mathbb{U}^i$ if $r_{m}\le \alpha_m$, or stop exploration on current relay link and choose it for data transmission if $r_{m}\ge \beta_m$, otherwise continue sending probe packet once more on the current relay link which is being explored to check the link quality if $\alpha_m < r_{m} < \beta_m$. Here $\alpha_m$, $\beta_m$ are appropriate constants.
		Also, we can say for the thresholds: $\cdots \ge \alpha_m \ge \alpha_{m-1} \ge  \cdots \ge \alpha_{1} $ and similarly  $\cdots \le \beta_m \le \beta_{m-1} \le  \cdots \le \beta_1$.

	\end{thm}
	\begin{proof}
		Proof given in Appendix \ref{A8}.
	\end{proof}

	 Using above result, we can say as $m \rightarrow \infty$, $\alpha_m$ and $\beta_m$ converges to some scalar $\overline{\alpha}$ and $\overline{\beta}$ respectively. This is because $\beta_m$ is bounded below and $\alpha_m$ is bounded above. Hence, for very large horizon length $M$, the optimal policy can be approximated by a stationary threshold policy with a time-invariant threshold $\overline{\alpha}$ and $\overline{\beta}$. In this case we can say that at time instant $m$, if  $r_{m}\le \overline{\alpha}$  then stop exploration on current link and start exploring other links from $\mathbb{U}^i$, if $r_{m}\ge \overline{\beta}$ then stop exploration on current link and choose it for data transmission, otherwise  if $\overline{\alpha} < r_{m} < \overline{\beta}$ then continue sending probe packet once more on the current link which is being explored to check the link quality.

	It is easy to check that $\Phi(r,w)$ is a
	non-decreasing function in $r$. Let us denote $r_0$ as the prior belief. We can say that, when prior belief satisfies $r_0>\Phi(r_0,w)=r_1$, then  $r_2=\Phi(r_1,0) < \Phi(r_0,0)=r_1$ and $r_2=\Phi(r_1,1) < \Phi(r_0,1)=r_1$. Proceeding in this way, we can show that $r_m$ strictly decreases with $m$ whenever we observe several successive ACK failures/successes. Similarly, when the prior belief satisfies $r_0<\Phi(r_0,w)=r_1$, then we can say that $r_2=\Phi(r_1,0) > \Phi(r_0,0)=r_1$ and $r_2=\Phi(r_1,1) > \Phi(r_0,1)=r_1$. Proceeding in this way, we can show that $r_m$ strictly increases with $m$ whenever we observe several successive ACK failures/successes. For getting $x$ successive ACK failures, we can define recursively a probability $\pi_x$ as: $\pi_1=\Phi(r_0,\overline{A})$,  $\pi_2=\Phi(\pi_1,\overline{A})$, $\cdots$, $\pi_x=\Phi(\pi_{x-1},\overline{A})$ with $\pi_0=r_0$.
	Similarly, for getting $x$ successive ACK successes, we can define recursively a probability $\pi_x'$ as: $\pi_1'=\Phi(r_0,A)$,  $\pi_2'=\Phi(\pi_1',A)$, $\cdots$, $\pi_x'=\Phi(\pi_{x-1}',A)$ with $\pi'_0=r_0$. Let $c$ and $d$ be the smallest integer such that $\pi_c \le \overline{\alpha}$ and $\pi_{d}' \ge \overline{\beta}$ respectively.  We can further simplify the stationary threshold policy as follows.
	\begin{cor}
		Using theorem \ref{thm_2}, we can simplify the optimal policy further as follows. Let $c$ and $d$ be the smallest integer such that $\pi_c \le \overline{\alpha}$ and $\pi_d' \ge \overline{\beta}$ respectively. If there are $c$ successive ACK failures, stop exploring on current relay link being explored and start exploring some other relay link. If there are $d$ successive ACK successes, stop exploring on current relay link being explored and choose it for data transmission. Otherwise, decision cannot be made and more probe packets are needed to be sent to learn the link quality.
	\end{cor}

	\section{Experiments and Results} \label{experiment}
	\begin{figure*}[h!]
		\centering
		\begin{minipage}[b]{.22\textwidth}
			\includegraphics[width=1\textwidth]{1.eps}
			\caption{$M$ vs \% of times no decision can be made}
			\label{fig-1}
		\end{minipage}\qquad
	\begin{minipage}[b]{.21\textwidth}
		\includegraphics[width=1\textwidth]{1_1.eps}
		\caption{$M$ vs \% of times no decision can be made when $M\ge3$}
		\label{fig-1_2}
	\end{minipage}\qquad
		\begin{minipage}[b]{.21\textwidth}
			\includegraphics[width=1\textwidth]{2.eps}
			\caption{No. of times relay are explored and switched with dynamic obstacles.}
			\label{fig-2}
		\end{minipage}\qquad
		\begin{minipage}[b]{.24\textwidth}
			\includegraphics[width=1\textwidth]{4.eps}
			\caption{Trade-off for average exploration time and average E2E delay.}
			\label{fig-4}
		\end{minipage}
	\end{figure*}

	Service region of dimension $100~m \times 100~m$ is divided into grid zones each of dimension $10~m\times10~m$. We have taken $\delta=100~ms$ and $\epsilon=1~ms$.
	UEs are using directional transmitter and receiver antennas for $60~GHz$ frequency with $G_r=G_t=6~dB$. LOS path loss exponent is $2.5$ and zero mean  log-normal shadowing random variable with standard deviation $3.5$ \cite{7974772,7109864}. Thermal noise density is $-174~dBm/Hz$ and devices are using $24~dBm$ transmit power.  Capacity of each link $(i,j)$ is $B\log_2(1+S_{ij})~bits/sec$, where $B=20~MHz$ \cite{al2014path} is bandwidth and $S_{ij}$ is the received signal to noise ratio. The packet length is fixed and of size $65535~bytes$. A source UE in a grid zone $i$ can transmit data to a maximum of $24$ grids surrounding it ($|\mathbb{U}^i| \le 24$).  There are maximum  of $16$ static and $D$ dynamic obstacles present in these $24$ grids surrounding the source UE, where $D\in\{0,4,8,12,16\}$. Each static obstacle is placed uniformly in the service region. Each dynamic obstacles is moving randomly and independently of each other and following a simple blockage model such that with probability $0.5$ it will block a given link otherwise it will not block the link.  A single source-destination pair is assumed for simplicity and all other devices in a given zone may act as relay. We have written our own C++ custom code and run them on a \texttt{GNU} $4.8$ compiler on Intel core $i7$ machine and average of 1000 runs are taken. The average exploration time and E2E delay for sending $100$ packets are the main parameters considered. We have compared our results with the received signal strength (\texttt{RSS Based}) approach which is most commonly used for relay selection during exploration after transmitter and receiver beams are perfectly aligned.

	Value of $M$ is derived experimentally in figure \ref{fig-1}. This value should not be too small, otherwise there won't be sufficient exploration to lead to a decision. Also $M$ should not be too big, otherwise a lot of time would be wasted in exploration. Figure \ref{fig-1} shows the percentage of time the decision could not be made on varying  values of $M$ for different $D$. It is clear that for almost all values of $D$, decision cannot be made when $M$ is $1$ and $2$. For $D=16$ and $M=3$, for around $3.5\%$ of cases decision cannot be made a shown in figure \ref{fig-1_2}. This figure shows that with dynamic obstacles, no. of times the decision of relay selection cannot be made, also increases. We have chosen $M=4$ for performing further simulations.

	Figure \ref{fig-2} shows the effect of dynamic obstacles on percentage of times new relays are explored and switched. With dynamic obstacles, the probability of link breakage increases and thus number of times new relays are explored and switched also increases. Our approach outperforms the RSS based approach because it  chooses the link by learning dynamic obstacle's presence while RSS based approach greedily chooses a link based on best RSS values which might have higher chance of blockage.	Figure \ref{fig-4} shows following results: first with dynamic obstacles, average E2E delay is increased which follows from previous result of figure \ref{fig-2}. Since higher number of dynamic obstacles will have higher chance of breaking a link and thus causes increase in E2E delay. Similarly, following the argument of previous result of figure \ref{fig-2}, our approach outperforms RSS based approach. Secondly, with dynamic obstacles, average exploration time is increased because higher link breakage with dynamic obstacles would cause more exploration and switching of relays. Here RSS based approach performs better since it always takes one time unit of exploration, whereas our approach takes additional exploration time units to learn the quality of link considering dynamic obstacles. Thirdly, it can be seen from above two results that, average exploration time is higher in our approach at the cost of reducing the cost of average E2E delay significantly compared to the RSS based approach. This describes the trade-off of additional exploration time  over the average E2E delay.

	\section{Conclusion}  \label{conclusion}
	The problem of selecting a given relay during exploration is investigated taking into account presence of dynamic obstacles. We have modeled this problem as a finite horizon POMDP	framework at each UE. Using this model, an optimal threshold policy is derived for each UE which is then simplified to a  stationary policy. This policy governs the source UE to  take decisions based on successive	ACK failures or success on the current relay link under exploration. Through simulations, the trade-off between average exploration time and E2E delay is shown, where our approach captures the effects of dynamic obstacles significantly compared to the RSS based approach.