Ke Wu

Gang Zheng

Solutions to Large Beam-deflection Problems by Taylor series and Padé approximant for Compliant Mechanisms

Keywords:

Solutions to large beam-deflection problems by Taylor series and Padé approximant for compliant mechanisms

Ke Wu, Gang

Introduction 1.Compliant mechanisms

Compliant Mechanisms (CMs) have presented several desired properties in mechanical applications, gradually becoming a hot research area in recent years [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. Different from rigid mechanisms, CMs demonstrate a new concept for energy (including force and motion) transferring that merely depends on the elastic deformation of the builtin flexible members [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF]. This novel concept therefore enables CMs to present several more promising features than rigid ones [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF], such as energy-saving actuation [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF], simplified fabrication process [START_REF] Farhadi Machekposhti | A review on compliant joints and rigid-body constant velocity universal joints toward the design of compliant homokinetic couplings[END_REF], simple mechanical maintenance [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF] and more precise motion under the same mechanism dimensions [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]. These mentioned advantages have been demonstrated in many CM-involved applications [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [START_REF] Farhadi Machekposhti | A review on compliant joints and rigid-body constant velocity universal joints toward the design of compliant homokinetic couplings[END_REF], for instance, MEMs [START_REF] Sridhar Kota | Design of compliant mechanisms: applications to mems[END_REF], aircraft engineering [START_REF] Sridhar Kota | Smart structures and materials 2003: industrial and commercial applications of smart structures technologies[END_REF], surgical applications [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], high-precision motion stage [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF] and so on. In particular, Howell and his coworkers have done some impressive work on Origami mechanisms using CMs [START_REF] Greenwood | Conceptualizing stable states in origami-based devices using an energy visualization approach[END_REF], which gradually becomes one of the hottest and newest areas in the academic field of CMs. Therefore, to properly design CMs, accurately and efficiently modeling slender beams is of great importance [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF] where there are 2 major aspects: first, the fundamental modeling theory should be properly chosen; second, and the corresponding analytical or numerical methods that can implement the mentioned theories are highly needed.

The essence of modeling slender beams

As stated in the previous section, we would like to first stress and clarify that the fundamental modeling theory used for characterizing beam deflection and the corresponding methods to solve the derived governing equations are totally two different concepts. For example, beam theory is the theory used for characterizing the elastic behavior of slender beams while elliptical integral method [START_REF] Bulirsch | Numerical calculation of elliptic integrals and elliptic functions[END_REF] is just one method to solve the governing equations behind the theory. Similarly, the so-called Finite Element Method (FEM) used in this field essentially means using FEM to solve the governing equations behind solid mechanics or beam theory where solid mechanics or beam theory is the used theory and FEM is the computational method [START_REF] Singiresu | The finite element method in engineering[END_REF]. Note that seeking for effective theories and corresponding efficient methods to model flexible slender beams is always essential since they can significantly improve the accuracy and efficiency in topology model-based optimization of CMs [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF] .

Therefore, in the following subsections, we mainly review and also recall the essence of Euler Bernoulli beam theory since it's the most commonly-dependent, simplest and straightforward theory in CMs. Likewise, the commonlyused methods to handle it will be briefly reviewed and discussed accordingly.

Euler Bernoulli beam theory

The famous Euler Bernoulli beam equation implies that the curvature d d ( ) of the deflected slender beam at some point , or equivalently point ( ( ), ( )) is proportional to the external moment ( ) exerted at the same point or equivalently point ( ( ), ( )):

d d = ( ), ∈ [0, ] (1)
where ( ) refers to the rotation angle of the point with respect to X axis as shown in Fig. 1; is the Young's modulus of the material; denotes the second moment of inertia of the cross-section area; denotes the beam length; , and are the constant forces and moment exerted at the beam-end point = . The graphical explanation is presented in Fig. 1.

Then, considering the varying cross-section along the beam axis and an initial curvature of the slender beam, we can arrive at:

d d = ( ) ( ) + 1 ( ) (2)
where ( ) denotes the second moment of inertia of the cross section area at point , which implies the cross section may vary along the beam axis; similarly, 1 ( ) represents the initial curvature of the slender beam implying that the initial curvature may vary along the beam axis as well.

According to the definition of ( ) in Eq. ( 1), ( ) is the exerted moment at any arbitrary point or equivalently ( ( ), ( )) on the beam due to all external loads on the beam. Therefore, ( ) can be formulated as:

( ) = ( ( ) -( )) + ( ( ) -( )) + (3)
Rearranging Eq. ( 2) and Eq. ( 3), we will end up with:

( ) d d = ( ( ) -( )) + ( ( ) -( )) + + ( ) ( ) (4)
Differentiating [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF] to eliminate the constant variables ( ), ( ) and , we can obtain:

d d d d + ( ) d 2 d 2 = -( cos ( ) + sin ( )) + ( d d ( ) ( ) - ( ) d d ( ) ( ) 2 ) (5)
where d d = cos and d d = sin are utilized. Note that [START_REF] Sridhar Kota | Design of compliant mechanisms: applications to mems[END_REF] implies Euler Bernoulli beam theory can handle slender beams with varying cross sections and varying curvature along the beam axis, and there are many decent works regarding the aforementioned [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF][13] [START_REF] Henning | Non-linear analytical modeling of planar compliant mechanisms[END_REF]. However, for standard manufacturing process in CMs, there exist many contributions dealing with the following common cases, which are the simpler cases of (5): slender straight beams or initially curved beams (ICBs) with constant cross-section or constant initial curvature. Therefore, we can just narrow down the scope of the studied problem by simplifying (5):

D.E. d 2 d 2 = -( cos ( ) + sin ( )) (6)
D.E. d 2 d 2 = - (cos ( ) + sin ( )) B.C. (0) = 0 d d ( ) = + 1 (7)
where 1 → 0 means the studied beam has no initial curvature (straight beam). Note that nearly all existing commonly used methods to handle beam-deflection problems are based on [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] [9] [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], and these methods are briefly reviewed in the next section. For a more comprehensive review about modeling methods for large deflection of slender beams in CMs, [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF] is recommended.

Commonly used modeling methods for large deflection of slender beams in CMs

To handle [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], CM researchers have found out many effective methods, which are mainly classified into 2 categories: analytical methods and numerical methods.

1. Analytical methods: elliptical integral method was introduced by Howell [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] to solve BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Theoretically, this method is the most accurate one since it's analytical. However, in practice, it is implemented in a semianalytical manner since some complex integrals need to be numerically calculated.

Remark 1. Pseudo rigid body model (PRBM) [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF], also proposed by Howell, essentially introduced a unique modeling concept where the deflection of a slender beam can be physically represented by that of a rigid mechanism. It serves as an effective modeling option at the early stage of design process. Note that it is definitely not the analytical solution to the original BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] although it has a closed-form formulation.

2. Numerical methods: BCM [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF], proposed by Awtar, serves as an effective method of solving beam-deflection problems only within intermediate range where the curvature d d defined in ( 7) is linearized and rearranged as

d 2 d 2
in a global frame for easier and more straightforward solution. Dr. Chen then extended this work and proposed Timoshenko Beam Constraint Model (TBCM) for the purpose of characterizing stubby beams (which takes the influence of shear force into account) [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF]. To handle large-deflection problems, several chained algorithms are available: chained BCM [START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[END_REF] and chained PRBM [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF][START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF][START_REF] Yu | A pseudo-rigid-body 2r model of flexural beam in compliant mechanisms[END_REF][START_REF] Chen | Finding the optimal characteristic parameters for 3r pseudo-rigid-body model using an improved particle swarm optimizer[END_REF][START_REF] Zhu | Pseudo-rigid-body model for the flexural beam with an inflection point in compliant mechanisms[END_REF][START_REF] Yu | A novel model of large deflection beams with combined end loads in compliant mechanisms[END_REF][START_REF] Jin | A cprbm-based method for large-deflection analysis of contact-aided compliant mechanisms considering beam-to-beam contacts[END_REF][START_REF] Jr | A 3-d chain algorithm with pseudo-rigid-body model elements[END_REF] (see Remark 2 for details about chained-PRBM). Note that using chained algorithms to model CMs inevitably leads to high dimensional models, which is logically time-consuming to solve. FEM is often used to analyze large beam-deflection problems in CM research community. Essentially speaking, it's a method for numerically solving differential equations in engineering [START_REF] Singiresu | The finite element method in engineering[END_REF]. Therefore, FEM can handle the governing BVP (7) derived from Euler-Bernoulli beam theory, and it can also deal with the governing equations behind solid mechanics that describe beam deflection. However, FEM develops high-dimensional computational models, which may result in high computational expenses, for example, using solid mechanics to handle 3D modeling. In particular, using FEM to solve planar linear beam theories may still present desired efficiency in some cases [13][14]. This is why many CM researchers tend to use FEM-based commercial software (such as ANSYS, COMSOL and ABAQUS) only as a verification method but not quite a commonly used option for model-based optimization of CMs [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF]. Mathematically speaking, [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] is essentially a BVP of an ordinary differential equation (ODE) so weighted residual methods (denoted as WRM in this paper), finite difference method and shooting method [13][14] can therefore be utilized to approximate the solution of BVP (7) [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF].

Remark 2. Chained-PRBM artfully transforms the static beam deflection into a kinematic motion of a rigid mechanism where a limited number of rigid links associated with stiffness joints to approximate the solution of [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Basically, in the current literature, there are mainly two types of Chained-PRBM.

(a) The first type is the classic one where a small number of rigid links are used always with the need of estimating model parameters. For example, [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] provides a closed-form formulation for analyzing large deflection, which is significantly useful in the process of first-stage design. However, it suffers from the limitation of predicting larger end slope of flexible beams. The method proposed in [START_REF] Venkatasubramanian | Pseudo-rigid-body models for circular beams under combined tip loads[END_REF] is also efficient since only one unknown needs to be solved (see Eq. 13 in [START_REF] Venkatasubramanian | Pseudo-rigid-body models for circular beams under combined tip loads[END_REF]) using Newton's method or its relatives. However, it also needs to estimate the model parameters in advance. To us, this method is more like using a large number of loading cases to train the best parameters for the model via optimization (see Eq. 15 in [START_REF] Venkatasubramanian | Pseudo-rigid-body models for circular beams under combined tip loads[END_REF]). This type of 'learning' strategies is common in the field of machine learning. Only if the number of training cases (loading cases) approach infinite, the final model will always converge to the correct or accepted solutions under given convergence criteria [START_REF] Carleo | Machine learning and the physical sciences[END_REF]. (b) [START_REF] Jin | Large-deflection analysis of general beams in contact-aided compliant mechanisms using chained pseudo-rigid-body model[END_REF] is another typical type. If the number of the rigid links is large enough or approaches infinite, this Chained-PRBM can accurately approximate the solution of (7) essentially since:

( ) = = (8)
where is the length of each rigid link and = with being the beam length. If the number approaches infinitely large, equivalently, approaches infinitely small, we have

( ) = = lim →0 = d d (9)
with the boundary conditions available. Then, this method directly goes back to numerically solve BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] with reliable accuracy. However, due to the large number of links discretized, it ends up with a large number of unknowns to be numerically solved. As a result, the group of equations derived from the principle of minimal potential energy or Euler equation (moment equilibrium) can be solved by Newton's method or its relatives. The high dimension of the Jacobian matrix may result in expensive computational cost when calculating its inverse [START_REF] Carl | Solving nonlinear equations with Newton's method[END_REF].

As summarized in the current literature, previous numerical frameworks, such as beam-theory-based chained-algorithms and solid-mechanics-based or beam-theory-based FEM, for characterizing the large deflection of slender beams need to solve a large number of unknowns [START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[END_REF][19], which can be regarded that they all aim to locally approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. To increase the accuracy of the results, we need to increase the number of discretization to a certain level, which however may result in more computational expenses. In the recently published work [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], the authors proposed to use several weighted residual methods to numerically solve Euler Bernoulli beam equation, which is essentially using a high-order polynomial with unknown coefficients to approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. In [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], 6th-order polynomials are used so still 7 unknowns need to be solved (although the number of unknowns has been reduced quite a lot). Normally, the unknowns that exist in a group of developed nonlinear equations are solved by Newton Raphson method and its relatives. In this process, the larger number of unknowns yields higher dimension of the Jacobian matrix and its inverse matrix, resulting in more computational expense [START_REF] Isaacson | Analysis of numerical methods[END_REF]. Especially, when applying Newton-type method to seek optimal values of those multi-variables might yield local minima, depending on the chosen initial conditions, which is a well-known issue when applying such an approach. One solution to overcome this problem is to reduce the number of unknowns to be solved. Motivated by this idea, in this paper, we propose to use high-order Taylor series method and rational functions (Padé approximant) but with only 1 unknown to approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. In particular, we have noticed some inspiring works where Taylor series method and Padé approximant are used to handle the relatives of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] but only certain simpler cases are considered also without taking practical mechanical applications like CMs into account [32][33]. Roughly speaking, we first construct a high-order polynomial as the approximate solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] via Taylor expansion at = 0. By using most of the information of ODE ( 6), all coefficients of such an approximation can be then written as functions of only 1 unknown. Similarly, we can also construct a rational function to do so as well. Finally, the only unknown variable of the 2 constructed functions can be solved simply by satisfying the boundary conditions of [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Compared to [32][33], this paper deals with more complex scenarios, such as varying curvature and varying cross sections, and also proves the feasibility of conducting mechanism synthesis for CMs using the proposed methods as well.

Main content of this paper

In this paper, we first present the derivation of using Taylor series method and rational functions (Padé approximant) to handle BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] respectively. Then, the advantages and limitations of the 2 presented methods are discussed in detail through a convergence analysis, followed by comparisons and verification using weighted residual methods on two beam-end loading cases (a slender straight beam and a slender circularly curved beam). Both in the convergence analysis and the two verification cases, BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] is nondimensionalized for more general comparisons. Note that in the verification part, weighted residual methods are used as the standard since they have been proved to be an accurate and efficient method to solve (7) [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF]. In the end, the proposed methods have been proved feasible and efficient to model CMs where three representative cases (two compliant parallelograms and one pre-buckled bistable mechanism) are studied with error analysis provided afterwards.

