Risk-Aware Motion Planning for a Limbed Robot with Stochastic Gripping Forces Using Nonlinear Programming

We present a motion planning algorithm with probabilistic guarantees for limbed robots with stochastic gripping forces. Planners based on deterministic models with a worst-case uncertainty can be conservative and inflexible to consider the stochastic behavior of the contact, especially when a gripper is installed. Our proposed planner enables the robot to simultaneously plan its pose and contact force trajectories while considering the risk associated with the gripping forces. Our planner is formulated as a nonlinear programming problem with chance constraints, which allows the robot to generate a variety of motions based on different risk bounds. To model the gripping forces as random variables, we employ Gaussian Process regression. We validate our proposed motion planning algorithm on an 11.5 kg six-limbed robot for two-wall climbing. Our results show that our proposed planner generates various trajectories (e.g., avoiding low friction terrain under the low risk bound, choosing an unstable but faster gait under the high risk bound) by changing the probability of risk based on various specifications.


I. INTRODUCTION
P LANNING complex motions for limbed robots is challenging because planners need to design footsteps and a body trajectories while considering the robot kinematics and reaction forces. Motion planning for limbed robots has been studied by a number of researchers. Sampling-based planning, such as the Probabilistic-Roadmap (PRM), samples the environment and generates the motion while satisfying static equilibrium and kinematics for a robot [1], [2]. Optimizationbased planning, such as Mixed-Integer Convex Programming (MICP) and Nonlinear Programming (NLP), solves the solution given constraints using optimization algorithms such as gradient descent [3]- [9].  While many papers discuss motion planning for the robot, few studies have investigated how planning is affected by stochastic gripping forces. One of the open problems in motion planning of a limbed robot equipped with grippers is the stochastic nature of gripping [10], [11]. For example, the gripping forces caused by spine grippers depend on the stochastically distributed asperity strength (Fig. 2). Thus, risk results from the randomness of the gripping force. By considering risk in a probabilistic manner, the planner can design a variety of trajectories based on various specifications.
The stochastic planning problem can be categorized into two approaches: robust approaches [5]- [7], [12] and riskbounded approaches [7], [8], [13]- [15]. In robust approaches, the planners design trajectories that guarantee the feasibility of the motion given the uncertainty bounds. A soft constraintsbased robust planning was investigated in [6], where the planner allows the solution to be at the boundary of stability. Tas showed the planner to remain collision-free for the worstcase uncertainty for automated driving [12].
On the other hand, the risk-bounded approach designs trajectories that guarantee the feasibility of the motion given the probability density function (PDF): it prevents the probability of violating state constraints (violation probability) from being higher than a pre-specified probability. Prete formulated a chance-constrained optimization problem of a bipedal robot c 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
arXiv:2006.02656v2 [cs.RO] 6 Jun 2020 by approximating a joint chance constraints with linear inequality constraints [7]. Planning on slippery terrain was in [8], where the planner utilizes the prediction of the coefficient of friction to design the motion of the body and footsteps, respectively. Our approach is similar to [8], but we model the stochastic contact force of the robot and formulate the planning algorithm considering the trajectory of a body and footsteps simultaneously. For tasks with a higher probability of failure (e.g., climbing on slippery terrain) [15], the risk-bounded approach has advantages over the robust approach. Because the robust approach often uses a much less informative deterministic model, it is likely to generate conservative solutions with the worstcase uncertainty bound. For demanding tasks, this may be infeasible, with such a planner generating no possible solution and failing to achieve specified goals. In contrast, because the risk-bounded approach can be more aggressive, the problem may be feasible, generating trajectories that carry a probability of failure through risk-taking alongside a non-zero chance of successfully achieving the goal. The violation probability provides a tuning knob to define a Pareto boundary on the risk between failure while finding a trajectory vs. failure while executing a found trajectory. This user-defined parameter can be task-and environment-specific, in contrast to the rigidity of the robust approach.
In this paper, we address a motion planning algorithm formulated as NLP for a limbed robot with stochastic gripping forces. Our proposed planner solves for stable postures and forces simultaneously with guaranteed bounded risk. In addition, chance constraints are introduced into the planner that restrict contact forces in a probabilistic manner. We employ a Gaussian Process (GP), a non-parametric Bayesian regression tool, to acquire the PDF of the gripping force. Our proposed motion planning algorithm is validated on an 11.5 kg hexapod robot with spine grippers for multi-surface climbing. While we focus on multi-surface robotic climbing with spine grippers in this paper, our proposed planner can be applied to other robots with any type of grippers for performing any task (e.g., planning of walking, grasping) as long as the robot has contact points with stochastic models.
The contributions of this paper are as follows: 1) We formulate risk-bounded NLP-based planning that considers the stochasticity of gripping forces. 2) We employ the Gaussian Process to model gripping forces as random variables. 3) We validate the algorithm in hardware experiments.

