Walid Amehri

Gang Zheng

Alexandre Kruszewski

FEM-based Exterior Workspace Boundary Estimation for Soft Robots via Optimization

Keywords: Soft robotics, Workspace boundary, Constrained optimization

This paper investigates the exterior workspace boundary of a soft robot with a certain configuration controlled by equipped bounded actuators. To achieve this, we implement an optimization-based approach on the studied soft robot which has been modeled by the Finite Element Method (FEM). Finally, we provide numerical simulations of various configurations to demonstrate the validity of the suggested technique, which, in comparison to the conventional forward method, may considerably minimize the complexity of exterior workspace boundary estimation.

I. INTRODUCTION

Soft robots are made of deformable materials that may be widely found in the nature, such as rubber, origami, silicone, and several deformable materials [START_REF] Iida | Soft robotics: challenges and perspectives[END_REF]. With their natural conformity and low stiffness, they are the safest device for human interaction [START_REF] Yim | Design and rolling locomotion of a magnetically actuated soft capsule endoscope[END_REF] and the ideal tool for dynamical environments operations [START_REF] Shepherd | Multigait soft robot[END_REF], resulting in their increasing use in many applications.

However, the process of designing soft robots [START_REF] Rus | Design, fabrication and control of soft robots[END_REF] is still governed by trial-and-error or bio-inspired notions. Given a particular soft robot's configuration, evaluating its reachable workspace is still an open subject, and is very essential for other soft robotics' main scientific challenges, such as control [START_REF] Della Santina | Dynamic control of soft robots interacting with the environment[END_REF], trajectory planning [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], and optimization design [START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF].

Concerning the modeling techniques, different methods have indeed been reported in the literature for soft robots' Walid Amehri, Gang Zheng and Alexandre Kruszewski are with University of Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France flexibility characterization. On the one side, the soft robotics extensively uses the modeling technique PCC (Piece-wise Constant Curvature) [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF], but again the constant curvature hypothesis is not always accurate, notably when external stresses are applied on the robot. On the other side, the PCS (Piece-wise Constant Strain) model [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF] can be regarded as an extension of the PCC method. This method accounts for both geometric non-linearity and material proprieties, in which a finite group of strain vectors is used to characterize the deformation of the investigated soft robot. The PCS technique, like the PCC method, can only be used to model continuously slender soft robot. In contrast, an appropriate method for modeling the general geometric shape of soft robots is the Finite Element Method (FEM) [START_REF] Zheng | Controllability pre-verification of silicone soft robots based on finiteelement method[END_REF], and may, though at a high dimension model, give a precise soft robots model.

In the literature, various approaches were explored for the workspace estimation of rigid-link mechanisms. In [START_REF] Haug | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF], a continuation technique based on bifurcation theory was provided to estimate the workspace boundaries of several rigid-link robots. Next, a method based on analytical study of Jacobian was established in [START_REF] Abdel-Malek | Analytical boundary of the workspace for general3-dof mechanisms[END_REF] to determine the investigated mechanism unique behavior. Following that, in [START_REF] Gouttefarde | Interval-analysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF], cablesdriven parallel rigid manipulator's wrench-feasible workspace was assessed using interval analysis techniques. Then, an optimization technique to determine the limit of rigid mechanism workspace was suggested in [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF]. Despite their excellence in defining the workspace of various types of rigid robots, the strategies noted above cannot be explicitly extended to evaluate the workspace of soft robots since the direct kinematic model explicitly characterizes the input/output (actuation/end-effector position) relation for the rigid case, but such a connection is typically represented in an implicit way by a complex static model for soft robots when modeling the deformation of the soft material, which will be discussed in Section III.