Solving large beam-deflection problems via Taylor series method and Padé approximant 2.1. Taylor series method

Taylor series is a simple and effective option to approximate a given analytic function within its radius of convergence [START_REF] Murray | Advanced mathematics[END_REF]. Noticing the form of BVP ( 7), it's handy and straightforward to calculate derivatives at = 0 considering the existence of trigonometric functions in the ODE along with one of the boundary conditions: (0) = 0.

Therefore, we construct a Maclaurin series to approximate the solution of BVP ( 7):

Θ( ) = ∑ =0 d d (0) ! (10)
where Θ( ) is noted as the approximate solution of [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] with order and all d d (0) derivatives are obtained via using chain rule on the governing differential equation of ( 7) in an iterative manner (see the appendix for the detailed deduction).

As it has been shown in the appendix that any th-order of d d (0) can be expressed as a function of one unknown variable d d (0), so for the sake of simplicity we note it as:

d d (0) = ( d d (0)) (11)
where the functions are given in the appendix. Consequently, we can arrive at the approximation of

d d ( ) expressed as a function of the unknown variable d d (0), i.e., dΘ d ( ) = d( ∑ =0 d d (0) ! ) d = ∑ =1 ( d d (0)) ( -1)! ( -1) (12)
Based on the above equation ( 12), 2 strategies can be adopted to implement Taylor series method to solve BVP (7) [START_REF] Abbasbandy | Analytic continuation of taylor series and the boundary value problems of some nonlinear ordinary differential equations[END_REF]:

1. "Explicit" strategy

In this strategy, the second boundary condition provided in [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]:

d d ( ) = + 1 is explicitly imposed, i.e., dΘ d ( ) = ∑ =1 ( d d (0)) ( -1)! ( -1) = + 1 (13)
which assumes that the series converges within [0, ] if → +∞. Logically, [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF] will be a sufficient condition to numerically determine d d (0) when is large enough. Besides, in [START_REF] Henning | Analysis of planar compliant mechanisms based on non-linear analytical modeling including shear and lateral contraction[END_REF] there is only one unknown: d d (0), and equivalently it is a low-dimensional model. Therefore, we can use Newton-Raphson method or other related methods to both directly and efficiently solve (13).

"Minimal" strategy

In this strategy, the auxiliary condition is obtained by imposing that the next or equivalently ( + 1)th term of the truncated Taylor series vanishes, i.e.,

+1 ( d d (0)) ( + 1)! ( +1) = 0 ( 14
)
where it is assumed that th-order Maclaurin series converges within [0, ], and then + 1th term possibly converges to 0, which is sufficient for numerically determining the value of d d (0). Obviously, we can solve this via Newton-Raphson method as well. Obviously, the complexity of ( 14) is less than that of (13).

Remark 3. This "minimal" procedure is not always feasible to all types of functions in terms of Taylor expansion (Maclaurin expansion). For example, if we expand sin( ) using a Maclaurin series, we will have the following formulation:

sin( ) = - 3 3! + 5 5! - 7 7! + 9 9! ...
Then, we will always have:

d (sin( )) d | =0 ! = 0 ( = 2, 4 , 6, 8, 10...)
Obviously, if the numerical solution of an ODE is sin( ) within its corresponding domain, it's very likely to achieve some condition like [START_REF] Henning | Non-linear analytical modeling of planar compliant mechanisms[END_REF] although the used Taylor series may not converge at all in its given domain.

As long as the value of d d (0) is obtained using either of the above strategies, we can logically approximate ( ) via (10).

Padé approximant (rational function method)

A Taylor series can often be accelerated quite dramatically (or turned from divergent to convergent) by being rearranged into a ratio of two such series. Supposing we have approximated the solution of BVP (7) using Maclaurin series shown in [START_REF] Bulirsch | Numerical calculation of elliptic integrals and elliptic functions[END_REF], we can then construct a rational function via:

Θ( ) = ∑ =0 d d (0) ! = ∑ =0 = ( ) ( ) = ∑ =0 ∑ =0 (15)
where 0), ( ) and ( ) are polynomial functions of respective order and . ( 15) is called the Padé approximant [ , ] of Θ( ). The unknown coefficients and of ( ) and ( ) are determined by [START_REF] Bulirsch | Numerical calculation of elliptic integrals and elliptic functions[END_REF]:

= 1 ! d d (
( ) = ( )Θ( ) + [ +1 ] (0) = 1 + = (16) equivalently, 0 + 1 + 2 2 + ... = 0 + 1 + 2 2 + ... 1 + 1 + 2 2 + ... ( 17
)
Multiplying up the denominator gives the following equivalent set of relations:

0 = 0 1 = 1 + 0 1 2 = 2 + 1 1 + 0 2 3 = 3 + 2 1 + 1 2 + 0 3 ... ( 18
)
With given (derived in Section. 2.1), each new line introduces 2 new unknowns, and . However, due to the constraint + = , then there will be as many equations as unknowns (ignoring all terms ( + +1 )). For instance, if = 2 and = = 1, then we have 0 = 0 , 1 = 1 + 0 1 and 0 = 2 + 1 1 , i.e., 3 unknowns with 3 equations. Consequently, we can solve all these unknown coefficients and as functions of , equivalently functions of d d (0).

Remark 4. We need to take a look at the formulation of ( 15) first. In general, we have the following properties for rational functions:

if > , lim →+∞ Θ( ) = lim →+∞ ∑ =0 ∑ =0 = +∞ if = , lim →+∞ Θ( ) = lim →+∞ ∑ =0 ∑ =0 = if < , lim →+∞ Θ( ) = lim →+∞ ∑ =0 ∑ =0 = 0
Obviously, Θ( ) will gradually converge to a fixed value (which is called the horizontal asymptote of the rational function) if = . Next, rearranging (1), (2), and (3), we can obtain:

d d ( ) = ( ( ) -( )) + ( ( ) -( )) + + 1
It is easy to notice that d d ( ) will normally decrease and asymptotically reach + 1 if → . Therefore, ( ), which is formulated as ( ) = ∫ 0 d d ( )d will be logically bounded and also asymptotically reach a certain value as well if → . Therefore, given the similar properties of both the physical phenomenon and rational functions, this paper sets = = 2 where is chosen to be even for Padé approximant.

Next, differentiating ( ) ( ) :

d d ( ( ) ( ) ) = ′ ( ) ( ) -( ) ′ ( ) 2 ( ) (19)
and explicitly satisfying the second boundary condition = :

′ ( ) ( ) -( ) ′ ( ) 2 ( ) = + 1 ( 20
)
where there is only one unknown d d (0), it can be therefore easily solved via Newton-Raphson method. Then, we can finalize the formulation of the rational function [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF] to approximate the solution of (7):

Θ( ) = ∑ =0 ∑ =0 (21)
Remark 5. Considering the characteristics of rational functions, (21) may present the following 2 undesired properties in terms of characterizing or approximating ( ): [START_REF] Yu | A pseudo-rigid-body 2r model of flexural beam in compliant mechanisms[END_REF] will have a vertical asymptote at = ̃ . Note that either of the above scenarios contracts the continuum hypothesis stated in solid mechanics. We can ignore the hole when numerically approximating ( ) in terms of the first. For the latter scenario, we need to change and (normally increase them) to avoid the vertical asymptote.

1. if ∑ =0 ̃ = 0 and ∑ =0 ̃ = 0, (21) will have a hole at = ̃ . 2. if ∑ =0 ̃ = 0 and ∑ =0 ̃ ≠ 0,

Inflections points and nonlinear (post-)buckling 2.3.1. Inflections points along a deflected beam

According to [START_REF] Kimball | Modeling of flexural beams subjected to arbitrary end loads[END_REF] and [START_REF] Zhang | A comprehensive elliptic integral solution to the large deflection problems of thin beams in compliant mechanisms[END_REF], the essence of inflection points in a beam is all about the sign change of d d . If no inflection point exists, the slope ( d d ) is either monotonically increasing or decreasing, it must always be positive or negative, respectively. If an inflection point exists, the sign of d d changes once along the range of the independent variable ∈ [0, ] . Mathematically speaking, if d d ( * ) = 0 and d 2 d 2 ( * ) ≠ 0, then = * is an inflection point. For example, in nonlinear (post-)buckling of slender straight beams, inflection points are common (see Fig. 2a). The proposed methods numerically approximate the solution of BVP ( 7) so the inflection points will be implicitly treated by the proposed numerical methods while seeking the solution of BVP (7).

Nonlinear (post-)buckling

Note that there is a special type of beam-deflection scenarios called nonlinear buckling and post buckling. As shown in Fig. 2a, a slender straight beam fixed at one end is subjected to an axial force . Its governing equation can be rearranged from BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] as follows:

D.E. d 2 d 2 = - (cos ( ) + sin ( )) = - sin ( ) B.C. (0) = 0 d d ( ) = 0 (22)
Physically speaking, if reaches the buckling limit, the straight slender beam will experience nonlinear buckling, possibly resulting in 3 different (deformed) beam shapes shown in Fig. 2a. For nonlinear post bucking of the studied slender beam, if we exert the same beam-end loading and after the slender beam buckles (for example, Geometry 1 and Geometry 3 in Fig. 2a), the slender beam will also end up with 2 possible equilibrium states (Geometry A and Geometry B in Fig. 2b). The corresponding governing BVP that describes the post buckling behavior of the slender beam can be formulated as:

D.E. d 2 d 2 = - (cos ( ) + sin ( )) B.C. (0) = 0 d d ( ) = (23)
which is exactly the same as the original BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Mathematically speaking, it is normal that BVPs may have multiple solutions given some specific conditions. Normally, in applied mathematics, the concept of bifurcation theory [START_REF] Chow | Methods of bifurcation theory[END_REF] needs to be introduced to solve these BVPs where some techniques are used, such as adjusting initial guesses, determining the boundary conditions by making most of the information of the governing equation, perturbation method and so on. For example, the proposed Taylor series method and Padé approximant can handle the 2 scenarios (buckling and post-buckling) via directly choosing the proper initial guess. This is quite handy than using weighted residual methods [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF][39] since we only need to deal with just one initial guess of the unknown. Apparently, this specific topic is out of the scope of this paper, and will be further studied in our future work.

Characterization of beam shape

Finally, we can characterize the deformed beam shape via:

( ) = ∫ 0 cos(Θ( ))d ; ( ) = ∫ 0 sin(Θ( ))d ; (24)
for ∈ [0, ] with being the length of the deformed beam.

Convergence analysis and numerical validation

According to [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], it is obvious that modeling large deflection of slender beams is essentially a boundary value problem. In the first part of this section, we focus on the convergence analysis of Taylor series method and Padé approximant. Next, numerical testing is conducted to prove the feasibility of the used methods.

Prerequisites

In this paper, we first nondimensionalize the governing equation ( 7) for more general comparison. Therefore, the dependent variable and the independent variable are nondimensionalized via:

̂ = 1 ; = 1 ̂ ̂ = ; = ̂ (25)
Then, we can nondimensionalize (7) using (25) accordingly:

d 2 ̂ d ̂ 2 = - 2 cos ̂ ( ̂ ) - 2 sin ̂ ( ̂ ) (26)
with the corresponding nondimensionalized boundary conditions:

̂ (0) = 0 d ̂ d ̂ (1) = + (27)
Reformulating ( 26) and ( 27), we can arrive at:

D.E. d 2 ̂ d ̂ 2 = -cos ̂ ( ̂ ) -sin ̂ ( ̂ ) B.C. ̂ (0) = 0 d ̂ d ̂ (1) = + 1 (28)
where

= 2 ; = 2 ; = ; = ; ̂ ∈ [0, 1] (29)
Logically, we will have the following:

̂ ( ̂ ) = ∫ ̂ 0 cos( ̂ ( ))d ; ̂ ( ̂ ) = ∫ ̂ 0 sin( ̂ ( ))d (30)
where ̂ ( ̂ ) and ̂ ( ̂ ) are the nondimensionalized coordinates along the beam axis.

Remark 6. It is noticed that ODE [START_REF] Carleo | Machine learning and the physical sciences[END_REF] essentially governs all sorts of large deflection problems under the framework of nondimensionalization where all physical parameters can be represented by nondimensionalized values (see ( 29) and ( 30)). This strategy provides both feasible and logical prerequisites for Section 3.2 and Section 3.

Convergence analysis

Mathematically speaking, we use a Taylor series or Padé approximant of limited order to approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Obviously, it is not necessary to use a very high order since that may increase the complexity of numerical calculation. Moreover, a relatively low order may still present desired accuracy and convergence as well as high efficiency in its given domain. In general, we can implement an adaptive algorithm to determine the proper order of the Taylor series and Padé approximant: if |Θ ( ) -Θ -1 ( )| ≤ where is the th option ( = 1, 2, 3...) for and Θ ( ) is the approximation using a polynomial of th order, then is settled: = ; if not, = +1 = + 1. Likewise, a similar strategy can be implemented for Padé approximant as well. However, as shown in [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], a polynomial with a well-chosen fixed order can already accurately approximate the solution of [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] in terms of modeling large deflection. Therefore, we are interested in determining the order of the feasible Taylor series or Padé approximant in advance to bypass the adaptive algorithm for more straightforward feasibility.

Besides, given the same order, Taylor series and Padé approximant may present different convergence properties (such as advantages and limitations) in terms of approximating the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Therefore, being aware of

Taylor series method

As mentioned before, we are interested in setting up the proper order of the used Taylor series in terms of modeling large deflection problems. The deflections of a slender straight beam and a slender circularly curved beam are discussed respectively. The convergence of the series is of great importance in this method so the results of different orders will be analyzed in detail.

1. Large deflection of a slender straight beam under beam-end loading: in Fig. 3, as the order of the used Taylor series increases, the final deformed beam shape converges to a fixed curve (marked as convergence line) where in particular 14th-order series and 16th-order series almost cover each other and both converge to the convergence line. Generally speaking, a higher order of Taylor series normally can better approximate the original function within its radius of convergence [START_REF] Dienes | The Taylor series: an introduction to the theory of functions of a complex variable[END_REF]. However, for the studied problem, we have some extra constraints: two algebraic constraints represented by two boundary conditions and the ODE within the range ̂ ∈ [0, 1]. It is observed that in this case, the Taylor series subjected to the mentioned constraints turns out to converge with gradually decreasing oscillation. And during this convergence, low order Taylor series may happen to present 'better' precision than high order one, as depicted in Fig. 3. The results imply that given a proper order, Taylor series method can converge in terms of solving BVP (28) in this case. 2. Large deflection of a slender circularly curved beam under beam-end loading:

(a) Curvature-increasing deflection: in Fig. 4a, we cannot find solutions using 6th, 8th, 10th and 14th Taylor series even though the convergence criteria is set to be less restrictive. Moreover, the solutions found by the rest series do not present the trend of convergence at all when the order increases. Concluded from the results, Taylor series of even up to 16th order struggle to converge regarding solving BVP (7) of large anti-clockwise deflections of slender circularly curved beams. (b) Curvature-decreasing deflection: in Fig. 4b, on the contrary, Taylor series method presents great convergence in this case. The convergence trend of Taylor series is observed when the order increases, and series of orders up to 8th can converge to the convergence line well.

Remark 7. It is observed that Taylor series present different convergence properties in terms of curvatureincreasing and curvature-decreasing deflections. Mathematically speaking, the following are the possible reasons that make the Taylor series of a random function ( ) fail to converge.

(a) If ( ) is not infinitely differentiable, the Taylor series of ( ) may not even be defined.

(b) If ( ) is infinitely differentiable but not analytic, i. the derivatives of ( ) at the expansion point may grow so quickly that the Taylor series may not converge. Equivalently, the radius of convergence is 0 in this case. ii. The series may converge to something other than ( ). (c) If ( ) is analytic, its Taylor series about expansion point converges to the function in some neighborhood for every expansion point in its domain. However, the truncation error tends to grow rapidly away from expansion point if the order of the Taylor series is not large enough. In the case of curvature-increasing deflection, the constructed Taylor series expanded around ̂ = 0 struggles to converge. Due to the fact that ̂ ( ̂ ) increases more sharply (equivalently, the overall curvature gets larger within its domain) as independent variable ̂ increases, it can be logically concluded that (b) and (c) could be the major contributing reasons. On the contrary, in the case of curvature-decreasing deflection, ̂ ( ̂ ) increases more slowly (equivalently, the overall curvature gets smaller within its domain) as independent variable ̂ increases. Logically, the constructed Taylor series can more easily approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], which has been also numerically validated as shown in Fig. 4b.

Padé approximant

Likewise, choosing the proper order of the used Padé approximant is also of great importance. The deflections of a slender straight beam and a slender circularly curved beam are discussed in the following. The corresponding convergence analysis is provided.

1. Large deflection of a slender straight beam under beam-end loading: as shown in Fig. 5, Padé approximant of orders larger than 6th present great convergence in terms of approximating the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Padé approximant of 6th to 12th orders almost converges to the red dash line (convergence line) and cover each other closely. However, in particular, the curve of 6th Padé approximant has an abrupt change about (0.65, 0.45) (the red box) in Fig. 5. This is due to the fact that in the form of the approximating rational function, there exists a vertical asymptote in its domain ̂ ∈ [0, 1]. See the detailed analysis in Remark 5. Overall speaking, Padé approximant of relatively low order can still present great convergence in this case. 2. Large deflection of a slender circularly curved beam under beam-end loading:

(a) Curvature-increasing deflection: in Fig. 6a, it can be easily noticed that with the increase of the order, Padé approximant gradually converges to the red dash line (convergence line). With order up to 12th, its corresponding curve almost covers the convergence line. Therefore, we can draw the conclusion that it is feasible to achieve convergence using Padé approximant regarding solving BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] of large deflections of slender circularly curved beams. In particular, the curve of 8th Padé approximant has an abrupt change about (0.1, 0.59) (the red box) in Fig. 6a. The reason is explained in Remark 5. (b) Curvature-decreasing deflection: in Fig. 6b, Padé approximant gradually converges to the convergence line as the order increases. These rational functions of orders larger than 8th almost all cover the convergence line. In particular, the curve of 8th Padé approximant has an abrupt change about (0.42, 0.12) (the red box) in Fig. 6b. The reason is explained in Remark 5.

Comparisons between the mentioned methods

Generally speaking, concluded from the above results, Padé approximant presents better convergence than the corresponding truncated Taylor series in terms of solving BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. Specifically, Padé approximant performs well in both approximating the solution of the governing BVP for large deflection of straight beams and circularly curved beams even with a small number of series order. On the contrary, Taylor series method has been proved feasible only for analyzing large deflection of straight beams but with high-order series, and it struggles to converge in some cases of large deflection of circularly curved beams even under a high order. This is because Padé approximant makes most of the information of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] where the "shape" of the solution is pre-used by Padé approximant (the details are stated in Remark 4). However, Taylor series cannot be ensured to rapidly converge in all cases [START_REF] Abbasbandy | Analytic continuation of taylor series and the boundary value problems of some nonlinear ordinary differential equations[END_REF] where the details are stated in Remark 7. Compared to Taylor series method, Padé approximant has its unique limitations due to its formulation. As stated in Remark 5, the approximating rational functions might experience holes and vertical asymptotes in their domains, which contracts the continuum hypothesis of both solid mechanics and Euler Bernoulli beam theory. To handle this particular phenomenon, we can just increase the order of Padé approximant. From an engineering point of view, we choose TSM to analyze straight beams and PA to analyze circularly curved beams for further mechanism analysis.

Numerical validation

Based on the analysis stated in the last section, we then need to validate the proposed 2 methods through modeling a slender straight beam and a circularly curved beam under several beam-end loading scenarios where the deflected beam shapes and beam-end coordinates are studied. The graphical and numerical results are verified by beam-theorybased finite element method and compared by weighted residual method implemented in a global manner [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF] and Chained Beam Constrained Model [START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[END_REF]. In the following of this section, FEM, CBCM, WRM, TSM and PA stand for finite element method, chained beam constrained model, weighted residual method, Taylor series method and Padé approximant respectively. Besides, ER, CE and NOE denote the error with respect to the results of FEM, computational expense and number of elements in FEM and CBCM.

A slender straight beam under beam-end loading:

Case 1: = 1.3; = 0.9; = 0.4;

Case 2: = 0.65; = 0.45; = 0.2;

Case 3: = -0.65; = -0.45; = -0.2;

Case 4: = -2.6; = -1.8; = -0.8;

As shown in Fig. 7, the deformed beam shapes under different beam-end loading are graphically presented. The graphical and coordinate results of TSM and PA have both been verified by FEM and compared by WRM and CBCM in Fig. 7 and Table 1. In particular, compared to TSM, PA of equal or even smaller order can still present slightly better results. TSM and PA are both slightly more time-saving than FEM, CBCM and WRM in terms computational expense (see Table 1). As shown in Fig. 8, the deformed beam shapes under different beam-end loading are graphically presented. The graphical and coordinate results of TSM and PA have both been verified by FEM and compared by CBCM and WRM in Fig. 8 and Table 2 where TSM is not used for analyzing curvature-increasing defection of the circularly curved beam due to the reasons stated in Remark 7. Similarly, compared to TSM, PA of equal or even smaller order can still present slightly better results. TSM and PA are both slightly more time-saving than FEM,CBCM and WRM in terms computational expense (see Table 2).

A slender circularly curved beam under beam-end loading:

Discussions

Overall speaking, PA and TSM are both handy to deal with large deflection of slender beams (essentially solving BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]) except that TSM struggles to handle curvature-increasing deflection of slender circularly curved beams (essentially struggling to approximate the solution of BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] when its first-order derivative gets larger). It is interesting to notice that the rational functions generally present better accuracy than their truncated Taylor series in the case of solving BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. However, the comparison between PA and TSM is still an ongoing subject in applied mathematics, and the convincing reasons as far as we can provide is stated in Remark 4. Besides, the less computational expense of TSM and PA is due to the fact that only one unknown is solved in these two methods where Newton-Raphson method only needs to deal with one-dimension Jacobin matrix (that is the slope of the function).

Mechanisms synthesis using Taylor series method and Padé approximant

In this section, we aim to model 3 representative compliant mechanisms to prove the feasibility of Taylor series method and Padé approximant: straight-beam-based compliant parallelograms, circularly-curved-beam-based compliant parallelograms and mixed-type compliant parallelograms. These CMs are all of distributed compliance in their flexures. Not only verifying the feasibility of Taylor series method and Padé approximant for mechanism synthesis, we also aim to explore some desired characteristics of the aforementioned CMs for further optimization in our future work.

General methodology

For mechanism synthesis, we need to follow 3 modeling principles of mechanics [START_REF] Pieter | Advanced strength of materials[END_REF] These above 3 relationships govern the deformation of a flexible body. In other words, to model a CM, we need to develop the above three types of relationships, the corresponding equations of which need to be solved then through Newton-Raphson method. Finally, the modeling results are compared and verified by FEM results (where 3-D solid mesh element are used) with corresponding errors presented accordingly for further analysis.

Detailed modeling process and numerical results

Straight-beam-based compliant parallelogram

In this section, a straight-beam-based compliant parallelogram is modeled where Taylor series method is used to handle the constitutive beam equation, BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. As demonstrated in Fig. 9, the compliant parallelogram consists of two slender straight beams parallel to each other (beam #1 and beam #2) and a rigid motion stage. Each of the slender beams is fixed at one end respectively (O 1 and O 2 ) and rigidly attached to a rigid motion stage on the other end. , and ℎ refer to the length, width and thickness of the 2 involved elementary beams respectively. 2 local coordinate systems are set up for each flexible beam (see 9). The global coordinate frame X-O-Y is at the mid-point of O 1 and O 2 . A combined load composed of , and is exerted at the reference point K (the mid-point of AB), resulting in the deformed shape of the studied parallelogram mechanism shown in Fig. 9.

X 1 -O 1 -Y 1 and X 2 -O 2 -Y 2 in Fig.
, and are the forces and moment applied at the tip A of beam #1 by the rigid motion stage. , and are the forces and moment applied at the tip B of beam #2 by the rigid motion stage. Besides, A * , B * and K * refer to the corresponding points after deformation accordingly. Φ is the rotational angle of the motion stage. Note that O 1 O 2 = AB = for this parallelogram mechanism.

1) Constitutive relationships