II. PROBLEM FORMULATION
This section describes the friction cone considering maximum gripping forces, model of a position-controlled limbed robot with multi-contact surfaces, and the modeling process of a gripping force through GP.

A. Friction Cone with Stochastic Gripping Forces
With grippers, the friction cone constraint can be relaxed on the contact point. For our spine-based gripper, even under a zero normal load, the spines insert into the microscopic gaps on the surface (Fig. 2), generating a significant amount of shear force ( Fig. 6) [19]. For a magnet-based gripper, the reaction forces includes the additional magnetic force imposed by the gripper itself, offsetting the friction cone as seen by the rest of the robot. Thus, we modify the regular friction cone, adding in an offset shear force when a normal force is zero to account for the gripping force. As the normal force increases, the maximal allowable shear force increase in the same way as a regular frictional force, with a coefficient of friction λ that is assumed to be a constant only depending on the property of the contact surface. This contact model is illustrated in Fig. 3, where f r is the reaction force between the surface and the gripper. f g,m is the maximum gripping force from grippers under a zero normal force. Note that f g,m is measured per gripper as a unit. In general, f g,m can have both normal and shear components. However, for our spine grippers, the normal component of f g,m is relatively small, so we assume that the gripper generates only shear adhesion. The gripper does not slip when f r is within this friction cone, as indicated by the shaded region in Fig. 3. Since the interaction between the micro-spines and the surface is highly random, f g,m is naturally modeled as a Gaussian random variable. However, the orientation of the spine and the number of spines in contact with the surface also change as the orientation of the gripper changes, which leads to a shift of the mean and standard deviation of f g,m . We learn this model from data by GP. With GP, our proposed planner is able to deal with the stochastic nature of gripping taking into account the gripper orientation.

B. Model of Reaction Force Using Limb Compliance
During multi-surface locomotion, the robot leverages the compliance from its motors in order to squeeze itself between multi-surfaces, as depicted in Fig. 2. One difficulty multi-limbed robots have is that reaction forces are statically indeterminate [16]. Consequently, reaction forces cannot be determined by static equilibrium equations when the robot supports its weight more than three contact points. Hence, in order to calculate the reaction force under this condition, the deformation of the robotic system should be considered.
From the standard elasticity theory, f r can be described as the spring force using the Virtual Joint Method [17]: where K is the stiffness matrix for H degree-of-freedom limb. k is a diagonal matrix that has k i diagonal elements, and k i is the spring coefficients of the position-controlled servos. J is a 3 × H Jacobian matrix. The deflection is imposed by terrain where δ wall is the displacement of the terrain and δ CoM is the body deflection, sag-down, due to the compliance as shown in Fig. 2.