In the soft robotics community, several approaches have been proposed to estimate the workspace of soft robots. A Monte-Carlo based approach was presented in [START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF] to characterize the workspace of concentric tube continuum robots. In [START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF], an approximation method was used to determine the workspace boundaries of continuum robots by sweeping the robot from the neutral configuration to the fully curved configuration, for the purpose of deducing the workspace surface from the robot's tip trajectory. The above two workspace estimation approaches are exhaustive since they consist of exploring the inputs domain (actuation domain) which dimension is not stable and may vary corresponding to the investigated soft robot's configuration. Next, an interval analysis based approach using FEM was proposed to estimate the workspace of soft robots [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF]. This workspace estimation approach, and even if it consists of discretizing the outputs domain (end-effector domain) which have a stable dimension regardless of the inputs dimension, is still exhaustive in a way that it consists of exploring the whole possible configurations of the workspace instead of mapping its boundary. On the other hand, the work of [START_REF] Amehri | Discrete cosserat method for soft manipulators workspace estimation: an optimizationbased approach[END_REF] provides an optimization-based approach that was implemented on the Discrete Cosserat model to estimate the workspace of soft slender manipulators by mapping their exterior workspace boundary. Being inspired by the above work, this paper extends the applicability of the optimization approach by implementing it on the adopted static FEM model [19, p. 428] to estimate the exterior workspace boundary of soft robots. Such an extension is motivated by the following pros. Firstly, FEM can be used to model any type of soft robot, while Cosserat method is limited to model slender type one. Secondly FEM can provide more accurate modeling precision than Cosserat method when modeling soft robots (even for slender type one), since Cosserat theory assumes that all cross-sections (which are perpendicular to the centroidal line of slender soft robots) are rigid and undeformable. Such an assumption is imposed for the sake of modeling simplicity, but inevitably decrease the modeling precision. However, to ensure the above pros (i.e., generality and precision), the price to be paid when using FEM is firstly the resulted high-dimensional model. Therefore there is a tradeoff between dimension and precision when fixing the number of elements used to model soft robots. Another price to be paid is the computation complexity: when solving the proposed optimization problem, it is necessary to calculate the gradient of the cost function. Analytical gradient can be deduced when using Cosserat, while it is not possible to have analytical one when using FEM to model soft robots. Consequently, numerical gradient has been used for FEM which increases again the computation time, compared to Cosserat.

II. PROBLEM STATEMENT

We use the structure of a tendon-driven trunk-like soft robot, as illustrated in Fig. 1, to explain clearly the examined topic. Provided the intended soft robot's specific shape, and attached actuators (i.e., tendons) whose amplitudes (i.e., tendons' tension) are in reality limited (noted as u P U), the point that should be obviously investigated is whether the soft robot's end-effector can achieve the yellow position in Fig. 1 given the existence of a valid input control u P U? Denote q as the nodal displacement vector, the end-effector position is noted as p E P R 3 , the end-effector displacement vector with respect to its position p E is noted as q E P R 3 . Then the definition of the workspace can be established by: Definition 1. Given a soft robot's modeled by : q " Φpu, q, 9 qq, with a bounded actuator u P U, the end-effector's p E " gpqq workspace W E is a subspace of R 3 , defined below: W E " tp E |Du P U, q, s.t.Φpu, q, 0q " 0, p E " gpqqu Consistent with Definition 1, the studied problem can be investigated via the ascertainment of the existence of a feasible input u P U allowing the containment of the yellow point (depicted in Fig. 1) inside the workspace W E . Hereafter, the examined topic will be tackled by firstly introducing the studied soft robot's FEM static model, and secondly by implementing two different model-based workspace estimation methods on the designed soft robot.

III. STATIC MODEL OF SOFT ROBOTS

The static model of soft robots is determined using the Finite Element Method in this study. We discretize the flexible field into tiny finite elements n e to obtain the structure's mesh which contains n p nodal displacements (equivalently DoFs). Diverse topologies presented in [START_REF] Reddy | An introduction to the finite element method[END_REF] can be used to obtain the mesh of the structure. We adopted the linear tetrahedron element for two primary reasons in this study. On the one hand, using linear tetrahedron elements ensures that the proposed workspace estimation approach can be applied to a variety of soft robot setups and shapes. On the other hand, efficient integration of element matrices can be ensured using linear tetrahedron elements. The selected mesh and the associated discretized representation of the domain are depicted in Fig. 2 for a soft robot in the concept of a trunk. The workspace, in agreement with Definition 1, comprises all equilibrium points. We may derive the following static model by applying this property to the robot's motion equation [19, p. 428], which is described by the following equation:

K N pqq " H T pqqu `F pqq (1)
where the nodal displacement vector is represented by q P R 3np , the internal force is noted as K N pqq P R 3np , the matrix to describe the actuator direction is noted as H T pqq P R 3npˆnu , the control input vector is noted as u P R nu , the dimension of the inputs is noted as n u , and the external force vector acting on the soft robot is presented by F pqq P R 3np .

Assumption 1. According to [19, p. 430], the internal forces of the structure K N pqq can be approximated as:

K N pqq " Kpqq q (2)
where Kpqq can be calculated by putting each element's stiffness matrix together as follows: Kpqq " ř ne e"1 K e pqq, the stiffness matrix of an element e depends on the elasticity matrix ∇, the strain matrix B e pqq, and the volume V e pqq of the element e. And it can be computed via the following formula [19, p. 429]: K e pqq " ş Vepqq B T e pqq ∇B e pqqdvpqq.

IV. EXTERIOR WORKSPACE BOUNDARY ESTIMATION

With the workspace being defined in Definition 1 and imposing Assumption 1 on the static model (1), the following is used to govern soft robot's workspace:

Kpqq q " H T pqqu `F pqq (3a) q E " Cq (3b) p E " q E `pp0q E (3c)
with the end-effector position being noted as p E , the endeffector's displacement vector with respect to its initial position p p0q

E is denoted by q E , and the matrix of selection related to the node coordinates of the end-effector is symbolized as C.

For all admissible inputs u P U, estimating W E in accordance with (3) can be achieved by calculating the frontiers of p E . However, based on (3), the expression Kpqqq varies depending on the material characteristics and is nonlinear, and the input-output mapping φ : U Ñ W E is implicit and nonlinear. To this aim, numerical techniques will be applied to estimate the mapping φ, and two different manners (the forward method and the optimization-based method) will be studied next.

A. Forward method

Firstly, the mapping φ can be iteratively approximated using a classical technique (namely the forward method) which consists of the discretizing the input bound U and finding the corresponding set of end-effector poses which represents the workspace W E of the studied soft robot (as illustrated by Fig. 3). By discretizing U, we can obtain U d . Next, given the information of q pj´1q and p pj´1q E which represent the nodal displacement and the position of the end-effector corresponding at the pj ´1qth iteration respectively (for j " 1, q p0q represents the initial nodal displacement and p p0q E represents the initial end-effector's position), then for any u P U d , we apply Newton-Raphson method to solve (3) and get q pjq by: q pjq " q pj´1q ´Ψq pu, q pj´1q q Ψpu, q pj´1q q, for any u P U d where Ψ q pu, qq denotes the gradient of Ψpu, qq with respect to the nodal displacement vector q, with Ψpu, q pj´1q q " Kpq pj´1q qq pj´1q ´HT pq pj´1q qu ´F pq pj´1q q Then, the position of the end-effector with respect to its nodal displacement q pjq is calculated as: p pjq E " Cq pjq `pp0q E . Finally, as illustrated in Fig. 3, the estimation of W E may be obtained by iterating the aforementioned approach for every possible input u P U d .

We would like to emphasize that the computation complexity of this approach is relatively high. In this approach, we have n nu d potential possibilities of input u since u P U is a subset of R nu and the bound for each entry of u is discretized into n d subsets. In this context, the forward method's computing complexity will increase as the dimension of inputs n u becomes larger, and will exponentially increase if the discretization precision n d increases, even if n u is low.

In the following, a novel technique by discretizing the endeffector domain will be proposed to reduce the computing cost when calculating the workspace of soft robots, rather than discretizing the actuators domain U with variable dimension.

B. Optimization-based method

The proposed method consists of mapping the surface enveloping soft robots' workspace using numerical algorithms: Step 1: Workspace boundary BW E definition

In accordance with Definition 1, for a specified positive radius ε, the neighborhood of p E can be defined as: B ε pp E q " tp E P R 3 | p E ´p E ă εu. Using the above statement, we can consistently define the workspace boundary BW E as: Definition 2. The workspace boundary BW E represents a subset of its workspace W E , and can be defined as follows:

BW E " tp E P W E | @ε ą 0, Dp E P B ε pp E q, s.t. pE R W E u
Step 2: A single point on BW E determination In this step, we explain how to find an initial boundary point p b E on the workspace boundary BW E (as depicted via a simplified illustration in Fig. 4). In fact, this step is necessary and needs to be treated before we proceed toward mapping the workspace enclosing surface. where the obtained actuation u b and its associated displacement vector q b represent the optimal and feasible solution (since they satisfy the bounded constraint and the nonlinear constraints) to minimize the distance between p E pqq and v. Finally, we deduce the related end-effector boundary position using the following: p b E " Cq b `pp0q E . The radiating point v should be outside the workspace, as indicated in Fig. 4. For this, we simply have to select v such that || Ov|| ą L max , where O " r0, 0, 0s T is the origin of the reference frame, since in practice the length L of a soft robot is bounded by L max .

Step 3: Workspace boundary BW E mapping

The final step consists of mapping the workspace boundary BW E by employing a numerical technique. To realize this, we can proceed by finding the boundary points corresponding to multiple successive directions (as described by a simplified illustration in Fig. 5). denoting the basic matrices of rotation with respect to the x, y and z axis respectively. Thus, the direction vector v pjq is equal: v pjq " R pjq δ v 0 . Then, given n d successive rays, with respective directions v pjq , j " 0, . . . , n d emanating from an initial radiating point v 0 at angular intervals of δ x , δ y and δ z (with δ x,y,z " 2π{n d ), the workspace boundary can be numerically mapped by solving the following optimization problem: ru bj , q bj s " arg min u,q }p E pqq ´vpjq } 2 2 s.t. u P U and (3) [START_REF] Della Santina | Dynamic control of soft robots interacting with the environment[END_REF] where the achieved actuation vector u bj and its associated displacement vector q bj are the optimal and feasible solution to minimize the distance between v pjq and p E . Accordingly, the associated end-effector boundary position can be deduced as:

p bj E " Cq bj `pp0q E .

Non-convexity issues and proposed solution

The process of mapping the workspace boundary BW E may sometimes be confronted with non-convexity difficulties which might lead to an overestimation of the workspace, as depicted by the following illustration (refer to Fig. 6). When this situation occurs ( }p

bn 1 E ´pbn 2 E } 2 ą r ), it is possible that either p bn 1 E or p bn 2 E
will be the only achievable solutions of ( 5) for all direction vectors v pjq between v pn1q and v pn2q , i.e., @v pjq P rv pn1q v pn2q s, (5) ùñ p bn The phenomenon yielding such a problem can be tracked down by recognizing an abrupt jump resulted from ( 5) between two successive iterations of at least one actuation u bj k P u bj . Accordingly, once such a situation is detected, this paper proposes a strategy that will allow us to completely identify the probably missed turquoise colored boundary points in Fig. 6. This technique focuses firstly on obtaining a list of radiating points v l by discretizing the space between p 7. Afterwards, we solve (4) for each point the obtained list v l , and get their related boundary points. Finally, we can successfully identify the boundary points lying on the non-convex part, and achieve the correct estimation of the workspace boundary (as illustrated in Fig. 7).

The proposed strategy can be summarized as follows:

Realization of the algorithm

The suggested optimization-based approach with the selected strategy to resolve the non-convexity issue is described in a concise algorithm (See Algorithm 1).