Beam #1: D.E.

d 2 1 d 2 = - (cos 1 ( ) + sin 1 ( )) B.C. 1 (0) = 0 d 1 d ( ) = (31)
Beam #2: D.E.

B.C. 2 (0) = 0 d 2 d ( ) = (32)
2) Force equilibrium

- - = 0; - - = 0 - - + ∕2 sin Φ - ∕2 cos Φ - sin(Φ) + cos(Φ) = 0 (33)
3) Geometric compatibility relationships The applied loads are set up as:

= ∕2 , * * = ∫ 0 cos 1 ( )d ∫ 0 sin 1 ( )d + 0 ∕2 = -∕2 , * * = ∫ 0 cos 2 ( )d ∫ 0 sin 2 ( )d + 0 -∕2 AB = [( -) 2 + ( -) 2 ] 0.5 = A * B * = [( * - * ) 2 + ( * - * ) 2 ] 0.5 = | * - * | = sin Φ, | * - * | = cos Φ, 1 ( ) = 2 ( ) = Φ = 1 2 + 1 2 , * * = 1 2 * * + 1 2 * * , Δ Δ = * * - (34) Deformed Φ B A B* A* K K* P K M K F K P B F B P A F A M A M B Fabricated O X Y X 2 Y 2 O 2 X 1 Y 1 O 1 Front view B A K L #1 #2 h T O 1 O 2 K w Oblique
⎡ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎦ = = ⎡ ⎢ ⎢ ⎣ 125 N 25 N 12.5 N.m ⎤ ⎥ ⎥ ⎦ ( = 1, 2, 3, 4, ⋯ , 10)
The results of and Δ Δ Φ are provided via using the methodology stated in Section. 4.1 where the constitutive equations are solved by Taylor series method. These results are verified by solid-mechanics-based FEM verification (see Fig. 10). It is observed that Taylor series method presents high accuracy with errors (compared to FEM) all below 2%. In particular, the next-step initial guess is set up as the last-step answer for faster convergence [START_REF] Singiresu | The finite element method in engineering[END_REF]. The graphical results by FEM are shown in Fig. 14a where the fabricated and the 10 -deformed mechanisms are graphically presented.

Circularly-curved-beam-based compliant parallelogram

In this section, a circularly-curved-beam-based compliant parallelogram is modeled where Padé approximant is used to handle the constitutive beam equation, BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF]. As demonstrated in Fig. 11, the compliant parallelogram consists of two slender circularly-curved beams (beam #1 and beam #2) and a rigid motion stage. In particular, these two circularly-curved beams share the same center of circle. Each of the slender beams is fixed at one end respectively (O 1 and O 2 ) and rigidly attached to a rigid motion stage on the other end. Similarly, , and ℎ refer Front view

O 2 X 2 Y 2 P A F A M A A A* K* F K P K M K O X Y O 1 X 1 Y 1 K O 1 O 2 Oblique view B A #1 #2 R h w K O R R R
B.C. 1 (0) = 0 d 1 d ( ) = + 1 (35
4 ) = 0 + - - ∕2 sin(Φ) - ∕2 cos(Φ) - sin( ∕4 -Φ) + cos( ∕4 -Φ) = 0 (
3) Geometric compatibility relationships The applied loads are set up as:

= 2 sin( ∕4) sin( ∕4) , * * = cos sin -sin cos 0 1 1 0 ∫ 0 cos 1 ( )d ∫ 0 sin 1 ( )d + 0 sin( ∕4) = 2 sin( ∕4) -sin( ∕4) , * * = cos sin -sin cos ∫ 0 cos 2 ( )d ∫ 0 sin 2 ( )d + 0 -sin( ∕4) AB = [( -) 2 + ( -) 2 ] 0.5 = A * B * = [( * - * ) 2 + ( * - * ) 2 ] 0.5 = | * - * | = sin Φ, | * - * | = cos Φ, 1 ( ) -∕2 = ∕2 -2 ( ) = Φ = 1 2 + 1 2 , * * = 1 2 * * + 1 2 * * , Δ Δ = * * - (38)
⎡ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎦ = = ⎡ ⎢ ⎢ ⎣ 200 N 100 N 10 N.m ⎤ ⎥ ⎥ ⎦ ( = 1, 2, 3, 4...9)
The results of and Δ Δ Φ are provided via using the methodology stated in Section. 4.1 where the constitutive equations are solved by Padé approximant. These results are verified by solid-mechanics-based FEM (see Fig. 12). It is observed that Padé approximant presents high accuracy with maximum errors (compared to FEM) around 2%. Likewise, the next-step initial guess is set up as the last-step answer for faster convergence [START_REF] Singiresu | The finite element method in engineering[END_REF]. The graphical results by FEM are shown in Fig. 14b where the fabricated and the 9 -deformed mechanisms are graphically presented.

Besides, the load-bearing ability is of great importance in terms of analyzing the performance of a compliant parallelogram. Here, we take the circularly-curved-beam-based compliant parallelogram for example to analyze its load-bearing ability via Padé approximant where the relationship between and Δ is studied: The numerical results are provided in Fig. 13 with the maximum errors (compared to FEM) around 3.5%.

⎡ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎦ = = ⎡ ⎢ ⎢ ⎣ 0 N 100 N 0 N.m

Pre-bucked Bi-stable Mechanisms

Here, a special branch of CMs called Pre-bucked Bi-stable Mechanisms (PBM) [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF] are modeled here to prove the feasibility of the proposed method in mechanism synthesis. As shown in Fig. 15, this PBM is composed of a slender straight beam fixed at one end and connected with a rigid passive joint. The translational and vertical motions of the beam end is realized by the cooperation of two rollers. The initiation of PBM is Scenario 1. As is increased, the mechanism will end up with three possible states: still staying in Scenario 1; buckling into Scenarios 2 and 3. This special physical phenomena can be also described by Euler Bernoulli beam theory, which means we can use the proposed methods to handle the governing BVP. However, PBMs cannot be modeled through the normal procedure as presented in Section 2.2 since the existence of rigid joints in PBMs changes the governing equation of [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] and its boundary conditions. Therefore, we need to modify BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] according to the practical constraints (see Fig. 15 and Fig. 16):

d d = ( ) = ( ( ) -( )) (39)
Then differentiating (39) leads to:

d 2 d 2 = -sin ( ) (40)

Conclusions

Compliant mechanisms have been used as a promising option in many robotic systems where the built-in elementary flexible members play a dominant role in these mechanisms. Therefore, solving the constitutive equations of the elementary flexible members, such as slender straight beams and slender circularly-curved beams are a major focus. In this paper, we use Taylor series method and Padé approximant to solve the governing constitutive equations. Its accuracy and efficiency have been compared with weighted residual method and also verified by Finite Element Method (FEM) respectively. The feasibility of the proposed method has also been proved in terms of synthesizing CMs where three representative cases are studied. Our future work will focus on efficient numerical methods to solve the 3-D deflection of slender beams.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix: Iterative deduction of d d (0)

According to the governing equation [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], we can easily obtain the 2nd-order derivative at 0 as:

Figure 1 :

 1 Figure 1: General loading beam-end loading of a slender beam

1

 1 0 are used. Besides, 2 more boundary conditions are needed: (0) = 0 implying that the tangent angle at = 0 is zero since it is fixed; being the constant initial curvature of the studied flexible beam (circularly curved beams) which is derived from Euler-Bernoulli beam equation at = . Therefore, we will end up with a boundary value problem (BVP) of ODE (6):

Figure 2 :

 2 Figure 2: Schematic graph of nonlinear buckling

Figure 3 :

 3 Figure 3: Convergence analysis on large-deflection problem of a slender straight beam using Taylor series

Figure 4 :

 4 Figure 4: Convergence analysis on large-deflection problem of a slender circularly curved beam using Taylor series

  approximates 4th order (2-2) Padé approximates 6th order (3-3) Padé approximates 8th order (4-4) Padé approximates 10th order (5-5) Padé approximates 12th order (6-6) Padé approximates

Figure 5 : 1 . 1 2 .

 5112 Figure 5: Convergence analysis on large-deflection problem of a slender straight beam using Padé approximant

Figure 6 :

 6 Figure 6: Convergence analysis on large-deflection problem of a slender circularly curved beam via Padé approximant

  [START_REF] Robert R Archer | An introduction to mechanics of solids[END_REF][43]: a) Constitutive relationships b) Force equilibrium or force compatibility relationships c) Geometric equilibrium or geometric compatibility relationships

Figure 9 :Remark 8 .

 98 Figure 9: Constitutive, force equilibrium and compatibility relationships of a straight-beam-based compliant parallelogram and its schematic geometry

Figure 10 :

 10 Figure 10: Numerical results via Taylor series method