C. Model of Gripping Force Using Gaussian Process
The objective of using GP is to predict the maximum gripping force f g,m in a probabilistic way.
There are many design decisions that go into the formulation of the GP problem, including choice of kernel, distance metric, and associated weighting between state variables [18]. We can start with the simplest formulation with all state variables equally weighted under the Euclidean distance metric using the squared exponential kernel as a starting point. In practice, this choice was observed to work well enough to not necessitate further design. A more general characterization of the effects of these hyperparameters can be found in [18]. In this work, we assume that the maximum gripping forces by spine grippers is a function of the gripper orientation and the coefficient of friction [10], [11], [19], [20]. This is because with a microscopic view, the spine-asperity interaction is different depending on how a spine is inclined with respect to the asperity as shown in Fig. 2. GP can handle the effects of other unmodeled parameters by treating them into uncertainty. Hence, the state s is a four-dimensional vector with s = [α, β, γ, λ] where α, β, γ are the Euler angles along x, y, z axis defined in Fig. 4.
Here, we assume that the maximum shear force follows Gaussian distribution. Given a data set S = {s 1 , · · · , s n } with the measured shear forces y g,m = [y g,m 1 , . . . , y g,m n ] , the maximum shear force f g,m by a gripper can therefore be modeled as: is a positive definite kernel. In this work, we employ the squared exponential kernel as follows: where σ 2 f represents the amplitude parameter and l defines the smoothness of the function f g,m .
Here, let D = [s 1 , · · · , s n ] be the matrix of the inputs. In order to predict the mean and variance matrix at D * , we obtain the predictive mean and variance of the maximum shear force by assuming that it is jointly Gaussian as follows: , and σ 2 n is the variance of the Gaussian observation noise with zero mean.
Our GP procedure can be generalizable to model other gripping forces as long as the gripping force changes continuously as the orientation of the gripper changes. For instance, the GP approach can be used to model the gecko gripper force [21] using the detachment angle as the state of the GP.

D. Spine Gripper for Wall Climbing
A three-finger spine-based gripper was designed (Fig. 4) using spine cells based on [20]. Each finger consists of a spine cell with 25 machine needles loaded with 5 mN/mm springs, and a slider mechanism holds the cell with one compliant plastic spring. The diameter of the needle at the tip is 0.93 mm, and it is made of carbon steel. The gripper center module includes one spine palm with the same spine configurations as the cells. The attachment component is fixed at the tip of the robot limb at 37 • from the limb axis to maximize the contact area. The finger, center, and attachment members are assembled with a one-slider, two linkage mechanism (Fig. 4). This linkage system is designed to provide a passive micro grip as the center palm presses up against a wall. The three fingers are located at 120 • apart from each other in z-axis and tilted about 15 • from z-axis.

E. Data Collection
To collect a dataset, maximum gripping forces f g,m were evaluated with a minimal normal force at varied orientations as summarized in Table I. We collected 20 data sequences for every orientation as a training dataset and 10 data sequences as a testing dataset. The coefficient of friction between spines and environments was measured by loading a constant mass on the gripper. A small activation force is necessary to compress spine springs and ensure that the spines are touching the wall, but is assumed to be negligible. These orientations were selected to cover possible gripper angles during regular wall climbing. The gripper was fixed to a linear slider at an orientation and pulled by a force gauge on 36-grit and 80-grit The predicted maximum gripping force PDF from GP, the training data PDF, and the testing dataset sandpapers that are commonly used to emulate rough surface with microscopic asperities [20], as shown in Fig. 5. The GP hyperparameters were optimized using the BFGS algorithm [22]. The obtained testing data with the predicted PDF of the maximum gripping force and the PDF of the training data is illustrated in Fig. 6. The predicted maximum gripping force and the training data are displayed as a mean ± with a 95 % confidence interval. Overall, we show that the GP prediction works well with different states.