Algorithm 1 Calculate BW E Require: v 0 , C, U, n d , n s , r, p p0q E δ Ð 2π{n d ; BW E Ð ∅; Ź Initialization for j Ð 1 to n d do v pjq Ð R pjq δ v 0 Ź Direction vector ru bj , q bj s Ð solve (5) pv pjq , C, U, p p0q E q Ź Optimal solution p bj E Ð Cq bj `pp0q E Ź Boundary point BW E Ð BW E ' p bj E Ź Append if }p bj´1 E ´pbj E } 2 ą r then Ź Non-convexity v l Ð dispp bj´1 E , p bj E , n s q Ź Discretize for s Ð 1 to n s do ru bs , q bs s Ð solve (4) pv psq l , C, U, p p0q E q p bs E Ð Cq bs `pp0q E Ź Boundary point BW E Ð BW E ' p bs E Ź Append end for end if end for

VI. VALIDATION AND ANALYSIS

In this section, we illustrate the presented results by evaluating the reachable workspace of a trunk-like soft robot with a Young's modulus E " 1.8 ˆ10 7 P a, and a Poisson ratio ν " 0.45. The following presented scenarios were conducted on an Intel Xeon(R) processor with 16 GB of RAM and a 3.50 GHz clock speed. We investigate three situations in which a trunk-like soft robot is actuated by various tendon routing configurations, and we use simultaneously the forward and optimization-based approaches to estimate the accessible workspace W E and the exterior workspace boundary BW E for every case. For each situation, the outcomes of both techniques are superimposed and shown by a three-view graphic (Oxy, Ozy, and Ozx plans). We suggest discretizing the actuation vector with a step size of 1 Newton (unit) for the forward method, and discretizing the angles with a discretization step size of 0.08 Radian for the optimization-based approach.

A. Scenario 1: 2 tendons-driven trunk-like soft robot

The first scenario is the trunk-like soft robot operated by two tendons (as illustrated in Fig. 8a). We exert a force with a magnitude within U p0q " U " r0 100s ˆr0 100s.

Using the explanation outlined in Section IV-A, we first use the forward approach with n d " 100 to obtain the related estimation of W E for this soft robot's configuration (represented by blue colored points in 8b, 8c and 8d). The forward technique needs n nu d " 100 2 operations in terms of operation complexity, and it takes about 1463 seconds to compute the complete estimation of this scenario workspace in terms of time computation.

For the same scenario, the optimization-based approach required a total number of 864 operations with a time computation of 984 seconds. The obtained boundary points are depicted by green colored points in Figs. 8b, 8c and8d, which represent the exterior workspace boundary BW E of this soft robot's configuration. The second scenario considers the trunk-like soft robot operated by three tendons (as illustrated in Fig. 9a). We apply a force magnitude within U p0q " U " r0 100s ˆr0 100s ˆr0 100s.

The forward approach is firstly applied with n d " 100 to obtain this configuration's workspace W E (illustrated in Figs. 9b,9c, and 9d by blue-colored points). The forward technique needs n nu d " 100 3 operations, and 1.4 ˆ10 5 seconds to achieve the entire estimation of this scenario's workspace. Next, using Step. 3 (Fig. 5) of the optimization-based approach with n dz " 72 and n dx " 12, we map the boundary points of this particular configuration. For this scenario, the proposed optimization-based approach required a total number of 864 operations with a time computation of 1034 seconds. The achieved boundary points are illustrated by green colored points in Figs. 9b, 9c, and9d, which represent the workspace boundary BW E of this particular configuration.

C. Scenario 3: 4 tendons-driven trunk-like soft robot

The third scenario considers the trunk-like soft robot operated by four symmetric tendons (as illustrated in Fig. 10a). We apply a force whose magnitude within a specific range U p0q " U " r0 100s 1 ˆ¨¨¨ˆr0 100s 4 .

The forward approach is firstly applied with n d " 100 to obtain this configuration's workspace W E (illustrated by blue colored points in Figs. 10b, 10c and10d). The forward technique needs n nu d " 100 4 operations, and 1.4 ˆ10 7 seconds to completely estimate this scenario's workspace. Then, we use Step. 3 (Fig. 5) of the optimization-based approach with n dz " 72 and n dx " 12 to map the boundary points of this particular configuration. For this scenario, the optimizationbased approach required a total number of 864 operations with a time computation of 1039 seconds. The obtained boundary points are depicted by green colored points in Figs. 10b, 10c and10d, which represent the workspace boundary BW E of the studied soft robot's configuration.

D. Notes on scenario results

Based on the achieved results, we can clearly see from the three-view drawing of each scenario that the workspace boundary estimation result achieved via the optimization-based method encloses the estimated workspace obtained via the forward approach. Furthermore, when matched to the forward technique, the suggested optimization-based methodology may achieve equivalent workspace estimation accuracy while decreasing the number of operations and the amount of time it takes to estimate the workspace for each scenario.

E. Notes on computation complexity

We summarize the operations complexity and time calculation for every scenario in the following tables (see Tables I andII), where the time computation of the forward method for configurations with n u ě 3 is merely estimated since they need a very lengthy time. In addition, we can also illustrate the complexity of operations and the time computation by Figs. 11a and Fig. 11b, where it can be seen how expanding the dimension of the actuators has an exponential effect on the forward method, but the suggested optimization-based approach's is essentially linearly stable.

F. Notes on limitation

This paper proposed an optimization-based approach to efficiently estimate the workspace boundary of soft robots. The proposed optimization-based approach consists of mapping the surfaces enclosing the workspace, i.e., this approach can only successfully identify the points that constitutes the exterior boundary of the workspace. Although this approach has such a limitation, however as shown in [START_REF] Amehri | Discrete cosserat method for soft manipulators workspace estimation: an optimizationbased approach[END_REF], the optimization-based approach can effectively and successfully map the exterior boundaries of multisection soft robots actuated by a generic dimension of tendons.

VII. CONCLUSIONS

This paper proposed two different workspace estimation approaches for soft robots modeled by the FEM method. On one side, we provided a classic strategy, namely the forward approach, which discretizes the actuators' space in order to estimate the workspace of a soft robot. In contrast, we have suggested a specific optimization-based technique for reducing the forward method's computing cost, by mapping the exterior workspace boundary of soft robots and avoiding the costly calculation of inside configurations. Ultimately, utilizing the forward technique as a comparison reference, a validation of the suggested strategy was achieved using several configurations of soft robots. In comparison to the simple forward approach, which exponentially expands as the actuator's size rises, we have demonstrated the optimization-based approach efficiency.

Figure 1 :

 1 Figure 1: Soft robot and its accessible workspace.

Figure 2 :

 2 Figure 2: Geometric model and its associated FEM.

Figure 3 :

 3 Figure 3: Forward technique.

1 )

 1 Firstly, a definition of the workspace boundary, denoted as BW E , is established. 2) Secondly, a ray is emanated from a radiating point v along a certain direction to find a boundary point p b E on BW E . 3) In the end, we may get a map of BW E by finding p b E for many consecutive rays with appropriate directions originating at angular intervals.