Figure 11 :

 11 Figure 11: Constitutive, force equilibrium and compatibility relationships of a circularly-curved-beam-based compliant parallelogram and its schematic geometry

  A, A * , B * , K and K * . Then, we can numerically solve Eqs. (31) to (34) using Newton-Raphson method. Besides, Δ Δ Φ refers to the translational and rotational displacements of K, the reference point of the modeled parallelogram mechanism. Modeling the studied revolute mechanism is all about the finding out the relationships between applied loads and the corresponding displacements Δ Δ Φ both at the reference point K. The material properties and geometry information are provided in the following: = 200 × 10 9 Pa; = 0.015 m; ℎ = 0.0045 m; = 0.25 m; = ℎ 3 12 ; = 2 ∕ ; = 2 sin( ∕4)

  [START_REF] Dienes | The Taylor series: an introduction to the theory of functions of a complex variable[END_REF], -39, ⋯ , 0, ⋯ ,[START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF][START_REF] Dienes | The Taylor series: an introduction to the theory of functions of a complex variable[END_REF]

Figure 13 :

 13 Figure 13: The relationship between and Δ

FabricatedFigure 14 : 3 Figure 15 :Figure 17 :

 1431517 Figure 14: Graphical results evaluated by solid-mechanics-based Finite element method

  = -(cos (0) + sin (0)) = -which in fact is constant.

  3.

	0.7	Undeformed beam shape Convergence line					Convergence line
		2th order Taylor series
		4th order Taylor series
	0.6	6th order Taylor series
		8th order Taylor series
		10th order Taylor series
	0.5	12th order Taylor series
		14th order Taylor series
		16th order Taylor series
	0.4
	0.3
	0.2
	0.1
						Undeformed beam shape
	0
		0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Table 1 :

 1 Comparison of beam-end-coordinates of a deformed straight beam

							Case Case 1 Case 2 Case 3 Case 4
								1.3	0.65	-0.65	-2.6
						Loading	0.9	0.45	-0.45	-1.8
								0.4	0.2	-0.2	-0.8
							̂ (1)	0.7157	0.9379	0.9767	0.8769
						TSM	̂ (1) ̂ (1)	0.6130 1.1633	0.3098 -0.1916 -0.4247 0.5324 -0.3357 -0.8255
							Order	16	8	4	16
	0.8	Taylor series method (TSM) Padé approximates (PA) Weighted Residual Method (WRM) Chained Beam Contraint Model (CBCM)		PA	̂ (1) ̂ (1) ̂ (1) Order	0.7145 0.6140 1.1663 10	0.9379 0.3098 -0.1918 -0.4264 0.9767 0.8757 0.5324 -0.3360 -0.8304 8 4 10
	0.6	Finite Element Method (FEM)				̂ (1)	0.7146	0.9377	0.9765	0.8756
						WRM	̂ (1)	0.6139	0.3104 -0.1924 -0.4266
							̂ (1)	1.1658	0.5336 -0.3371 -0.8311
	0.4		Case 1	Case 2	CBCM	̂ (1) ̂ (1) ̂ (1) NOE	0.7146 0.6139 1.1659 16	0.9376 0.3104 -0.1925 -0.4267 0.9765 0.8756 0.5336 -0.3372 -0.8311 16 16 16
	0.2						̂ (1)	0.7149	0.9377	0.9766	0.8756
						FEM	̂ (1) ̂ (1)	0.6139 1.1657	0.3104 -0.1925 -0.4269 0.5335 -0.3372 -0.8315
				Undeformed shape		NOE	15	15	15	15
	0						̂ (1)	0.11%	0.02%	0.01%	0.15%
					Case 3	ER TSM	̂ (1) ̂ (1)	0.15% 0.21%	0.19% 0.21%	0.47% 0.44%	0.52% 0.08%
	-0.2					ER PA	̂ (1) ̂ (1) ̂ (1)	0.06% 0.02% 0.05%	0.02% 0.19% 0.21%	0.01% 0.36% 0.36%	0.01% 0.12% 0.08%
	-0.4			Case 4		CE (s)	TSM PA WRM	0.008 0.010 0.020	0.009 0.009 0.027	0.009 0.010 0.012	0.011 0.010 0.017
	0	0.2	0.4	0.6	0.8	1	CBCM	0.20	0.30	0.25	0.35
							FEM	11	5	4	8
	Figure 7: Comparison of graphical results of a
	deformed straight beam

Table 2 :

 2 Comparison of beam-end-coordinates results of a deformed circularly-curved beam

  )

	Beam #2: D.E.	d 2 d 2 2	= -	(cos 2 ( ) +	sin 2 ( ))
			B.C. 2 (0) = 0				(36)
					d 2 d	( ) =	+	1
	+	cos(	4	) +	cos(	4	) -	cos(	4	) -	cos(	4	) = 0
	+	sin(	4	) +	sin(	4	) -	sin(	4	) -	sin

if > 0 and

if < 0, which clearly corresponds to Scenarios 2 and 3 as shown in Fig. 16. Here, we only take (45) for example to finish the rest deduction and logically solving (46) follows the same procedure. Assume ( ) as:

then we can obtain the value of by solving ( ) = 0 under Newton-Raphson method. So far, we have obtained all the information needed of the boundary conditions:

As stated in Section 2.2, we can also define the approximate solution of (40):

Therefore, we can arrive at the following relationships according to (48):

where similarly there is only one unknown d d (0) in the above equations. Then, we can use Newton's method to solve (50) directly. After obtaining the value of d d (0), we will be able to approximate the solution of ( 40) via (49). Then, we can calculate the coordinates of the beam end via:

In this modeling, we are interested in how the external force actuates the studied mechanism so the relationships of -( ), -( ) and -are studied (see Fig.

. Error analysis

Here, we aim to analyzing where the errors come from in the numerical models. According to Fig. 10, Fig. 12 and Fig. 13, Taylor series method and Padé approximant have proved effective in terms of mechanism synthesis (essentially dealing with the constitutive equations of slender straight beams and slender circularly-curved beams). Apart from the common numerical errors and model errors stated in [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], Taylor series method and Padé approximant have their unique numerical errors due to their intrinsic strategies of solving ODEs. Basically, both of them utilize the information of the studied ODE (where in this case BVP [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] is concerned) at one boundary to construct the approximate Taylor series or Padé approximant. Then, by satisfying another boundary, we will be able to find the proper series to approximate the solution to the studied ODE. However, in this process, the order of the used Taylor series or Padé approximant plays a rather important role. If it's not high enough, the solution will not converge to the correct answer (as shown in Fig. 3, Fig. 4, Fig. 5 and Fig. 6). Therefore, the numerical errors may also come from the improperly chosen order of the used Taylor series or Padé approximant. Logically, we can reduce these errors by increasing the order of the used Taylor series or Padé approximant.

Complexity of modeling

As we have stated in the appendix and according to [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], for a given order , we can always deduce the analytic formulation of the th derivative of at = 0, and all those high-order derivatives are the functions of only one unknown variable d d (0). The whole modeling process tends to be complex due to the deductions of the final nonlinear equation to be solved. However, it is worth noting that those tedious deductions can be automatically obtained via symbolic computation software, such as Maple. For practical coding, Maple provides the interface to transform the formulation into the commonly used coding languages, like Matlab, Python, C sharp, fortran, java and so on for further mechanism synthesis. Therefore, the complexity of modeling via the proposed Taylor series method and Padé approximant could be logically deduced using proper mathematical tools.

Then, by differentiating [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF] with respect to and take the value = 0, we can get the 3rd-order derivative at 0 as:

which indeed is a function of one unknown variable d d (0). Calculating again the twice differentiation of ( 7) with respect to and take = 0, we obtain the 4th-order derivative at 0 as:

which can be implicitly regarded as a function of one unknown variable d d (0), since it has been shown in the former calculation that d 2 d 2 (0) is a function of the unknown variable d d (0). By iterating the above procedure, we can obtain the successive high-order derivatives d d (0), which gives the 5th-order derivative at 0 as:

the 6th-order derivative at 0 as:

the 7th-order derivative at 0 as: the 8th-order derivative at 0 as: (0) = ......

In summary, according to [START_REF] Theodosia | Surgical applications of compliant mechanisms: A review[END_REF], for a given order , we can always deduce the analytic formulation of the th derivative of at = 0, and all those high-order derivatives are the functions of only one unknown variable d d (0). It is worth noting that those tedious deductions can be automatically obtained via symbolic computation software, such as Mathematica and Maple.


FEM-Based Nonlinear Controller for a Soft Trunk Robot

Ke Wu, Gang

FEM-based Nonlinear Controller for a Soft Trunk Robot

Ke Wu 1 and Gang Zheng 2 Abstract-Soft robots have been attracting much research attention in the past decades owing to their unique properties that are inherently different from conventional rigid robots. The intrinsic characteristics of soft robots, such as compliance and flexibility, have made them a promising supplement for many existing robotic applications. For instance, soft robots provide safer human-machine interaction and better adaptability to unexplored environments than the conventional rigid ones. However, their "soft" feature would conversely cause complex material deformation due to their high-DoF nature, which has made their kinematic or dynamic models more difficult to deduce. Logically, the process of designing model-based controllers would therefore suffer from the lack of efficient and accurate modeling methods. In this paper, through adopting Finite Element Method (FEM), we demonstrate a FEM-based nonlinear control method for a soft trunk robot to achieve trajectory tracking objectives. Basically, the proposed controller is developed by utilizing the essential information from the dynamic FE model of the studied robot. The details of deriving this nonlinear control method are demonstrated, then followed by 4 practical experiments on the soft truck robot to verify the feasibility of the control method. The corresponding experimental results are presented accordingly.

I. INTRODUCTION

Soft robots, distinguished from conventional rigid robots, have nowadays displayed many promising properties over conventional rigid robots, gradually bringing about thriving development in different industries such as manufacturing and specialized robotics [START_REF] Rus | Design, fabrication and control of soft robots[END_REF]. They are inherently flexible, compliant and lighter than conventional rigid robots, empowering themselves to change their deformation to transform into a more suitable configuration for a given task or an unfamiliar work environment [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Due to the mentioned advantageous features, these soft robots can achieve safe interaction with operators and environments [START_REF] Lee | Soft robot review[END_REF], [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], and they can also be powered by remote actuation with simple structure [START_REF] Gravagne | Kinematic transformations for remotely-actuated planar continuum robots[END_REF]. Obviously, these advantageous attributes have revealed the potential of soft robots to be integrated into many more industrial applications.

However, the beneficial properties due to their compliant nature can also lead to difficulties in modeling. The soft robot's complex material deformation results in nearly infinite degrees of freedom (DoF) and is essentially a highdimensional nonlinear system which is problematic to handle. The large number of designs including geometry and actuation have made every single soft robot unique, making it difficult to find a universal method to model them. From gang.zheng@inria.fr the perspective of controller design, the nonlinear hysteresis of these robots owing to compliance and flexibility also limits high-frequency control [START_REF] Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF].

To realize fast, reliable, accurate and energy-efficient control, many representative studies have been completed which can be mainly categorized into two types: modelbased controllers relying on analytical models and modelfree controllers based on machine learning techniques or empirical methods.

One of the most notable studies is a modified Denavit Hartenberg method called Constant Curvature (CC) assumption proposed in [START_REF] Hannan | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF] for mapping its configuration space to task space, which has been adopted in many static model-based controllers. In terms of multi-section continuum robots, the Piece-wise Constant Curvature (PCC) model logically becomes a recommended model partially owing to the success of CC model, which has as well been validated in many cases [START_REF] Jones | Kinematics for multisection continuum robots[END_REF]. Apart from the above mentioned, some more complex modeling approaches adopting beam theory [START_REF] Gravagne | Large deflection dynamics and control for planar continuum robots[END_REF] and Cosserat rod theory [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF] were investigated for soft/continuum robots as well. Logically, dynamic controllers cannot be developed purely based on the mentioned static models as they are not able to capture the dynamics of soft robots. In the current literature, the studies about developing dynamic model-based controllers for soft robots are rather limited. The authors of [START_REF] Gravagne | Uniform regulation of a multisection continuum manipulator[END_REF] have completed one of the first remarkable contribution regarding dynamic control of a planar multi-section continuum robot where a special exponential trajectory tracking has been achieved. Under the framework of PCC assumption [START_REF] Jones | Kinematics for multisection continuum robots[END_REF], the contributors of [START_REF] Della Santina | Dynamic control of soft robots interacting with the environment[END_REF] achieved dynamically control of a pneumatic robot in a closed-loop manner where they used an approximated dynamic model to design the dynamic controller. It is also worth noting that some representative studies on developing the kinematics and dynamics of pneumatically actuated continuum robots have been completed via model shape functions (MSF) [START_REF] Godage | Modal kinematics for multisection continuum arms[END_REF], [START_REF]Dynamics for variable length multisection continuum arms[END_REF]. Machine learning or empirical methods have been playing a gradually important role in rigid robots [START_REF] Nguyen-Tuong | Model learning for robot control: a survey[END_REF]. Likewise, these data-dependent methods have also been adopted for soft robotics. Representatively, the authors of [START_REF] Giorelli | A feed-forward neural network learning the inverse kinetics of a soft cable-driven manipulator moving in three-dimensional space[END_REF] took the lead in developing a static model-free controller, which can correctly predict the desired cable tensions to reach a target in simulations. Different from modelbased methods, model-free approaches can bypass the need of obtaining information from the robot itself, but depends on a large number of training data, performing well in systems that are highly nonlinear and complicated, easily influenced by gravity [START_REF] Rolf | Efficient exploratory learning of inverse kinematics on a bionic elephant trunk[END_REF], or operate nicely in unexplored environments where modeling is barely possible [START_REF] Yip | Model-less feedback control of continuum manipulators in constrained environments[END_REF]. However, due to Fig. 1: The studied silicone-made trunk robot [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] the uncertain nature of model-free approaches, stability and convergence analysis is challenging.

In this paper, Finite Element Method (FEM) [START_REF] Rao | The finite element method in engineering[END_REF] is chosen to model soft robots. The mentioned popular models are quite restricted by the geometry of soft robots and the models' quasi-static assumption. Notably, FEM has gradually turned out to be one of the desired approaches of modeling and developing controllers for soft robots [START_REF] Zhang | Vision-Based Sensing of External Forces Acting on Soft Robots Using Finite Element Method[END_REF] [22] [START_REF] Thieffry | Reduced order control of soft robots with guaranteed stability[END_REF] [24] as FEM does not require soft robots to have certain types of geometries, structures and quasi-static assumptions, serving as quite a universal strategy for modeling soft robots.

In our previous work [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF], we proposed a gain-scheduling control method in local manner where the gains of the linearized controller are scheduled from one sub-workspace to another to achieve operation within the whole workspace. Besides, only path following is concerned without stressing the dynamics or dynamic control that much in the process. However in this paper, as an extension of the former work, we propose to conduct FEM-based nonlinear dynamic control of a soft trunk robot in a global manner within its workspace to track trajectories as well as achieving more globally spatial and dynamic aims. The proposed method has been proved to be valid for enabling the robot to track the trajectory even under various disturbances, showing great robustness. The theoretical derivation and the experimental testing on a soft trunk-like robot are presented accordingly.

II. PROBLEM STATEMENT

A. Investigated soft manipulator

This paper aims to study the control of a silicone-made trunk robot with its exact geometry shown in Fig. 1. Particularly, to place a position sensor at the robot's end effector, we made a long and uniform hole along the robot's body from its base to its tip and mounted 4 inextensible cables for tendon-driven actuation. 15 3D-printed rigid rings are mounted along its body (see Fig. 1) to reduce the friction between the cables and the soft silicone. The material we used for shaping the robot is isotropic silicone rubber with the following material properties:

B. Control objective

The main work of this paper aims to control the endeffector of the trunk-like soft robot (see Fig. 1) to relatively quickly and stably track a given desired trajectory with (or without) the interfere of various bounded disturbances. Fig. 2: The mesh elements of the trunk robot [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] III. MODELING As is known to all, it is always difficult to dynamically or statically model soft robots since they are able to present unexpected irregular deformation due to the strong flexibility and compliance [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF]. Here, FEM is utilized to derive the dynamics of the studied soft trunk-like robot. Via synthesizing its model information from a FEM simulator (SOFA Framework), we can design a FEM-based dynamic controller for this robot.

In terms of deriving the FE dynamical model of the robot, readers are referred to [START_REF] Rao | The finite element method in engineering[END_REF] and [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF]. Essentially speaking, the procedure is similar to the classic FEM. First, the number and type of mesh elements need to be chosen to discretize the geometry of the modeled deformable robot (demonstrated in Fig. 2). By doing this, the displacement of any arbitrary point within mesh elements can be determined through shape functions, as well as the corresponding strain and stress. Finally, Euler-Lagrange equations are utilized to derive the FE dynamical model (Eq. 1) of the studied soft robot. M (q)q + D( q, q) q + K(q)q = H(q) T u

where M (q) 2 R 3n⇥3n stands for mass matrix; D( q, q) 2 R 3n⇥3n stands for damping matrix; K(q) 2 R 3n⇥3n stands for stiffness matrix; H(q) 2 R 4⇥3n stands for the input control matrix; q 2 R 3n⇥1 , q 2 R 3n⇥1 and q 2 R 3n⇥1 stand for the displacement, velocity and acceleration of all the n mesh nodes respectively in the robot's FE model; u 2 R 4⇥1 is an input vector that contains the four tension forces exerted at the four cables respectively.

Here, the function of H(q) is to only distribute u to all the elements containing the cable nodes (the red nodes in Fig. 2) since the friction between the cables and trunk robot is neglected in the derived FE model. M (q) is the corresponding mass matrix at q and can be reduced as a constant diagonal matrix due to its inherent seldom varying characteristics so that it will be replaced by M in the following sections for simplicity.

Remark 1: Using the above model ( 1), the hysteresis and other nonlinear effects can be also modeled by K(q), the dynamic friction could be modeled by D( q, q) and the static friction could be considered as a matched external disturbance f added to the channel u (i.e., replace u in (1) by u + f ). Logically, the accuracy of the developed model plays a very important role in simulation and openloop control. However, we would like to emphasize that it does not need a rather accurate model to achieve the control objective as long as the controller has access to essential model information, since those nonlinear effects might be eliminated or attenuated via robust controllers. We will discuss this issue in Section IV. In the robot's FE model, tetrahedron mesh elements are chosen and the robot's body is spatially discretized into n = 1484 nodes. Logically speaking, the model will have a larger number of mesh elements if finer mesh elements are chosen to discretize the robot, which leads to more accurate results but more computational time. The FE model of the robot is demonstrated in Fig. 2. Note that the dynamic position control of the robot's end effector (noted as q e 2 R 3⇥1 ) is our main focus. Therefore, under the framework of the derived FE model, this problem can be transformed into an output regulation problem of the following nonlinear dynamic system: M (q)q + D( q, q) q + K(q)q = H(q) T u

q e = Cq (2)
where C 2 R 3⇥3n is used to pick out the end-effector node from all n = 1484 nodes defined in the FE model. Logically, we can easily define the workspace of the robot's end effector as follow:

W e = {q e 2 R 3⇥1 | q = q = 0, q e = Cq} (3)
constraint by [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Identically, it can also be formulated as:

W e = {q e 2 R 3⇥1
| q e = Cq, K(q)q = H(q) T u, 8u 2 U } (4) where U denotes the set of all admissible inputs.

Given two points inside W e , noted as the initial position q init e 2 W e and the final position q fin e 2 W e , there exist different techniques in the literature, such as differential flatness [START_REF] Sira-Ramirez | Differentially flat systems[END_REF], interpolation via polynomials [START_REF] Lewis | Robot manipulator control: theory and practice[END_REF], to generate smooth trajectories. Therefore, this paper will not present the details on how to plan a feasible end-effector trajectory q r e (t) in W e , but just assumes that there exists such a smooth feasible end-effector trajectory q r e (t) over the predefined time interval [0, T f ] such that q r e (0) = q init e and q r e (T f ) = q fin e . The objective of this paper is then to design a robust controller u such that q e (t) converges to q r e (t) in the presence of external disturbances.

IV. NONLINEAR CONTROLLER DESIGN

As is stated above, the control problem of the robot's end-effector is actually a problem of output regulation, as described in [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Therefore, we can re-formulated it into the following projected space:

qe + ↵(q, q) = CM 1 H(q) T u (5)
with

↵(q, q) = CM 1 D( q, q) q + CM 1 K(q)q
With the planned smooth feasible end-effector trajectory q r e (t) over [0, T ], let us define the following tracking error:

e(t) = q r e (t) q e (t) (6)
where q r e (t) stands for the given desired end-effector trajectory and q e (t) represents the measurable end-effector position.

With the introduced notations, ( 6) can be rewritten as:

ë = ↵ CM 1 H(q) T u + qr e (7)
In order to achieve the convergence of e to 0, the following nonlinear controller is proposed:

u = u + u o ( 8
)
where u o is the initial control input for the soft robot at the rest configuration and

u =(CM 1 H(q) T ) 1 ✓ qr e (t) + Z t 0 (K 1 ë(s) +K 2 ė(s) + K 3 e(s) + K 4 Z s 0 e(⌧ )d⌧ )ds ◆ ( 9
)
where K 1 , K 2 , K 3 and K 4 are the controller gains to be designed. Inserting the above nonlinear controller into [START_REF] Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF], we can finally obtain:

ë(t) =↵ Z t 0 (K 1 ë(s) + K 2 ė(s) + K 3 e(s) + K 4 Z s 0 e(⌧ )d⌧ )ds CM 1 H T (q)u o
Now, let us note

↵ ⇤ (t) = ↵ CM 1 H T uo = CM 1 (q) ⇥ D( q, q) q + K(q)q H T (q)uo ⇤ ( 10
)
and regard ↵ ⇤ (t) as a bounded disturbance. Then, the above equation becomes

ë + Z t 0 ✓ K 1 ë + K 2 ė + K 3 e + K 4 Z s 0 e(⌧ )d⌧ ◆ ds = ↵ ⇤ (t) (11
) which is equivalent to: e (4) + K 1 e (3) + K 2 e (2) + K 3 e (1) + K 4 e = ↵⇤ (t) [START_REF] Gravagne | Uniform regulation of a multisection continuum manipulator[END_REF] where e (i) means the i-th derivative of e with respect to time. It is worth noting that ( 12) is a classical 4-order linear system with a bounded disturbance ↵⇤ (t), therefore we can always find matrices K 1 , K 2 , K 3 and K 4 to guarantee that the tracking error e(t) will converge to a neighborhood of 0 whose size depends on the magnitude of ↵⇤ (t) [START_REF] Hespanha | Linear systems theory[END_REF].

The above analysis shows that the performance of the proposed controller (8) depends on the characteristics of ↵ ⇤ (t). It is clear that if ↵ ⇤ (t) contains high frequency, then the magnitude of ↵⇤ (t) will become larger than that of ↵ ⇤ (t). Due to the same reason, the magnitude of ↵⇤ (t) might be smaller if the highest frequency of ↵ ⇤ (t) is less than 1 2⇡ . Moreover, the smaller ↵⇤ (t) is, the closer e converges to 0, i.e., q e (t) will converge to a smaller neighborhood of q r e (t). Consequently, in order to get a satisfactory performance for the proposed controller [START_REF] Hannan | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], the planned trajectory should not vary too quickly. In practice, such a constraint can be easily satisfied by choosing a relatively larger value of T f , i.e., the time interval of the planned trajectory.