III. RISK-AWARE MOTION PLANNING
In this section, we present a complete risk-aware motion planning algorithm formulated as (8a)-(8k). The objective of our proposed planner is to find the optimal trajectory for the Center of Mass (CoM) position, its orientation, the foot position, and the reaction force for each foot in order to arrive at the destination while satisfying constraints. Our proposed planner enables the robot to find feasible trajectories that consider risk from the grippers under various environments.
We define one round of movement made by a robot when its body and all of its limbs have moved onto the next footholds. Note that for each round, the planner investigates several critical instants between two postures with pre-defined gait as explained in detail in Section IV. At j-th round, Γ is the decision variables that are given as: where p i,j is the foot i position, P CoM,j is the position of the body, Θ CoM,j is the orientation of the body, θ i,j are the joint angles for the limb i, and f r i,j is defined in Section II-A. In this study, f g,m i,j is treated as a random variable based on the model |∆Θ CoM | ≤ ∆Θ Th (angular stride) (8c) . Equation (8a) is the cost function that depends on the robot's state. Equation (8b), (8c), and (8d) bound the range of travel between rounds. Equation (8e) represents the forward kinematics constraints. In (8f), it ensures that p i,j is within the feasible terrain where the robot is able to put its limb. In this paper, we assume that the robot generates a quasistatic motion. Hence, the planner has the static equilibrium constraints expressed by (8g) and (8h), where F tot and M tot is the external force and moment, respectively. In this work, only gravity is considered as the external force. Equation (8i) and (8j) ensure that the motor torque is lower than the maximum motor torque where J (θ i ) is a Jacobian matrix. The reaction force f r i is constrained by (8k), which describes the friction cone constraints to prevent the robot from slipping where λ i,j (p i,j ) denotes the coefficient of friction at p i,j . Note that this constraints (8k) is also stochastic constraints due to f r i . Equation (8k) can be converted into deterministic constraints, which is explained in Section III-B.
Compared to sampling-based approaches such as RRT, NLP is able to formulate relatively complicated constraints such as friction cone constraints (8k), which are typically difficult for the sampling-based approaches to handle in terms of computation. In addition, MICP approaches such as [3], [5], [9] can increase the computation speed by decoupling the pose state from wrench states. However, they potentially limit the robot's mobility. The robot may not choose the trajectory on the low friction terrain in case the planner first solves the pose problem and then solves the wrench problem later since the pose optimization problem does not consider the wrench information. Although MICP can plan the trajectories considering both wrench and pose state simultaneously, it needs to sacrifice the accuracy by assuming an envelope approximation on bilinear terms [3] or allow relatively expensive computation as the number of the integer variables increases, which is intractable for high degree-of-freedom (DoF) robots (e.g., our robot has 24 DoF). In contrast, NLP can simultaneously solve the trajectory reasoning both the pose and the wrench with relatively less computation [4].
1) Kinematics: Forward kinematics (8e) is given as: where R(Θ CoM ) is the rotation matrix from the world frame to the body frame, p b i is the foot position relative to the body frame.
2) Feasible Contact Regions: We utilize NLP to formulate the planning algorithm so that any nonlinear terrain (i.e., nonflat terrain), such as tube and curve, can be directly described.
If a robot traverses on the flat terrain, the footstep regions are convex polygons as follows:

B. Chance Constraints
Here, we show that the friction cone constraints in (8k) can be expressed using chance constraints, which allow the planner to convert the stochastic constraints into deterministic constraints.
One key characteristic of robotic climbing is that climbing is a highly risky operation: a robot can easily fall without planning its motion correctly. Hence, it needs to restrict reaction forces using the friction cone constraints given as: To decrease the computation of solving for the NLP solver, we simplify the (12) by linearizing them as follows: where ζ i,j , ξ i,j are any tangential direction vectors on the wall plane.
Regarding (8k) formulated as (11), (13), (14), we rearrange the equations and the joint chance constraint is given by: where α k i,j are coefficient vectors, and β k i,j are coefficient scalars that consist of the convex polytopes defined in (11), (13), (14). In (15), M denotes the number of constraints defining the polytopes. ∆ is the user-defined violation probability, where the probability of violating constraints is under the ∆. We can regard ∆ as relating to the likelihood that gripper slip will be responsible for the failure of the robot. For example, if ∆ is high, the planner can explore a larger space because the feasible region expands in optimization. As a result, the robot tends to plan a trajectory with a high violation probability by assuming that the gripper generates enough force. For a robotic climbing task, these chance constraints enable the robot to perform challenging motions that would be infeasible without considering the gripping force. In contrast, if ∆ is small, the planner tends to generate more conservative motions because the robot assumes that the gripper does not output enough force to support the weight of the robot.
Imposing (15) is computationally intractable. Thus, using Boole's inequality, Blackmore [13], showed that the feasible solution to (15) is the feasible solution to the following equations: for all j = 1, . . . , N, k = 1, . . . , M . The violation probability for each constraint per round ∆ j,k is constrained in (17), in order not to exceed the given ∆. Because non-uniform risk allocation (17) is also computationally expensive [14], we use the following relation: Thus, the stochastic constraints (16) can be then converted into a deterministic constraint as given by: where Φ is the cumulative distribution function of the standard normal distribution. It can be transformed further as follows: where Φ −1 is the inverse function of Φ.

C. Cost Function
The cost function consists of intermediate costs and a terminal cost. In this work, the target mission is to arrive at the destination. Thus, the terminal cost is the distance from the position of the last pose to the destination.
where W D is the weighting matrix and q N = p 1,N , . . . , p L,N while q d is the configuration at the destination. The intermediate costs restrict a large amount of shifting in terms of linear and rotational motion of a body and the foot position as follows: where W BPos , W Foot , and W BRot are the weighting matrix.