Finally, as

  opposed to the forward technique, the suggested optimization-based strategy offers a number of benefits: 1) Since the end-effector space is always less than or equal to 3 (since we are only considering position access), it can reduce the dimension of the space to be discretized; 2) It results in reducing the total number of iterations required to estimate the workspace since it only focuses on finding then mapping the boundary points p b E that encloses W E . The process for the suggested optimization-based approach is summarized in the steps outlined below: [Step 1] Establishing the workspace boundary definition BW E ; [Step 2] Determining one point on BW E ; [Step 3] Achieving a map the workspace boundary BW E . V. REALIZATION OF THE OPTIMIZATION-BASED APPROACH This section will explain the implementation of the proposed optimization-based strategy.

Figure 4 :

 4 Figure 4: Calculate related actuation u b for a radiating v using (4), then derive associated boundary point p b E using (3).A boundary point p b E in a specific direction emerging from a radiating point v that is external to the workspace W E (as represented by Fig.4) may be found in accordance with Definition 2 by solving the constrained optimization problem defined below: ru b , q b s " arg min u,q

Figure 5 :

 5 Figure 5: Calculate related actuation u bj for a radiating v pjq using (5), then derive associated boundary point p bj E using (3). Precisely, with v 0 being the initial radiating vector that can be calculated based on the technique proposed in Step 2, and the rotation matrix R pjq δ is established by: R pjq δ " R δ pjq x R δ pjq y R δ pjq z , with R δ pjq x , R δ pjq y , and R δ pjq z

1 E or p bn 2 E

 12 As a result, we may miss the identification of the turquoise colored boundary points in Fig.6

Figure 6 :

 6 Figure 6: Non-convexity BW E : sudden jump from p bn 1 E to p bn 2 E .

  [Step. A] Calculate the set of v l by discretizing the space between p B] Compute the related feasible solution by solving (4) for each discretized point v l , then deduce its corresponding boundary point. To clearly understand the above steps, we present a successive iteration in Fig.7to illustrate the proposed strategy.

Figure 7 :

 7 Figure 7: Non-convexity solution: for v l , calculate its related actuation via (4) and deduce its associated boundary via (3).

  (a) Soft Robot with nu " 2. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 8 :

 8 Figure 8: Scenario 1: Forward method to get W E (blue points) and optimization-based method to get BW E (green points) .

  (a) Soft Robot with nu " 3. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 9 :

 9 Figure 9: Scenario 2: Forward method to get W E (blue points) and optimization-based method to get BW E (green points) .

  (a) Soft Robot with nu " 4. (b) Oxy plan view. (c) Ozy plan view. (d) Ozx plan view.

Figure 10 :

 10 Figure 10: Scenario 3: Forward method to get W E (blue points) and optimization-based method to get BW E (green points) .

Figure 11 :

 11 Figure 11: Comparison of operations complexity and time computation for forward method and optimization method.

Table I :

 I Comparison of operations complexity.

		Operations
	Inputs Approach	Forward	Optimization
	nu " 2	10 4	864
	nu " 3	10 6	864
	nu " 4	10 8	864

Table II :

 II Comparison of computational time.

		Time (seconds)
	Inputs Approach	Forward	Optimization
	nu " 2	1463	984
	nu " 3	« 1.4 ˆ10 5	1034
	nu " 4	« 1.4 ˆ10 7	1039