Remark 2: When considering nonlinear effects and frictions discussed in Remark 1, according to [START_REF] Gravagne | Large deflection dynamics and control for planar continuum robots[END_REF], it is clear that those issues are included into ↵ ⇤ (t). Therefore, the controller proposed in this paper is robust even if we have no precise model of those nonlinearities. The same assertion is true when the measurement is corrupted with noise. To show this, let us model the output of sensor as q e = Cq + ✏(t) where ✏(t) represents the noise in measurement channel. Following the similar procedure, we will end up with:

e (4) + K 1 e (3) + K 2 e (2) + K 3 e (1) + K 4 e = ↵⇤ (t) ✏(t)
Again, by regarding ↵ = ↵ ⇤ (t) ✏(t), we can draw the same robustness conclusion, i.e., if ↵ is bounded, then e(t) will converge to the neighborhood of 0.

The above discusses how to design the nonlinear controller for the trunk robot, the next section presents several experimental tests regarding trajectory tracking. The overall experimental setup of the testing platform is shown in Fig. 3, mainly including the mechanical setup, the Polhemus magnetic position sensor (which is placed at the end effector through the long internal hole of the studied robot as described in Fig. 1) and the actuation system.

In our work, a front circle trajectory (Fig. 5a) and a side circle trajectory (Fig. 5b) are selected for validation where a common PID controller and our proposed nonlinear controller are compared. Besides, the robustness (against temporary and permanent disturbances) of the proposed controller has also been experimentally proved. In the current literature, the experimental validation for most trunk-like robots is quite limited to front trajectories whereas we have added the side circle trajectory (Fig. 5b) in a highly nonlinear zone to further prove the validity of our proposed controller.

A. Trajectory tracking tests 1) Trajectory tracking using a common PID controller: In this test, the trunk robot was expected to track two different circle trajectories in X-Y plane and Y-Z plane respectively (see Fig. 5a and Fig. 5b) using a common PID controller.

It is well-known that tuning criteria for PID exists for SISO (single-input-single-output) systems. However, it is generally difficult, and there does not exist universal methods to tune those paramters when dealing with MIMO (multiinput-multi-output) systems. In this paper, we adopted our former work [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] to tune PID controller.

The results are demonstrated in Fig. 4. As shown in Fig. 4a, the end-effector tended to track the X-Y plane circle trajectory but was not able to do so. This is because the common PID cannot effectively deal the hysteresis when the tracked trajectory is time-varying (in this case, the average speed of the trajectory is around 1.5 cm/s). Fig. 4b indicates that the robot was out of control when tracking the Y-Z plane circle trajectory, and the common PID is invalid in this area of the workspace.

2) Trajectory tracking using the proposed nonlinear controller: In this test, we used the proposed nonlinear controller to control the robot, and the trunk end-effector was expected to track the two mentioned circle trajectories (Fig. 5). The corresponding results are displayed in Fig. 6a and Fig. 6b. It can be easily noticed that the proposed nonlinear controller enabled the robot's end-effector to track both X-Y plane circle trajectory and Y-Z plane circle trajectory with almost no hysteresis. This means the proposed nonlinear controller can more rapidly react to real-time trajectory error, presenting dominant advantages over the common PID controller.

B. Robustness Analysis

To test the robustness of the proposed controller, two more trajectory tracking tests were conducted under two different scenarios (temporary disturbance and permanent disturbance) shown in Fig. 7. 5a) under temporary disturbances exerted by a pen (Fig. 7a) . During the testing, the robot was interfered by 7 temporary disturbances when tracking the trajectory. According to the results in Fig. 8, the end-effector quickly reacted to the disturbances and caught up with the trajectory, which presents great robustness of the proposed controller against temporary disturbances.

VI. CONCLUSIONS

Soft robotics has attracted much research attention, presenting great potential as well in a wide variety of practical applications due to its unique characteristics compared to conventional rigid robots. Nevertheless, we still need to face many tough challenges regarding dynamically controlling them nowadays. In this paper, we propose a nonlinear controller for a soft trunk robot to track user-defined global trajectories. The feasibility has been proved via 4 different tests, and the robustness of the proposed method has been illustrated by rejecting the temporary and permanent external disturbances. However, as stated previously, the proposed controller also has some limitations, for example, the planned trajectory should not vary too rapidly. Therefore, in our future work, we will propose to employ a new control strategy to deal with control objectives regarding faster trajectory tracking.

  Young's modulus: E = 6.5 ⇥ 10 6 Pa Poisson ratio: ⌫ = 0.45 Mass: m = 0.045 kg

Fig. 3 :Fig. 4 :

 34 Fig. 3: Overall experimental setup

Fig. 5 :

 5 Fig. 5: Circle trajectory tracking

Fig. 6 :

 6 Fig. 6: Trajectory tracking using the proposed nonlinear controller

  (a) Temporary disturbance by touching the robot with a pen (b) Permanent disturbance by adding an extra piece of rubber on the robot.

Fig. 7 :

 7 Fig. 7: External disturbance tests

Fig. 8 :Fig. 9 :

 89 Fig. 8: Temporary disturbance results

2) Permanent disturbance: In this test, the trunk robot was expected to track the X-Y circle trajectory while carrying an extra piece of 0.033 kg rubber (almost the same weight as the trunk robot, see Fig.7b). According to the results demonstrated in Fig.9, the robot quickly caught up with the trajectory at the beginning of the trajectory, and then stably tracked the circle trajectory, which therefore proved the robustness of the controller against permanent disturbances.

*This work was supported by INRIA