D. Two Step Optimization for a Position-Controlled Robot
Although our proposed motion planner works for any limbed robot, there is a drawback for a position-controlled robot when wall-climbing. For the position-controlled robot, it is necessary to compute how much δ wall is necessary to generate the planned reaction forces. Therefore, the planner needs to include additional constraints from (1), (2) to realize the planned trajectory. However, we observed that the nonlinear solver has a numerical issue with (2), so it is intractable for the solver to solve our proposed NLP in (8a)-(8k) with (1), (2). To avoid this problem, we decouple the optimization problem into two-step problems shown in (23a)-(23l) and (24a)-(24d): the first planner is in charge of the pose and the reaction force of the robot, and the second planner finds δ wall,i,j , δ CoM,i,j , which are the control inputs to a position-controlled robot. In (24a)-(24d), the constraint (24b) ensures that δ wall,i,j is bounded under a certain threshold. We argue that this decoupling is reasonable because the first planner solves the "essential" problem (e.g., How much reaction force is necessary? What is the footstep trajectory?) to plan the force and pose trajectory. The second planner only computes the control input to the motors, and it does not have a significant effect on the entire motion planning. As explained, if the robot is force controlled, the planner does not need to consider (1), (2). As a result, the second optimization is not necessary for a force-controlled robot, and the whole motion is planned only based on the first optimization problem.

IV. RESULTS
In this section, we evaluate our proposed planner by testing the robot's performance in three different tasks: energyefficient climbing, climbing on non-uniform terrains, and climbing with a tripod gait.
We utilize Ipopt solver [23] to solve the planning problem on an Intel Core i7-8750H machine. The derivative of constraints are provided by CasAdi [24]. The optimizer is initialized with the default configuration of the robot (Fig. 1, bottom  configuration), and the specifications of the computation for Section IV-B is summarized in Table II.
We implement the results of our proposed planning algorithm (i.e., the motion plan), on a six-limbed robot, each limb of which has three DoF. Each joint uses pairs of Dynamixel MX-106 motors, providing a maximum torque at 27 Nm. The robot is equipped with a battery, computer, and IMU. The robot runs a PID loop to regulate its body orientation. No other sensor is used to control its linear position. The robot weighs 11.5 kg. The width of the robot's body is 442 mm while its height at its standing state is 180 mm. In each experiment, the robot climbs between two walls at a distance of 1200 mm, where the wall is covered with sandpapers of different grit size to adjust the coefficient of friction. All hardware demonstrations can be viewed in the accompanied video 1 .

A. Energy Efficient Planning
The objective of this task is to assess the consumed energy of climbing with two different violation probabilities. While the robot can grip the wall with a low violation probability (e.g., ∆ = 0.0005), there is a disadvantage of consuming more energy. On the other hand, the robot may perform an energyefficient motion with a higher violation probability (e.g., ∆ = 0.1). Here, we set N = 7, M = 6 to compute ∆ j,k . To show the trade-off between the consumed energy and the violation probability, we let the robot climb on the walls with one leg gait where the robot first lifts its right front limb, puts it on the next position, pushes its body up, lifts its right middle limb, and so on. Within each round, the planner investigates 12 critical instants for one leg gait: 6 instants after the robot   lifts one limb, and 6 instants after the robot places the limb on the next position and pushes its body. The planner solves the optimization problem for these 12 instants. We measure the current I i,t and the voltage V i,t of each limb i online when the robot climbs on the wall covered by the 36-grit sandpapers with one leg gait and estimated the power per one limb at every sampling time t. The power P i,t is estimated as follows: We plot the consumed power for two consecutive limbs from the hardware experiment in Fig. 7. Fig. 7 shows that the consumed power of a limb decreases when the limb is in the air while the other limbs increase the consumed power to increase the reaction force. Furthermore, the robot consumes more power with smaller ∆, which means that the robot needs to push the wall to increase f r . In contrast, if ∆ = 0.1, the solution requires less power, but has a larger probability of slipping. In Fig. 8, the total consumed energy from these limbs was calculated by integrating their power over time spent climbing. In our robot, the robot could decrease the energy by 46.5 % under ∆ = 0.1 compared with the energy under ∆ = 0.0005. In the left panel, the robot puts its feet on low and high friction terrain by taking a high risk bound. In the right panel, the robot puts its feet only on high friction terrain.

B. Climbing on Non-uniform Walls
This scenario demonstrates that the robot designs different trajectories under the different violation probability to climb on walls with varying coefficients of friction. The planned trajectories are shown in Fig. 9. In this example, the robot climbs between two walls where the terrain shown in black is covered by 36-grit sandpapers (λ = 2.3), the terrain shown in green is covered by 80-grit sandpapers (λ = 1.1), and the terrain shown in red is covered by the material with λ = 0 as shown in Fig. 9. The varying coefficients of friction are modeled by a parabola function, which encourages the solver to converge on a solution. In the left panel of Fig. 9, the violation probability ∆ is 0.1 while in the right panel, the violation probability ∆ is 0.001 for M = 6 and N = 7.
The left panel of Fig. 9 illustrates that the robot avoids the red area (zero friction) and puts its foot mostly in the black area (high friction), but sometimes also in the green area (low friction) to minimize the trajectory length. In the right panel of Fig. 9, the violation probability is decreased, and the robot footsteps completely remain inside the high friction area. As a result, our proposed NLP-based planner operates the pose and forces together and makes a trade-off between a shorter but more risky trajectory and a longer but safer trajectory. This cannot be achieved if the planner decouples the footstep and force planning, such as in [9]. Fig. 10 shows the trajectory with higher risk bound ∆ = 0.1 and compares the foot location at t = 146 with the foot location with ∆ = 0.001 in the hardware experiment. We notice that at t = 146 s, the foot touches the white area where the coefficient of friction is 0, which never happened with ∆ = 0.001. Since the robot only controls its body orientation based on IMU feedback and does not control its linear position, the implemented trajectory does not strictly follow the planned one. We observe that lower risk bound is beneficial in this situation to avoid failure since it compensates for the tracking error by the imperfect controller.

C. Climbing with Less Stable Gait: Tripod Gait
In this scenario, we demonstrate that the robot can conduct a tripod gait, when it lifts three legs simultaneously, by setting the violation probability much higher. Before installing the gripper on the current six-limbed robot, it was almost infeasible to climb on the walls with the tripod gait because of the torque limits of the motors. With the grippers installed,  however, the robot has a greater chance to climb on the walls with a tripod gait. If we set ∆ = 0, the problem is infeasible since the constraints under the worst-case uncertainty are conservative. This result would be equivalent to the results of other robust algorithm such as [12], where the optimizationbased robust approach with the worst-case uncertainty is proposed. However, by utilizing the chance constraints and increasing the violation probability, the planner generates a feasible solution. In our trial, we set the violation probability ∆ = 0.4 for M = 6 and N = 3, and allowed the robot to climb on a wall covered by 36-grit sandpapers. The planner investigates 4 critical instants: 2 instants after the robot lifts three limbs, and 2 instants after the robot places them down and pushes its body. The planned trajectory is illustrated in the left panel of Fig. 11. As shown in the right under the condition, the robot succeeded in climbing on the walls with the tripod gait and its climbing velocity was 2.5 cm/s, which is three times faster than the one leg gait.

V. CONCLUSION AND FUTURE WORKS
In this paper, we presented a motion planning algorithm for limbed robots with stochastic gripping forces. Our proposed planner exploits NLP to simultaneously plan a pose and force with guaranteed bounded risk. Maximum gripping forces are modeled as a Gaussian distribution by employing the GP, which provides the planner with the mean and the covariance information to formulate the chance constraints. We showed that under our planning framework, the robot demonstrates rich -sometimes drastically different -behaviors, including planning a risky but energy-efficient motion versus a safe but exhausting motion, avoiding danger zones like low friction environments, and choosing fast but less stable motions (i.e., a tripod gait) based on the different violation probabilities ∆ in hardware experiments.
The current limitation in this work is that the actual probability of failure is not strictly equal to pre-defined ∆ because other sources of uncertainty exist, such as sensor noises. In future work, we will extend our planner to take into consideration these sources.