Modeling and Optimal Control of an Octopus Tentacle
MMODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE
SIMONE CACACE, ANNA CHIARA LAI, AND PAOLA LORETI
Abstract.
We present a control model for an octopus tentacle, based on the dynamics of an inextensiblestring with curvature constraints and curvature controls. We derive the equations of motion together with anappropriate set of boundary conditions, and we characterize the corresponding equilibria. The model resultsin a system of fourth-order evolutive nonlinear controlled PDEs, generalizing the classic Euler’s dynamicelastica equation, that we approximate and solve numerically introducing a consistent finite differencescheme. We proceed investigating a reachability optimal control problem associated to our tentacle model.We first focus on the stationary case, by establishing a relation with the celebrated Dubins’ car problem.Moreover, we propose an augmented Lagrangian method for its numerical solution. Finally, we address theevolutive case obtaining first order optimality conditions, then we numerically solve the optimality systemby means of an adjoint-based gradient descent method. Introduction
Octopus tentacles are extremely complex and fashinating biological structures, whose study is motivatedby both natural science interests as well as applications to robotics, in particular, in the framework of softmanipulators [25]. In this paper we are interested in an optimal control problem associated to the dynamicsof an octopus tentacle. The aim is to model the voluntary muscle contractions of the tentacle and to deriveoptimality conditions for a reachability problem. To this end, we propose a continuous model, based on thedynamics of an inextensible string with curvature constraints and curvature controls. More precisely, ourmodel in an extension of the
Euler’s dynamic elastica equation (see equation (4) below) in the followingsense: we add to the exact inextensibility constraint and the bending momentum, a term representing acurvature constraint and a control term. The curvature constraint prevents the tentacle from bending over afixed treshold; while the control term, which we consider the main novelty of this model, forces pointwise thecurvature of the tentacle. After characterizing explicitly the equilibria of the tentacle system as functions ofthe controls, both in terms of shape configuration and tension, we introduce a finite difference scheme forits numerical approximation, and we present some simulations validating the results. Then, we address thefollowing optimal control problem: steer the tentacle tip to a target point minimizing the activation of thetentacle muscles. More precisely, we want to optimize the distance of the tentacle tip from a fixed targetpoint and a quadratic cost associated to the controls. This is done in two different setting, a stationarycase and a dynamic case. In the stationary case, we obtain a characterization of the reachable set for thetentacle tip, by establishing a relation with a classical motion planning problem first posed by Markov in1889 [21], also known as
Dubins car problem . Moreover, we provide an augmented Lagrangian methodfor the numerical solution of the optimal control problem and present some simulations. In the dynamiccase, the objective functional accounts for the whole evolution of the tentacle on a time interval [0 , T ], plusa term representing the kinetic energy at the final time T . The aim is then also to stop the tentacle assoon as possible. To this end we introduce an adjoint state and derive first order optimality conditions.The resulting optimality system is then solved numerically via an adjoint-based gradient descent methodand some simulations are presented. We remark that, to the best of our knowledge, the present theoreticapproach, in particular the dynamic case based on the optimal control of PDEs (see e.g. [27]), appears tobe new.We now outline the main features of our model, while postponing a more formal justification to Section 2.The unknowns are a curve q ( s, t ) ∈ R describing the shape of the tentacle, and the associated inextensibility Key words and phrases.
Octopus tentacle modeling, optimal control, inextensible elastic rods, soft manipulators. a r X i v : . [ m a t h . O C ] N ov S. CACACE, A. C. LAI, AND P. LORETI multiplier σ ( s, t ) ∈ R , playing the role of the tension of the tentacle. We remark that the variable s ∈ [0 , q , and we denote by q s , q ss , q tt partial derivatives inspace and time respectively. The tentacle model includes:- the inextensibility constraint, given by the equation | q s | = 1 and encompassed in the dynamics viaa scalar Lagrange multiplier σ ;- the bending momentum, i.e. an intrinsic resistance of the tentacle to bending. It is described byan elastic potential of the form ε ( s ) | q ss | , where the function ε represents a non uniform bendingstiffness;- a curvature constraint, given by the inequality | q ss | ≤ ω ( s ), for some threshold positive function ω . It is imposed using an unilateral potential of the form ν ( s ) (cid:0) | q ss | − ω ( s ) (cid:1) , where ν is a nonuniform elastic constant and ( · ) + denotes the positive part of its argument;- a control term forcing the curvature of the tentacle, given by the equality constraint q s × q ss = ω ( s ) u ( s, t ), for a control map u valued in [ − , µ ( s ) ( ω ( s ) u ( s, t ) − q s × q ss ) , where µ is the corresponding non uniform elastic constant and q s × q ss represents the signed curvature (see Section 2 below for the definition of the vector product × in R ).The tentacle dynamics is then described by the following system of fourth order, evolutive, non-linear,controlled PDEs: for ( s, t ) ∈ (0 , × (0 , T )(1) (cid:26) ρq tt = (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss | q s | = 1where G [ q, ν, ε, ω ]( s, t ) := ε ( s ) + ν ( s ) (cid:0) | q ss ( s, t ) | − ω ( s ) (cid:1) + ,H [ q, µ, u, ω ]( s, t ) := µ ( s ) ( ω ( s ) u ( s, t ) − q s ( s, t ) × q ss ( s, t )) , the function ρ ( s ) represents the non uniform mass distribution of the tentacle, and the symbol v ⊥ denotesthe counterclockwise orthogonal vector to v . The system is completed with appropriate boundary and initialconditions.To put the tentacle model into perspective, we now briefly discuss its relation with some well-known equa-tions in the literature, depending on specific choices of the parameters, summarized in the following table:Inextensibility Bending Curvature Curvaturemomentum constraints controlsEuler-Bernoulli beam × (cid:88) × × Whip equation (cid:88) × × ×
Dynamic elastica (cid:88) (cid:88) × ×
Tentacle equation (1) (cid:88) (cid:88) (cid:88) (cid:88)
Euler-Bernoulli beam . This model describes the vibration of a beam, depending on the given boundaryconditions. A typical choice is the cantilevered case, in which one endpoint of the beam is clamped andthe other one is free. To retrieve it, we neglect the inextensibility constraint, we set ν ( s ) = µ ( s ) ≡ ε ( s ) = const , so that (1) reduces to(2) ρq tt = − εq ssss . Assuming analytical initial data and uniform mass distribution ρ , Equation (2) admits a unique classical(vector) solution given by q ( s, t ) = (cid:88) k ∈ N a k e − iϕ k t (cosh( β k s ) − cos( β k s ) + γ ( β k )(sinh( β k s ) − sin( β k s ))) + sq s ( t, , where the coefficients a k ∈ C depend only on the initial data, the sequence of real numbers { β k } is composedby the positive solutions of the equation cosh( β ) cos( β ) + 1 = 0 , ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 3 while ϕ k and γ ( β k ) are given, for each k , respectively by ϕ k = β k (cid:114) ερ , γ ( β k ) = cosh β k + cos β k sinh β k + sin β k . Whip equation . Assuming ν ( s ) = µ ( s ) = ε ( s ) ≡
0, Equation (1) reduces to a second-order differentialequation modeling the motion of an inextensible string(3) (cid:26) ρq tt = ( σq s ) s | q s | = 1Typical boundary conditions prescribe one fixed endpoint and zero tension at the other endpoint. Theproblem results in a nonlinear system with a boundary layer at the free end, where the motion degeneratesinto a first order equation in space. To the best of our knowledge, one of the most complete results for thissystem is given in [24], proving existence and uniqueness of solutions in suitable weighted Sobolev spaces. Euler’s dynamic elastica . Assuming ν ( s ) = µ ( s ) ≡ ε ( s ) = const , we obtain the following equation,also known as the nonlinear Euler-Bernoulli beam equation(4) (cid:26) ρq tt = ( σq s ) s − εq ssss | q s | = 1Classical boundary conditions prescribe the position and the tangent of both endpoints. In this case, theproblem is well-posed and admits a smooth solution, representing the curve of fixed length that connectsthe two endpoints minimizing the integral of the squared curvature.The problem of modeling and control soft inextensible bodies, including octopus tentacles, has beenattacked with several approaches, including tools coming from calculus of variations, control theory, motionplanning for robotics. This resulted in a wide scientific literature which is almost impossible to exhaustivelysummarize. We then conclude this introduction by reporting only the papers that mostly inspired our work.For general references on the modeling of octopus tentacles we mention the series of papers [29, 30], wherea rich discrete control model is proposed: note that the main difference with the present paper is that wediscuss a simpler but fully continuous model. For the modeling framework, we refer to the survey [28] for anintroduction on the derivation of equations for inextensible strings (with no curvature constraints) and theirdiscrete counterpart, the chains. The papers [8, 9] treat the same topic, stressing the numerical aspects ofthe problem. The paper [22] provides a nice investigation of boundary conditions for the Euler’s dynamicelastica and presents some numerical simulations — also, see the paper [2] for a discussion of the discretecase. The papers [16, 17] contain a discussion in a control and number theoretic setting of the discrete,stationary case with an exponential decay in the mass distribution. The paper [19] includes more general,self-similar mass distributions; we finally refer to [18] for a first investigation of the continuous counterpartof these models.The problem related to curvature constraints and controls is widely investigated in the framework of softrobotics: the problem is in general to prescribe the curvature of a flexible device in a distributed fashion,in order to perform locomotion or grasping, see [26, 20]. Among many others, we refer to [12] for the studyof the kinematics of a soft-robot prototype. In [14, 13], after deriving the dynamics of a discrete tentacle-like actuator, a control strategy is prescribed a priori, and its efficiency is either numerically confirmed ormeasured by actual prototypes, respectively. An interesting approach based on machine learning can befound in [6].One our main results consists in showing that equilibria of (1) solve an ordinary differential equation re-curring in motion planning problems. For an overview on the aspects that are more related to our approach,we refer to the paper [1] and the references therein. In order to characterize the reachable set of the equilibria,we rely on the results in [4], see also [10]. Finally, we remark that the stationary case studied in this paperis very close to the framework of Euler’s elastica , i.e., the study of static equations for elastic rods – see [15]for an hystorical excursus on the problem, and the papers [7, 3, 5] and the references therein for an overview.The paper is organized as follows. Section 2 is devoted to our octopus tentacle model, including thederivation and the numerical approximation of its equations of motion, and also the characterization its
S. CACACE, A. C. LAI, AND P. LORETI equilibria. In Section 3, we present and numerically solve the optimal control problem in the stationarycase. Moreover, we characterize the reachable set of the tentacle tip. Finally, in Section 4, we address theoptimal control problem in the dynamic case. We derive first order optimality conditions, and we numericallysolve the corresponding optimality system. For each section, some numerical simulations complete thepresentation. 2.
A control model for an octopus tentacle
This section is devoted to our octopus tentacle model. Starting from a discrete particle system describinga non uniform N -pendulum, we encode all the geometrical constraints into a discrete Lagrangian for thesystem, then we obtain the continuous model as a formal limit for the number of particles N going toinfinity. We proceed by applying the least action principle to obtain the equation of motions, and we studyits equilibrium configurations depending on the given control map. Finally, we discretize the equations bymeans of a finite difference scheme and we present some simulations for validating the model.Let us fix some notations. The scalar product between two vectors q, w ∈ R is denoted by q · w , whereasthe Euclidean norm is denoted by | q | = √ q · q . Moreover, we define the counterclockwise orthogonal vector q ⊥ := (cid:18) − (cid:19) q and, with a little abuse of notation, we denote by q × w the vector product in R , namely the scalar quantity q × w := q · w ⊥ . We use the symbol ( · ) + to denote the positive part function: ( x ) + = x if x > x ) + = 0otherwise. We also use the symbol 0 to denote the null vector in R .2.1. From a discrete particle system to a continuous tentacle.
Our first step of modeling is toconsider the tentacle as a three-dimensional body with an axial symmetry, with a fixed endpoint and char-acterized by longitudinal inextensibility: the evolution of the curve on the plane representing the symmetryaxis is the object of our investigation. Roughly speaking, we neglect the torsion of the curve and reduce theproblem in a two-dimensional setting, approximating the whole tentacle with a planar inextensible string.We proceed by introducing a discretization, which is, as a matter of fact, the cornerstone of both the desiredcontinuous model and the numerical scheme for its simulation. We sample the shape of the tentacle by N particles (or joints) with positions q k ∈ R and masses m k >
0, for k = 1 , ..., N . We pose an additionalparticle q fixed at the origin, representing the anchor point of the tentacle. Moreover, for formal compu-tations, we also add two ghost particles q − and q N +1 at the endpoints: the first particle is aligned to thevertical axis, i.e. q − = q + (cid:96)e , where (cid:96) > N → ∞ , while e denote the unit vector (0 , q N +1 = q N + ( q N − q N − ). We introduce the dependency of the joints on time and,with a little abuse of notation, we denote the ( N + 3)-tuple of vectors ( q − ( t ) , q ( t ) , ..., q N +1 ( t )) by q ( t ), orsimply q depending on the context. Similarly, we denote the time derivative by ˙ q .We now impose the inextensibility constraint, assuming that the distance between consecutive joints isexactly (cid:96) , i.e., for k = 1 , ..., N | q k − q k − | = (cid:96) . We remark that, in most cephalopods species, a tentacle is thicker when close to the body and thinner whenreaching the free endpoint. This morphological aspect can be easily encoded by assuming that the sequenceof masses m k is decreasing in k . Note that, at this stage, the considered particle system can be viewed asjust a non uniform N -pendulum, or a chain.To deal with the ability of the tentacle to bend (as if it were a beam), we introduce some additionalconstraints involving its curvature. We assume that, for k = 1 , ..., N , the angle formed by the three consec-utive joints q k − , q k , q k +1 is bounded by a maximum angle α k . This relation can be easily expressed by thefollowing inequality constraint: ( q k +1 − q k ) · ( q k − q k − ) ≥ (cid:96) cos( α k ) . On the other hand, we take into account a bending momentum, namely we assume that the tentacle has anintrinsic resistance to bending. The position at rest corresponds to null relative angles, as described by the
ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 5 following equality constraint: ( q k +1 − q k ) × ( q k − q k − ) = 0 . Finally, we endow our model with a control term, pointwise prescribing the curvature of the tentacle. Moreprecisely, we can force the exact angle between the joints q k − , q k , q k +1 by means of the following equalityconstraint: ( q k +1 − q k ) × ( q k − q k − ) = (cid:96) sin( α k u k ) , where u k ∈ [ − ,
1] is the control associated to k -th joint. In this way, we take into account the fact that anoctopus tentacle is composed by independent muscular tissues. The local contraction of a muscle, which isthe way an octopus controls its limbs, prescribes indeed the local curvature.The next step is to form a suitable Lagrangian for our discrete model. To this end, we impose the threecurvature constraints discussed above via penalization, namely we define respectively the following quadraticpotentials: for k = 1 , ..., NG k ( q ) := ν k (cid:18) cos( α k ) − (cid:96) ( q k +1 − q k ) · ( q k − q k − ) (cid:19) ,B k ( q ) := ε k (cid:16) ( q k +1 − q k ) × ( q k − q k − ) (cid:17) ,H k ( q ) := µ k (cid:18) sin( α k u k ) − (cid:96) ( q k +1 − q k ) × ( q k − q k − ) (cid:19) , where the penalty parameters ν k , µ k , ε k play the role of elastic constants. This choice is motivated by thefact that the tentacles are soft elastic structures, allowing for small violations of their constraints, at leaston their radial components. We can incorporate in the elastic constants the morphological property thatthin joints bend easier than thick joints. In this respect, we can assume that the sequences ν k , µ k , ε k aredecreasing in k , whereas the sequence of angles α k is increasing in k .On the other hand, we impose the inextensibility constraint exactly, by defining for k = 1 , ..., NF k ( q ) := σ k (cid:0) | q k − q k − | − (cid:96) (cid:1) , where σ k is a Langrange multiplier representing the tension exerted along the tentacle. We must saythat, actually, longitudinal extension of the tentacles is also observed in many species (e.g. in the octopusvulgaris ), so one may wonder about the choice of an exact constraint. This is essentially motivated by threeaspects. First, we recall that the muscolar system of an octopus is formed by volume preserving modulescalled hydrostats : when a hydrostat is contracted, it shortens. The variability of the length of the tentacle ismainly due to this compensation mechanism, rather than mere elastic forces. We guess that a more accuratemodel should include a curvature depending local inextensibility constraint. Therefore, at this stage of ourinvestigation, we opted for a constant exact constraint which may be a first step towards this more accuratemathematical model. The second reason is given by the fact that octopus models are massively applied tosoft robotics, and a large class of these devices are indeed inextensible. Finally, as shown below, this choiceallows to slightly simplify the expression of the continuous Lagrangian of the system and, consequently, thecorresponding equations of motion.Now we are ready to introduce the Lagrangian associated to the considered particle system:(5) L N ( q, ˙ q, σ ) := N (cid:88) k =1 (cid:18) m k | ˙ q k | − (cid:96) F k ( q ) − (cid:96) G k ( q ) − (cid:96) B k ( q ) − (cid:96) H k ( q ) (cid:19) , where the first term represents the kinetic energy, while the other terms account for the constraints, suitablyrescaled in view of the limit N → ∞ we are going to perform. To this end, we assume that the total lengthof the chain is equal to 1, and we set (cid:96) = 1 /N . Moreover, we take the masses and the maximum bendingangles of the form m k = (cid:96)ρ k and α k = (cid:96)ω k respectively, for given mass distributions ρ k and curvatures ω k . Finally, we consider smooth functions ν, µ, ε, ρ, ω : [0 , → R + and, for T >
0, smooth functions q : [0 , × [0 , T ] → R , σ : [0 , × [0 , T ] → R and u : [0 , × [0 , T ] → [ − ,
1] such that, for all N , k = 1 , ..., N and t ∈ [0 , T ] ν k = ν ( k(cid:96) ) , µ k = µ ( k(cid:96) ) , ε k = ε ( k(cid:96) ) , ρ k = ρ ( k(cid:96) ) , ω k = ω ( k(cid:96) ) , S. CACACE, A. C. LAI, AND P. LORETI q k ( t ) = q ( k(cid:96), t ) , σ k ( t ) = σ ( k(cid:96), t ) , u k ( t ) = u ( k(cid:96), t ) . We remark that the space variable, denoted by s ∈ [0 , t by q ( s, t ). We employ standard notations q s , q ss for partial derivatives in space, and q t , q tt forpartial derivatives in time. Moreover, we set q ks := q s ( (cid:96)k ), q kss := q ss ( (cid:96)k ) and consider following Taylorexpansions: q k − q k − = (cid:96)q ks + o ( (cid:96) ) = (cid:96)q ks − (cid:96) q kss + o ( (cid:96) ) ,q k +1 − q k = (cid:96)q ks + 12 (cid:96) q kss + o ( (cid:96) ) ,q k +1 − q k + q k − = (cid:96) q kss + o ( (cid:96) ) . Then, as N → ∞ (hence (cid:96) → (cid:96) F k ( q ) = 1 (cid:96) σ k (cid:0) | q k − q k − | − (cid:96) (cid:1) = 1 (cid:96) σ k (cid:0) (cid:96) | q ks | + o ( (cid:96) ) − (cid:96) (cid:1) = (cid:96)σ k ( | q ks | −
1) + o ( (cid:96) ) . Since the inextensibility constraint is exact, we can use the relation | q k − q k − | = (cid:96) , for all k = 1 , ..., N , tosimplify the curvature potential given by G k . Indeed, we get1 (cid:96) G k ( q ) = 1 (cid:96) ν k (cid:18) cos( (cid:96)ω k ) − (cid:96) ( q k +1 − q k ) · ( q k − q k − ) (cid:19) = 1 (cid:96) ν k (cid:18) − (cid:96) ω k + o ( (cid:96) ) + 12 (cid:96) (cid:0) | q k +1 − q k + q k − | (cid:1) − (cid:19) = 1 (cid:96) ν k (cid:18) − (cid:96) ω k + o ( (cid:96) ) + 12 (cid:96) (cid:0) (cid:96) | q kss | + o ( (cid:96) ) (cid:1)(cid:19) = 14 (cid:96)ν k (cid:0) | q kss | − ω k (cid:1) + o ( (cid:96) ) . Moreover, recalling that q ks × q ks = 0, q kss × q kss = 0 and q ks × q kss = − q kss × q ks by definition of the vectorproduct, we obtain 1 (cid:96) B k ( q ) = 1 (cid:96) ε k (cid:16) ( q k +1 − q k ) × ( q k − q k − ) (cid:17) = 1 (cid:96) ε k (cid:16) ( (cid:96)q ks + (cid:96) q kss + o ( (cid:96) )) × ( (cid:96)q ks − (cid:96) q kss + o ( (cid:96) )) (cid:17) = 1 (cid:96) ε k (cid:16) (cid:96) q kss × q ks − (cid:96) q ks × q kss + o ( (cid:96) ) (cid:17) = (cid:96)ε k (cid:0) q ks × q kss (cid:1) + o ( (cid:96) ) . Similarly, 1 (cid:96) H k ( q ) = 1 (cid:96) µ k (cid:18) sin( (cid:96)ω k u k ) − (cid:96) ( q k +1 − q k ) × ( q k − q k − ) (cid:19) = 1 (cid:96) µ k (cid:16) (cid:96)ω k u k − (cid:96)q ks × q kss + o ( (cid:96) ) (cid:17) = (cid:96)µ k (cid:0) ω k u k − q ks × q kss (cid:1) + o ( (cid:96) ) . Putting together all the terms, summing on k = 1 , ..., N and taking the formal limit as N → ∞ , we obtainthat L N in (5) converges to L ∞ ( q, q t , σ ) := (cid:90) (cid:18) ρ | q t | − σ ( | q s | − − ν (cid:0) | q ss | − ω (cid:1) − ε ( q s × q ss ) − µ ( ωu − q s × q ss ) (cid:19) ds . Note that | q s | = 1 is the continuous counterpart of the inextensibility constraint. Again, since this constraintis imposed exactly via the Lagrange multiplier σ , we can employ it to simplify the bending momentum term ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 7 in L ∞ . Indeed, differentiating | q s | = 1 with respect to s , we get q s · q ss = 0, namely q s and q ss areorthogonal. This implies that ( q s × q ss ) = | q s | | q ss | = | q ss | , and we can finally obtain the followingequivalent Lagrangian for the continuous tentacle model: L ( q, q t , σ ) := (cid:90) (cid:18) ρ | q t | − σ ( | q s | − − ν (cid:0) | q ss | − ω (cid:1) − ε | q ss | − µ ( ωu − q s × q ss ) (cid:19) ds . (6)2.2. Equations of motion.
Here we derive the equations of motion for the continuous octopus tentaclemodel, by applying the classical least action principle to the Lagrangian (6). To this end, we first assigninitial data, i.e. we impose to the tentacle initial position and velocity, by means of smooth functions q , v : [0 , → R :(7) q ( s,
0) = q ( s ) , q t ( s,
0) = v ( s ) , for s ∈ [0 , . Moreover, we impose two boundary conditions at one endpoint of the tentacle. Following the discrete model,we defined the first two particles of the approximating chain as q − = q + (cid:96)e and q = 0. In the firstrelation we recognize a backward finite difference discretization of the derivative q s (0), while the secondrelation simply prescribes the anchor point of the tentacle. Hence, we set(8) q (0 , t ) = 0 , q s (0 , t ) = − e , for t ∈ [0 , T ] . We now define the action(9) S ( q, σ ) = (cid:90) T L ( q, q t , σ ) dt . Note that the unknowns of the problem are both the curve q and the tension σ . According to the leastaction principle, for T > , T ] and between the states( q ( s, , σ ( s, q ( s, T ) , σ ( s, T )) is the one for which S ( q, σ ) is stationary, namely we require (cid:104) δ q S ( q, σ ) , w (cid:105) := lim α → α ( S ( q + αw, σ ) − S ( q, σ )) = 0 , (cid:104) δ σ S ( q, σ ) , χ (cid:105) := lim α → α ( S ( q, σ + αχ ) − S ( q, σ )) = 0 , for all admissible test functions w : [0 , × [0 , T ] → R and χ : [0 , × [0 , T ] → R . We remark that the notionof admissible test is related to the boundary conditions in space and time. More precisely, the pair ( w, χ ) ischosen so that the corresponding variation of ( q, σ ) keep the same boundary conditions. In particular, for s ∈ (0 , w ( s,
0) = w ( s, T ) = 0, χ ( s,
0) = χ ( s, T ) = 0, and, in view of (7), also w t ( s,
0) = 0.Moreover, in view of (8), we impose w (0 , t ) = w s (0 , t ) = 0 for t ∈ (0 , T ). Note that no condition on ( w, χ )is required for s = 1, since it corresponds to the free end of the tentacle. Nevertheless, as we show below,additional boundary conditions will appear at s = 1 for both q and σ , as a consequence of the stationarityof the action.To proceed, we need additional assumptions on some elastic constants appearing in the Lagrangian (6).For some ε >
0, we require(10) ε ( s ) ≥ ε , for s ∈ [0 , ,µ (1 , t ) = µ s (1 , t ) = 0 , for t ∈ [0 , T ] . The first assumption prevents a degeneracy of the equations of motion, in particular at the free end, since ε has a regularizing effect as discussed below. The second assumption is more technical, and allows one toconsiderably reduce the expression for the boundary conditions. In order to simplify the presentation, wealso define(11) G [ q, ν, ε, ω ]( s, t ) := ε ( s ) + ν ( s ) (cid:0) | q ss ( s, t ) | − ω ( s ) (cid:1) + ,H [ q, µ, u, ω ]( s, t ) := µ ( s ) ( ω ( s ) u ( s, t ) − q s ( s, t ) × q ss ( s, t )) . S. CACACE, A. C. LAI, AND P. LORETI
Let us start by computing the variation of S with respect to σ in direction χ , and impose stationarity:we immediately get (cid:104) δ σ S , χ (cid:105) = (cid:90) T (cid:90) ( | q s | − χ ds dt = 0 , for every admissible test χ . This implies the inextensibility constraint, i.e. | q s | = 1 almost everywhere in (0 , × (0 , T ) . We now compute the variation of S with respect to q in direction w : we have (cid:104) δ q S , w (cid:105) = (cid:90) T (cid:90) (cid:110) ρq t · w t − σq s · w s − Gq ss · w ss − H ( − w s × q ss − q s × w ss ) (cid:111) ds dt . Integrating by parts each term, we get (cid:90) T (cid:90) ρq t · w t ds dt = (cid:90) [ ρq t · w ] T ds − (cid:90) T (cid:90) ρq tt · w ds dt , − (cid:90) T (cid:90) σq s · w s ds dt = − (cid:90) T [ σq s · w ] dt + (cid:90) T (cid:90) ( σq s ) s · w ds dt , − (cid:90) T (cid:90) Gq ss · w ss ds dt = − (cid:90) T [ Gq ss · w s ] dt + (cid:90) T [( Gq ss ) s · w ] dt − (cid:90) T (cid:90) ( Gq ss ) ss · w ds dt and, recalling that q × w = q · w ⊥ = − q ⊥ · w , − (cid:90) T (cid:90) H ( − w s × q ss − q s × w ss ) ds dt = (cid:90) T (cid:90) ( Hq ⊥ ss · w s − Hq ⊥ s · w ss ) ds dt = (cid:90) T (cid:2) Hq ⊥ ss · w (cid:3) dt − (cid:90) T (cid:2) Hq ⊥ s · w s (cid:3) dt + (cid:90) T (cid:2)(cid:0) Hq ⊥ s (cid:1) s · w (cid:3) dt − (cid:90) T (cid:90) (cid:0) Hq ⊥ ss (cid:1) s · w ds dt − (cid:90) T (cid:90) (cid:0) Hq ⊥ s (cid:1) ss · w ds dt . Summing up, we obtain (cid:104) δ q S , w (cid:105) = (cid:90) T (cid:90) (cid:8) − ρq tt + ( σq s ) s − ( Gq ss ) ss − (cid:0) Hq ⊥ ss (cid:1) s − (cid:0) Hq ⊥ s (cid:1) ss (cid:9) · w ds dt + (cid:90) (cid:104) ρq t · w (cid:105) T ds + (cid:90) T (cid:104) (cid:8) − σq s + ( Gq ss ) s + Hq ⊥ ss + (cid:0) Hq ⊥ s (cid:1) s (cid:9) · w − (cid:8) Gq ss + Hq ⊥ s (cid:9) · w s (cid:105) dt . Using the admissibility conditions w ( s,
0) = w ( s, T ) = w (0 , t ) = w s (0 , t ) = 0 and imposing stationarity, wefinally get (cid:104) δ q S , w (cid:105) = (cid:90) T (cid:90) (cid:8) − ρq tt + ( σq s ) s − ( Gq ss ) ss − (cid:0) Hq ⊥ ss (cid:1) s − (cid:0) Hq ⊥ s (cid:1) ss (cid:9) · w ds dt + (cid:90) T (cid:16) (cid:8) − σq s + ( Gq ss ) s + Hq ⊥ ss + (cid:0) Hq ⊥ s (cid:1) s (cid:9) · w (1 , t ) − (cid:8) Gq ss + Hq ⊥ s (cid:9) · w s (1 , t ) (cid:17) dt = 0 . Now, choosing w with compact support in (0 , × (0 , T ), we have in particular w (1 , t ) = w s (1 , t ) = 0, sothat the first term in the previous expression should vanish. By the arbitrariness of w we obtain, almosteverywhere, the equations of motion ρq tt = (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss . On the other hand, choosing w such that w (1 , t ) = 0 and w s (1 , t ) is arbitrary, we get the following boundarycondition at s = 1: Gq ss + Hq ⊥ s = 0 . ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 9
Note that H vanishes due to the assumption µ (1 , t ) = 0. Moreover, G is strictly positive due to theassumption ε ( s ) ≥ ε >
0. Then, we obtain q ss (1 , t ) = 0 for t ∈ (0 , T ) . Similarly, choosing w such that w s (1 , t ) = 0 and w (1 , t ) is arbitrary, we also get − σq s + G s q ss + Gq sss + Hq ⊥ ss + H s q ⊥ s + Hq ⊥ ss = 0 , where we expanded, for convenience, the derivatives of the two products. Using again the assumption µ (1 , t ) = 0, and also µ s (1 , t ) = 0, we observe that both H and H s vanish at s = 1. Moreover, we alreadyknow that q ss (1 , t ) = 0, hence by definition of G we get(12) − σq s + εq sss = 0and, dot multiplying by q s , also(13) − σ | q s | + εq s · q sss = 0 . We now employ the inextensibility constraint | q s | = 1 obtained above. Differentiating twice with respectto s , we obtain q s · q ss = 0 and q s · q sss = −| q ss | respectively. Prolonging these relations by continuity upto s = 1, and plugging them into (13), we get − σ − ε | q ss | = 0 , which implies σ (1 , t ) = 0 for t ∈ (0 , T ), using again q ss (1 , t ) = 0. Coming back to (12), by the assumption ε ( s ) >
0, we conclude that q sss (1 , t ) = 0 for t ∈ (0 , T ) . Putting together all the initial and boundary conditions with the necessary conditions derived above, weend up with the following strong formulation of the system of motion for the tentacle model:(14) ρq tt = (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss in (0 , × (0 , T ) | q s | = 1 in (0 , × (0 , T ) q ( s,
0) = q ( s ) s ∈ (0 , q t ( s,
0) = v ( s ) s ∈ (0 , q (0 , t ) = 0 t ∈ (0 , T ) q s (0 , t ) = − e t ∈ (0 , T ) q ss (1 , t ) = 0 t ∈ (0 , T ) q sss (1 , t ) = 0 t ∈ (0 , T ) σ (1 , t ) = 0 t ∈ (0 , T )We remark that the relations at the free endpoint are known respectively as the zero bending momentum,zero shear stress and zero tension boundary conditions. They provide, with the anchor point and the fixedtangent at s = 0, the so called cantilevered boundary conditions for the classical nonlinear Euler-Bernoullibeam discussed in the Introduction. We also observe that the equation is of order four in space. The positiveterm ε in G plays the role of a regularization and prevents the equation to degenerate to the third order or,worse, to the second order if also H = 0 somewhere or, worse than worse, to the first order at s = 1 where σ = 0.To conclude this section, let us spend few words about dissipation for the tentacle model. We can introducetwo kinds of dissipation mechanisms: the first one is an environmental viscous friction proportional to thevelocity, modeling the fact that the tentacle can be immersed in a fluid; the second one is an internal viscousfriction, associated to the muscles contractions and hence proportional to the change in time of the tentaclecurvature. It is very well known that such dissipative forces do not fit the standard variational frameworkof Lagrangian mechanics, they can enter in the equations of motion only formally, introducing a suitableRayleigh dissipation functional of the form (cid:82) R ( q t ( s )) ds , and defining, ad hoc, the variation of the action(9) (note that the variation of R is performed with respect to q t and not q !): (cid:104) δ q S , w (cid:105) := (cid:90) T ( (cid:104) δ q L , w (cid:105) − (cid:104) δ q t R , w (cid:105) ) dt . Here we choose R ( q t ) = (cid:90) (cid:0) β ( s ) | q t | + γ ( s ) | q sst | (cid:1) ds , where β and γ are positive functions playing the role of viscous friction coefficients, typically proportionalto the local surface area. Since our tentacle model is thick at the base and thin at the tip, we can assume, asfor the bending parameters ν, ε, µ , that also β and γ are decreasing in s . Imposing stationarity of the actionabove, using integration by parts and the boundary conditions in (14), we easily obtain the correspondingfrictional forces in the right hand side of the equations of motion:(15) ρq tt = (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss − βq t − γq sssst . This dissipative version of the model will be employed in the following numerical simulations, in order toobserve the equilibrium configurations of the system.2.3.
Equilibria.
Here we characterize the equilibria of (14), namely we provide, up to an ordinary integra-tion in space, an explicit solution ( q, σ ) to the stationary problem(16) (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss = 0 in (0 , | q s | = 1 in (0 , q (0) = 0 q s (0) = − e q ss (1) = 0 q sss (1) = 0 σ (1) = 0in terms of a given stationary control map u : [0 , → [ − , Proposition 1.
Let u ∈ C [0 , and assume µ (1) = µ s (1) = 0 , ε ( s ) > . Then the stationary problem (16) admits a unique solution ( q, σ ) ∈ C ([0 , × C ([0 , , which is given, for s ∈ [0 , , by (17) q ( s ) = (cid:90) s (cid:16) sin( θ ( ξ )) , − cos( θ ( ξ )) (cid:17) dξ , σ ( s ) = ε ( s ) (cid:16) ¯ ω ( s ) u ( s ) (cid:17) , where θ ( ξ ) := (cid:90) ξ ¯ ω ( z ) u ( z ) dz and ¯ ω ( s ) := µ ( s ) ω ( s ) µ ( s ) + ε ( s ) . Proof.
First of all, we observe that a solution to (16) is, by construction, a minimum point of the potential(18) V ( q ) = (cid:90) (cid:18) ν (cid:0) | q ss | − ω (cid:1) + 12 ε | q ss | + 12 µ ( ωu − q s × q ss ) (cid:19) ds . Differentiating the inextensibility constraint | q s | = 1, we obtain the orthogonality condition q s · q ss = 0 andhence, for every s ∈ (0 , α ∈ R such that q ss = αq ⊥ s , since q s and q ⊥ s define a local base for R . Moreover, we have | q ss | = | α | and q s × q ss = α . It follows that the minimization of V in (18) can beperformed pointwise, namely minimizing, for every s ∈ (0 , α : f ( α ) = 14 ν (cid:0) α − ω (cid:1) + 12 εα + 12 µ ( ωu − α ) . We easily obtain that f (cid:48) ( α ) = αν (cid:0) α − ω (cid:1) + + εα − µ ( ωu − α ) = 0if and only if α = µωuµ + ε + ν ( α − ω ) + . Recalling that ε > u ∈ [ − , | α | < ω and then (cid:0) α − ω (cid:1) + = 0. We conclude that theunique minimum point of f is α = µωuµ + ε = ¯ ωu . ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 11
Hence, the stationary problem (16) is equivalent to the following second order problem: find q such that(19) q ss = ¯ ωuq ⊥ s in (0 , | q s | = 1 in (0 , q (0) = 0 q s (0) = − e Note that, due to the assumption µ (1) = µ s (1) = 0, the conditions q ss (1) = q sss (1) = 0, although irrelevant,are now embedded in the first equation of (19). We proceed by introducing polar coordinates such that q s ( s ) = (cid:16) sin( θ ( s )) , − cos( θ ( s )) (cid:17) , where θ : [0 , → R is unknown. It follows that | q s | = 1, while q s (0) = − e is equivalent to θ (0) = 0.Moreover, differentiating in s , we obtain q ss ( s ) = θ s ( s ) (cid:16) cos( θ ( s )) , sin( θ ( s )) (cid:17) = θ s ( s ) q ⊥ s ( s ) , then we end up with the following Cauchy problem in θ : (cid:26) θ s = ¯ ωu in (0 , θ (0) = 0The unique solution is given by θ ( s ) = (cid:90) s ¯ ω ( z ) u ( z ) dz , and, integrating q s , we conclude q ( s ) = (cid:90) s (cid:16) sin( θ ( ξ )) , − cos( θ ( ξ )) (cid:17) dξ , where, for s = 0, we also recover the condition on the anchor point q (0) = 0.We now recall the differential equation in the stationary problem (16), (cid:0) σq s − Hq ⊥ ss (cid:1) s − (cid:0) Gq ss + Hq ⊥ s (cid:1) ss = 0 . Integrating in s and expanding the second term, we get σq s − Hq ⊥ ss − H s q ⊥ s − Gq sss − G s q ss = C , for some constant vector C ∈ R . By continuity, for s →
1, using the boundary conditions σ (1) = H (1) = H s (1) = 0 and q ss (1) = q sss (1) = 0, we easily get C = 0. Moreover, dot multiplying by q s and using | q s | = 1, q s · q ss = 0 and q s · q sss = −| q ss | , we get σ − Hq s × q ss + G | q ss | = 0 . Finally, using the definitions of
G, H, ¯ ω and substituting q s × q ss = ¯ ωu , | q ss | = (¯ ωu ) , we conclude σ = 2 µ ( ωu − ¯ ωu )¯ ωu − ε (¯ ωu ) = 2 µ (cid:18) µ + εµ − (cid:19) (¯ ωu ) − ε (¯ ωu ) = ε (¯ ωu ) , and this completes the proof. (cid:3) Discretization.
In this section we give some details on the numerical approximation and solution ofthe system of motion (14), in particular we numerically validate the model and the characterization of itsequilibria provided by Proposition 1. We introduce a uniform discretization of the space-time [0 , × [0 , T ],namely we consider integers N, M and define grid nodes s k = k ∆ s for k = 0 , ..., N and t n = n ∆ t for n = 0 , ..., M , where the discretization steps are given respectively by ∆ s = 1 /N and ∆ t = T /M . For ageneric scalar or vector function χ defined on [0 , × [0 , T ] we adopt the standard notation χ nk ≈ χ ( s k , t n )to identify the corresponding approximation. In particular, we suitably define the approximations of all thedata functions ρ, ε, ν, ω, µ, u, q , v appearing in (14).Space discretization for the derivatives of the unknowns q, σ is performed by employing standard finitedifference operators. Denoting by D − , D + respectively first order backward and forward differences, and by D c second order central differences, we introduce the following approximation of the accelerations in theequations of motion: a ( q nk , σ nk ) := 1 ρ k (cid:0) D + (cid:0) σ nk D − q nk − H nk D c q n ⊥ k (cid:1) − D c (cid:0) G nk D c q nk + H nk D − q n ⊥ k (cid:1)(cid:1) with G nk = ε k + ν k (cid:0) | D c q nk | − ω k (cid:1) + , H nk = µ k (cid:0) ω k u nk − D − q nk × D c q nk (cid:1) . Note that the tangent q s is approximated by backward differences, whereas other first order derivatives areapproximated by forward differences. This choice clearly gives an overall approximation of order one in spacefor the equations of motion and naturally comes from the discrete model. More precisely, the approximations a ( q nk , σ nk ) coincide with the accelerations of the particle system discussed in Section 2.1: they can be derivedfrom the least action principle directly applied to the discrete Lagrangian L N in (5).Let us take a look at the boundary conditions. Introducing two ghost nodes s − = − ∆ s and s N +1 =( N + 1)∆ s at the ends of the interval [0 , q nN +1 − q nN + q nN − = 0 zero bending moment q nN +1 − q nN + 3 q nN − − q nN − = 0 zero shear stress σ nN = 0 zero tension q n = 0 anchor point q n − q n − = − e ∆ s fixed tangentWe remark that, from the first two (linear) equations, we readily get the values of q nN and q nN +1 in terms ofthe internal nodes, while from the last equation we obtain a fixed value for the node q n − .Finally, we proceed with the time integration, by employing a classical Velocity Verlet scheme. Thischoice is mainly motivated by the symplectic nature of the scheme, in particular its ability to preserveenergy regardless of the duration of the simulation. For k = 1 , ..., N − n = 1 , ..., M −
1, starting from q k = q ( s k ) and v k = v ( s k ), we discretize the equations of motion and the inextensibility constraints as(20) (cid:40) q n +1 k = q nk + v nk ∆ t + 12 a ( q nk , σ nk )∆ t | D − q n +1 k | = 1and we update the velocities by setting v n +1 k = v nk + 12 (cid:16) a ( q nk , σ nk ) + a ( q n +1 k , σ nk ) (cid:17) ∆ t . In the dissipative case, the frictional forces appearing in (15) can be included in a ( q nk , σ nk ) using the followingapproximations: − β ( s k ) q t ( s k , t n ) ≈ − β k v nk , − γ ( s k ) q sssst ( s k , t n ) ≈ − γ k D c D c v nk . Note that, for every time step n = 0 , ..., M −
1, the system (20) consists in 3( N −
1) equations in the unknowns( q n +1 k , σ nk ) k =1 ,...,N − . Due to the definition of a ( q nk , σ nk ), the first 2( N −
1) equations are linear, whereas thelast N − s = 2 × − (corresponding to N = 50 joints),∆ t = 10 − and ρ ( s ) = exp( − s ) , ε ( s ) = 10 − (1 − . s ) , ν ( s ) = 10 − (1 − . s ) , ω ( s ) = 2 π (1 + s ) ,µ ( s ) = (1 − s ) exp( − . s / (1 − s )) , β ( s ) = 4 − s, γ ( s ) = 10 − (4 − s ) . We observe that ε satisfies the assumption ε (1) = 10 − >
0, whereas µ is just a small perturbation of alinear function satisfying the assumption µ (1) = µ s (1) = 0. We choose the fully extended initial profile q ( s ) = − e s , the initial velocity v ( s ) ≡ u ( s, t ) ≡
1, corresponding to a fullcontraction of the tentacle, constant in time. Figure 1 shows the evolution computed by the numerical schemein the time interval [0 , ρ , togive a qualitative idea of the non uniform mass distribution. Moreover, we show the different configurationsin gray scales, from light gray to black as the time increases. ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 13 ≤ t ≤ .
35 0 . ≤ t ≤ . . ≤ t ≤ . . ≤ t ≤ Figure 1. tentacle dynamics for a full contraction controlThe effect of the curvature constraints is clearly visible, while the presence of frictional forces guaranteesconvergence to the equilibrium configuration given in Proposition 1. This is confirmed in Figure 2, showingthe asymptotic behavior of E q ( n ) = (cid:107)| D c q n |− ¯ ω (cid:107) ∞ and E σ ( n ) = (cid:107) σ n − ε ¯ ω (cid:107) ∞ , respectively the error betweenthe curvature at time t n = n ∆ t and the equilibrium curvature, and the error between the tension at time t n and the equilibrium tension. ( a ) ( b ) Figure 2.
Asymptotic behavior of curvature error E q (a) and tension error E σ (b)3. The stationary optimal control problem
In this section we address the static optimal control problem discussed in the Introduction. The goal isto touch a point with the tentacle tip and minimum effort. More precisely, given q ∗ ∈ R and τ >
0, weconsider the following optimization problem:(21) min (cid:26) (cid:90) u ds + 12 τ | q (1) − q ∗ | (cid:27) subject to q ss = ¯ ωuq ⊥ s in (0 , | q s | = 1 in (0 , | u | ≤ , q (0) = 0 q s (0) = − e where the first term in the functional quantifies the activation of the tentacle muscles in terms of the control,while the second term penalizes the distance of the tentacle tip from the target point q ∗ . We observe that,by virtue of Proposition 1, we have an explicit characterization of the constraints appearing in (21) in terms of the control. In the context of optimal control, formula (17) provides the so called control to state map,and it can be employed to completely remove q from the optimization problem above, yielding the followingoptimization in the control only:min | u |≤ (cid:40) (cid:90) u ds + 12 τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) (cid:18) sin (cid:18)(cid:90) s ¯ ω ( ξ ) u ( ξ ) dξ (cid:19) , − cos (cid:18)(cid:90) s ¯ ω ( ξ ) u ( ξ ) dξ (cid:19)(cid:19) ds − q ∗ (cid:12)(cid:12)(cid:12)(cid:12) (cid:41) . Note that the inextensibility constraint is now hidden in the integral formulation but, from a numerical pointof view, we found that this problem is quite involved, mainly due to the fact that, in polar coordinates,the tentacle tip q (1) is obtained by integration on the whole curve. This results in a non local dependencyof the functional on the control, so that the corresponding Euler-Lagrange equation is in turn an integralequation in u . Then, we decided to follow the opposite approach, namely remove the dependency on thecontrol, using the relation | q ss | = ¯ ω | u | . This gives the following equivalent formulation of the problem:min (cid:26) (cid:90) ω | q ss | ds + 12 τ | q (1) − q ∗ | (cid:27) subject to | q s | = 1 in (0 , | q ss | ≤ ¯ ω in (0 , q (0) = 0 q s (0) = − e (22)The constraints are still there, and they should be enforced by suitable multipliers. Nevertheless, the problemis more tractable, as we will show below introducing an augmented Lagrangian method for its numericalsolution. But first, we present an investigation on the reachability of target points q ∗ , i.e., we study the setof q ∗ for which the equation q (1) = q ∗ admits a solution, under the inextensibility and curvature constraints.In particular, we establish a connection with the celebrated Dubins car problem, a classical example studiedin motion planning [21].3.1.
Reachability.
Let us start by making an important assumption on the model, i.e., that the elasticconstants µ, ε and the curvature bound ω are balanced in the following way:¯ ω ( s ) = µ ( s ) ω ( s ) µ ( s ) + ε ( s ) ≡ ¯ ω for s ∈ [0 ,
1) and lim s → ¯ ω ( s ) = ¯ ω , for some ¯ ω >
0. Incidentally, notice that this implies that ω is unbounded near 1, since µ (1) = 0. Withthis assumption standing, we recover the same results of Proposition 1, by restricting our set of admissiblecontrols to those satisfying u (1) = u s (1) = 0. More precisely we define the set of admissible controls asfollows A := { u : [0 , → R | u ∈ C ([0 , , | u | ≤ , u (1) = u s (1) = 0 } and we rewrite the problem (19) in the following way:(23) q s = (sin( θ ) , − cos( θ )) in (0 , θ s = ¯ ω u in (0 , q (0) = 0 θ (0) = 0which is exactly the classical Dubins car control system [21]. Then, we define the reachable set as R := { q (1) | q = q [ u ] is a C solution of (23) , u ∈ A} . Moreover, we consider the wider control set A (cid:48) := { u : [0 , → R | u piecewise continuous , | u | ≤ } and let R (cid:48) := { q [ u ] | q [ u ] is a C solution of (23) , u ∈ A (cid:48) } . In [4] the set R (cid:48) is characterized via the following control set: A (cid:48)(cid:48) := { u : [0 , → { , ± } | u has a finite number of discontinuities } . ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 15
In order to properly cite this result, we adopt the following notations: q ∈ CL (respectively q ∈ CC ) if q = q [ u ] is a C solution of (23) with(24) u ( s ) = (cid:40) ¯ u s ∈ [0 , ¯ s ]0 (respectively, − ¯ u ) s ∈ (¯ s, s ∈ [0 ,
1] and ¯ u ∈ {± } . In other words, q [ u ] ∈ CL if it is composed by a circular arc with radius1 / ¯ ω and by a line segment, while q [ u ] ∈ CC if it is composed by two tangent circular arcs with radius 1 / ¯ ω .We are now in position to recall Cockayne and Hall’s result: Proposition 2. [4]
For every ¯ ω > , if q ∗ ∈ R (cid:48) , then there exists u ∈ A (cid:48)(cid:48) such that the Carath´eodorysolution q [ u ] of (23) satisfies q [ u ](1) = q ∗ . Moreover, R (cid:48) is a compact set and q [ u ](1) belongs to its boundaryonly if u is a piecewise constant function of the form (24) , i.e., ∂R (cid:48) ⊆ { q (1) | q ∈ CC ∪ CL } . Using this result, we can prove the following proposition.
Proposition 3.
For every ¯ ω > R = int ( R (cid:48) ) , where int ( · ) denotes the interior of a set. In particular, ∂R ⊆ { q (1) | q ∈ CC ∪ CL } . Proof.
First of all, we remark that, by definition, R ⊂ R (cid:48) . By Proposition 2, the boundary of R (cid:48) is reachableonly by piecewise constant controls, then we have R ⊆ int( R (cid:48) ). To prove the other inclusion, let q ∗ ∈ int( R (cid:48) ):again by Proposition 2 there exists a control u ∈ A (cid:48)(cid:48) be such that q [ u ](1) = q ∗ . Now, note that A is densein A (cid:48)(cid:48) with respect to the L norm, we then may consider an approximating control sequence u n ∈ A suchthat (cid:90) | u n − u | ≤ n ∀ n ∈ N . Setting θ n ( s ) = (cid:82) s ¯ ω u n ( ξ ) dξ and θ ( s ) = (cid:82) s ¯ ω u ( ξ ) dξ we also get the estimate | θ n ( s ) − θ ( s ) | = ¯ ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) s ( u n ( ξ ) − u ( ξ )) dξ (cid:12)(cid:12)(cid:12)(cid:12) ≤ ¯ ω (cid:90) | u n − u | ≤ ¯ ω n ∀ s ∈ [0 , . The role of the function θ (and similarly of θ n ) becomes clearer by remarking that from (23), we get q ss [ u ] = ¯ ω u (cos( θ ) , sin( θ )) a.e. in [0 , . By the assumption | u | , | u n | ≤
1, it follows that q ss [ u n ] · q ss [ u ] = ¯ ω u n u cos( θ − θ n ) ≥ ¯ ω u n u (1 −
12 ( θ − θ n ) ) ≥ ¯ ω u n u (cid:18) − ¯ ω n (cid:19) ≥ ¯ ω u n u − ¯ ω n . Now, recalling that q [ u ](0) = q [ u n ](0) = 0 and q s [ u ](0) = q s [ u n ](0) = − e , we have | q [ u n ](1) − q ∗ | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) (cid:90) s ( q ss [ u n ]( ξ ) − q ss [ u ]( ξ )) dξds (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) (cid:90) | q ss [ u n ] − q ss [ u ] | = (cid:90) (cid:90) | q ss [ u n ] | + | q ss [ u ] | − q ss [ u n ] · q ss [ u ] ≤ ¯ ω (cid:90) (cid:90) (cid:18) ( u n − u ) + ¯ ω n (cid:19) ≤ ¯ ω + ¯ ω n → n → + ∞ . We deduce from the arbitrariness of q ∗ that R (cid:48) ⊆ cl ( R ) and consequently int ( R (cid:48) ) ⊆ R . The proof iscomplete. (cid:3) We conclude this section by employing Proposition 3 to compute the boundary of the reachable set R .To this end, we choose a uniform grid on [0 ,
1] and we discretize the integral representation of the tentacletip q [ u ](1) = (cid:90) (cid:18) sin (cid:18)(cid:90) s ¯ ω u ( ξ ) dξ (cid:19) , − cos (cid:18)(cid:90) s ¯ ω u ( ξ ) dξ (cid:19)(cid:19) ds by means of a rectangular quadrature rule. Then we let u exhaustively attain all the extremal configurationsof type (24) projected on the grid. Figure 3 shows the results, depending on the constant parameter ¯ ω .¯ ω = π ¯ ω = π ¯ ω = π + 1 ¯ ω = 2 π ¯ ω = π Figure 3. reachable set depending on ¯ ω Starting from the extremal value ¯ ω = 0, for which (0 , −
1) is clearly the unique reachable point, the set R grows up to form a cardiod, enclosing a region of unreachable points. As ¯ ω increases, the hole disappears,while for ¯ ω → ∞ , R converges to the full circle of radius 1, which is exactly the reachable set of aninextensible string without curvature constraints.Numerical evidence shows that Proposition 3 holds also in the case of non constant curvature ¯ ω , as for thesimulation presented in the next section. Nevertheless, a complete proof of such result, based on Pontryagin’smaximum principle, is still under investigation.3.2. An augmented Lagrangian method.
Here we propose a numerical method for solving the stationaryoptimal control problem (22), and we present some simulations. First of all, we relax the inequality constrainton the curvature in (22), by introducing a so called slack variable , namely we define z := ¯ ω − | q ss | , so that | q ss | ≤ ¯ ω if and only if z ≥
0. Note that we still have an inequality constraint in z , but it is much simplerto treat, as we show below. We introduce the following augmented Lagrangian L ( q, σ, z, λ ) = J ( q ) + 12 (cid:90) σ ( | q s | − ds + 12 (cid:90) λ ( | q ss | − ¯ ω + z ) ds + 14 ρ λ (cid:90) ( | q ss | − ¯ ω + z ) ds , where(25) J ( q ) = 12 (cid:90) ω | q ss | ds + 12 τ | q (1) − q ∗ | ,σ is again an exact Lagrange multiplier for the inextensibility constraint, while λ and ρ λ > | q ss | − ¯ ω + z = 0. We apply theclassical method of multipliers [11, 23] to compute a solution of our minimization problem, namely we iterateon i ≥ ( q ( i +1) , σ ( i +1) , z ( i +1) ) = arg min q, σ, z ≥ L ( q, σ, z, λ ( i ) ) λ ( i +1) = λ ( i ) + 1 ρ λ ( | q ( i +1) ss | − ¯ ω + z ( i +1) )where the optimization sub-problem for fixed λ ( i ) is solved as follows, by imposing first order optimalityconditions. ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 17
Taking the variation of L with respect to σ in direction χ , we immediately recover the inextensibilityconstraint. Indeed, (cid:104) δ σ L , χ (cid:105) = (cid:90) ( | q s | − χds = 0 ∀ χ = ⇒ | q s | = 1 a.e. in (0 , . On the other hand, taking the variation of L with respect to z in direction v and imposing the constraint z ≥
0, we get the variational inequality (cid:104) δ z L , v (cid:105) = 12 (cid:90) (cid:18) λ ( i ) + 1 ρ λ ( | q ss | − ¯ ω + z ) (cid:19) ( v − z ) ds ≥ ∀ v ≥ , from which we readily obtain the optimal value z = max (cid:110) − λ ( i ) ρ λ − | q ss | + ¯ ω , (cid:111) a.e. in (0 , . Now, taking the variation of L with respect to q in direction w and substituting the above expression forthe optimal z , we get (cid:104) δ q L , w (cid:105) = (cid:90) (cid:16) Λ( q ss , λ ( i ) ) q ss · w ss + σq s · w s (cid:17) ds + 1 τ ( q (1) − q ∗ ) · w (1) = 0 ∀ w with Λ( q ss , λ ( i ) ) = 1¯ ω + max (cid:26) λ ( i ) + 1 ρ λ ( | q ss | − ¯ ω ) , (cid:27) . Integrating by parts and imposing boundary conditions, we end up with the following strong formulation ofthe optimality conditions:(26) (cid:0) Λ( q ss , λ ( i ) ) q ss (cid:1) ss − ( σq s ) s = 0 in (0 , | q s | = 1 in (0 , q (0) = 0 , q s (0) = − e q ss (1) = 0 , q sss (1) = 0 σ (1) q s (1) + τ ( q (1) − q ∗ ) = 0We remark that, apart the nonlinear term in Λ related to the curvature constraint, the first equation isexactly the stationary Euler’s elastica equation, completed with cantilevered boundary conditions. On theother hand, the last equation quantifies the balance between the tension of the tentacle and the boundaryforce resulting from the potential energy τ | q (1) − q ∗ | associated to the target point. More explicitly,dot-multiplying this equation by q s and using the inextensibility constraint, we obtain σ (1) = − τ ( q (1) − q ∗ ) · q s (1) , i.e., the tension at the free end is exactly opposite to the tangent projection of the boundary force.We discretize the optimality system (26) by means of a finite difference scheme as in Section 2.4, obtaining(27) D c (cid:0) Λ( D c q k , λ ( i ) ) D c q k (cid:1) − D + ( σ k D − q k ) = 0 k = 1 , ..., N − | D − q k | = 1 k = 1 , ..., N − q = 0 , q − = q + e ∆ sq N +1 − q N + q N − = 0 q N +1 − q N + 3 q N − − q N − = 0 σ N D − q N + τ ( q N − q ∗ ) = 0This is again a nonlinear system in the pair of unknowns ( q k , σ k ) k =1 ,...,N , and it can be easily solved bymeans of a quasi-Newton method. More precisely, at each iteration, the Newton step for (27) is computedby freezing the nonlinearity in the term Λ( D c q k , λ ( i ) ) at the previous iteration, namely by replacing the fullJacobian of D c (cid:0) Λ( D c q k , λ ( i ) ) D c (cid:1) with Λ( D c q k , λ ( i ) ) D c D c only. This is a reasonable approximation as weapproach a solution, indeed, by definition, Λ does not depend on q whenever the curvature constraint issatisfied. Finally, to evaluate convergence, we discretize the functional J in (25) by means of a simple rectangularquadrature rule:(28) J (cid:93) ( q ) = 12 ∆ s N (cid:88) k =0 ω k | D c q k | + 12 τ | q N − q ∗ | . Summing up, we have the following augmented Lagrangian algorithm:
Algorithm 1 Assign q (0) , λ (0) , a target point q ∗ , penalties τ, ρ λ > tol >
0. Compute J (cid:93) ( q (0) ) andset i = 0 repeat Compute the solution ( q ( i +1) , σ ( i +1) ) of (27) Update the multipliers λ ( i +1) k = max (cid:26) λ ( i ) k + 1 ρ λ ( | D c q ( i +1) k | − ¯ ω ) , (cid:27) for k = 1 , ..., N − . Compute J (cid:93) ( q ( i +1) ) and set i = i + 1 until (cid:12)(cid:12) J (cid:93) ( q ( i ) ) − J (cid:93) ( q ( i − ) (cid:12)(cid:12) < tol Note that this algorithm consists in two nested loops, an external loop for the update of the multipliers, andan internal loop for the computation of the solution. In practice, we realized that a single loop is enough toobtain and speed up the convergence, by performing both the update and the Newton step simultaneously.We choose the parameters for the simulations as in Section 2.4. Moreover, we set q (0) = − e s , λ (0) = 0, τ = 10 − , ρ λ = 10 and tol = 10 − . In Figure 4, we show the corresponding reachable set, and some optimalsolutions computed by Algorithm 1 for different target points. For each configuration we also report thegraph of the signed curvature κ = D − q × D c q (thick line), compared with the bounds [ − ¯ ω, ¯ ω ] (thin lines).In particular, we observe that the reachable set is a cardiod also in the case of non constant curvature ¯ ω .Moreover, we recognize a Dubins path in the top-right solution, formed by a CL -like curve as in (24), witha circular part of variable maximal radius 1 / ¯ ω plus a straight segment. -10-505100 0.2 0.4 0.6 0.8 1-10-505100 0.2 0.4 0.6 0.8 1-10-505100 0.2 0.4 0.6 0.8 1 -10-505100 0.2 0.4 0.6 0.8 1 Figure 4. solutions of the stationary optimal control problem for different target points
ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 19 The dynamic optimal control problem
In this section we tackle the dynamic optimal control problem discussed in the Introduction. Moreprecisely, given q ∗ ∈ R and T, τ >
0, we want to minimize the functional(29) J ( q, u ) = 12 (cid:90) T (cid:90) u ds dt + 12 τ (cid:90) T | q (1 , t ) − q ∗ | dt + 12 (cid:90) ρ ( s ) | q t ( s, T ) | ds , among all the controls u : [0 , × [0 , T ] → [ − , J accountrespectively for the activation of the tentacle muscles and the tip-target distance during all all evolution,whereas the last term corresponds to the kinetic energy at the final time T . The problem is then to reachand keep the tip close to the target point with minimum effort, but also to stop the whole tentacle as soonas possible. After deriving first order optimality conditions, we propose and implement an adjoint-basedgradient descent algorithm for its numerical solution.4.1. First order optimality conditions.
We introduce the adjoint state (¯ q, ¯ σ ), with ¯ q : [0 , × [0 , T ] → R ,¯ σ : [0 , × [0 , T ] → R , and we form the following Lagrangian¯ L ( q, σ, u, ¯ q, ¯ σ ) := J ( q, u ) + (cid:90) T (cid:90) ¯ q · (cid:8) ρq tt − ( σq s − Hq ⊥ ss ) s + ( Gq ss + Hq ⊥ s ) ss (cid:9) ds dt + 12 (cid:90) T (cid:90) ¯ σ ( | q s | − ds dt . By construction, the stationarity of ¯ L with respect to the adjoint variables immediately gives back theequations of motion and the inextensibility constraint. We then focus on the stationarity of ¯ L with respectto q , σ and u , which involves long and technical computations. Here we report the key steps, in particularthe derivation of the boundary conditions for the adjoint system. To this end, we recall, for the reader’sconvenience, the definitions of G and H given in (11), G [ q, ν, ε, ω ]( s, t ) := ε ( s ) + ν ( s ) (cid:0) | q ss ( s, t ) | − ω ( s ) (cid:1) + ,H [ q, µ, u, ω ]( s, t ) := µ ( s ) ( ω ( s ) u ( s, t ) − q s ( s, t ) × q ss ( s, t )) . and we denote the corresponding linearizations by¯ G [ q, ¯ q, ν, ω ]( s, t ) = g [ q, ν, ω ]( s, t ) q ss ( s, t ) · ¯ q ss ( s, t ) , ¯ H [ q, ¯ q, µ ])( s, t ) = − µ ( s ) (¯ q s ( s, t ) × q ss ( s, t ) + q s ( s, t ) × ¯ q ss ( s, t )) , where g [ q, ν, ω ]( s, t ) = 2 ν ( s ) ( | q ss ( s, t ) | − ω ( s )) and ( · ) stands for the Heaviside function, i.e. ( x ) = 1 for x ≥ ( x ) = 0 otherwise. As usual, we assume enough regularity to perform the computations, and weemploy the given boundary conditions to define the admissible tests ( w, χ ) associated to ( q, σ ). Moreover,we recall the assumptions (10) on the functions ε, µ and some useful properties of G, H . We summarize allthese relations in the following list: for ( s, t ) ∈ [0 , × [0 , T ](30) q ( s,
0) = q ( s ) ⇒ w ( s,
0) = 0 q t ( s,
0) = v ( s ) ⇒ w t ( s,
0) = 0 q (0 , t ) = 0 ⇒ w (0 , t ) = 0 q s (0 , t ) = − e ⇒ w s (0 , t ) = 0 q ss (1 , t ) = 0 ⇒ w ss (1 , t ) = 0 q sss (1 , t ) = 0 ⇒ w sss (1 , t ) = 0 σ (1 , t ) = 0 ⇒ χ (1 , t ) = 0 ε ( s ) ≥ ε > ⇒ G (1 , t ) = ε (1) > µ (1 , t ) = µ s (1 , t ) = 0 ⇒ H (1 , t ) = H s (1 , t ) = 0Let us start by computing the variation of ¯ L with respect to σ in direction χ : (cid:104) δ σ ¯ L , χ (cid:105) = − (cid:90) T (cid:90) ¯ q · ( χq s ) s ds dt = − (cid:90) T [¯ q · q s χ ] dt + (cid:90) T (cid:90) ¯ q s · q s χ ds dt . Using (30) and imposing stationarity, by the arbitrariness of χ , we get for almost every ( s, t ) ∈ (0 , × (0 , T )(31) ¯ q (0 , t ) · q s (0 , t ) = 0 and ¯ q s ( s, t ) · q s ( s, t ) = 0 , that we prolong by continuity to all ( s, t ) ∈ [0 , × [0 , T ].We now compute the variation of ¯ L with respect to q in direction w : (cid:104) δ q ¯ L , w (cid:105) = (cid:90) T ( q (1 , t ) − q ∗ ) · w (1 , t ) dt + (cid:90) ρq t ( s, T ) · w t ( s, T ) ds + (cid:90) T (cid:90) ρ ¯ q · w tt ds dt − (cid:90) T (cid:90) ¯ q · ( σw s + ¯ H [ q, w ] q ⊥ ss + Hw ⊥ ss ) s ds dt + (cid:90) T (cid:90) ¯ q · ( ¯ G [ q, w ] q ss + Gw ss + ¯ H [ q, w ] q ⊥ s + Hw ⊥ s ) ss ds dt + (cid:90) T (cid:90) ¯ σq s · w s ds dt , Using again (30), after a very long integration by parts, we obtain(32) (cid:104) δ q ¯ L , w (cid:105) = (cid:90) T (cid:90) (cid:8) ρ ¯ q tt − (cid:0) σ ¯ q s − H ¯ q ⊥ ss + ¯ σq s − ¯ Hq ⊥ ss (cid:1) s + (cid:0) G ¯ q ss + H ¯ q ⊥ s + ¯ Gq ss + ¯ Hq ⊥ s (cid:1) ss (cid:9) · w ds dt − (cid:90) ρ ¯ q t · w ( s, T ) ds + (cid:90) ρ ( q t + ¯ q ) · w t ( s, T ) ds + (cid:90) T ε ¯ q ss · w s (1 , t ) dt + (cid:90) T (cid:26) ¯ σq s − ( G ¯ q ss ) s + 1 τ ( q − q ∗ ) (cid:27) · w (1 , t ) dt − (cid:90) T (cid:8) g (¯ q · q ss ) q ss + G ¯ q + µ (¯ q · q ⊥ s ) q ⊥ s (cid:9) · w sss (0 , t ) dt + (cid:90) T (cid:8) H ¯ q ⊥ − G s ¯ q − ¯ q · (2 µq ⊥ ss + µ s q ⊥ s ) q ⊥ s − ¯ q · ( gq ss ) s q ss − g (¯ q · q ss ) q sss + g (¯ q s · q ss ) q ss + G ¯ q s + µ (¯ q s · q ⊥ s ) q ⊥ s (cid:9) · w ss (0 , t ) dt . Imposing stationarity, by the arbitrariness of w , it follows that all the integrands in the expression aboveshould vanish for almost every ( s, t ) ∈ (0 , × (0 , T ). In order, we get the adjoint equation ρ ¯ q tt = (cid:0) σ ¯ q s − H ¯ q ⊥ ss + ¯ σq s − ¯ Hq ⊥ ss (cid:1) s − (cid:0) G ¯ q ss + H ¯ q ⊥ s + ¯ Gq ss + ¯ Hq ⊥ s (cid:1) ss , and the final conditions ¯ q ( s, T ) = − q t ( s, T ) , ¯ q t ( s, T ) = 0 . Moreover, since ε (1) > q ss (1 , t ) = 0 , and then we also get(33) ¯ σ (1 , t ) q s (1 , t ) − ε (1)¯ q sss (1 , t ) + 1 τ ( q (1 , t ) − q ∗ ) = 0 . ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 21
Now consider the relation ¯ q s · q s = 0 in (31). Differentiating twice in s , we obtain q s · ¯ q sss +2 q ss · ¯ q ss + q sss · ¯ q s = 0,and using the boundary conditions q ss (1 , t ) = q sss (1 , t ) = 0, it follows that q s (1 , t ) · ¯ q sss (1 , t ) = 0. Dot-multiplying (33) by q s , we then conclude¯ σ (1 , t ) = − τ ( q (1 , t ) − q ∗ ) · q s (1 , t ) . Substituting in (33) and projecting ( q − q ∗ ) on the base { q s , q ⊥ s } , we also get¯ q sss (1 , t ) = 1 τ ε (1) (cid:0) ( q − q ∗ ) · q ⊥ s (cid:1) q ⊥ s (1 , t ) . We now proceed with the stationarity of the last two integrands in (32), related to the boundary conditionsat s = 0. The first one gives(34) g (¯ q · q ss ) q ss + G ¯ q + µ (¯ q · q ⊥ s ) q ⊥ s = 0 . Using the orthogonality relation ¯ q (0 , t ) · q s (0 , t ) = 0 in (31), for every t ∈ [0 , T ] we can write ¯ q = αq ⊥ s forsome α ∈ R . On the other hand, due to the relation q s · q ss = 0, we can also write q ss = βq ⊥ s for some β ∈ R .Substituting in (34), we get α ( gβ + G + µ ) q ⊥ s = 0 . Since
G > g, µ ≥
0, we deduce α = 0 and, consequently¯ q (0 , t ) = 0 . Substituting in (32), several terms cancel out and the last integrand reduces to g (¯ q s · q ss ) q ss + G ¯ q s + µ (¯ q s · q ⊥ s ) q ⊥ s = 0 . Note that this expression is exactly (34) with ¯ q s in place of ¯ q . Then, using the relation ¯ q s (0 , t ) · q s (0 , t ) = 0in (31), we can repeat the same argument above to conclude¯ q s (0 , t ) = 0 . We finally compute the variation of ¯ L with respect to the control, assuming for a moment that u isunconstrained, and considering smooth admissible tests v : [0 , × [0 , T ] → R . We get (cid:104) δ u ¯ L , v (cid:105) = (cid:90) T (cid:90) uv ds dt + (cid:90) T (cid:90) ¯ q · { ( µωvq ⊥ ss ) s + ( µωvq ⊥ s ) ss } ds dt . Integrating by parts and using µ (1 , t ) = µ s (1 , t ) = 0, ¯ q (0 , t ) = ¯ q s (0 , t ) = 0, we easily get(35) (cid:104) δ u ¯ L , v (cid:105) = (cid:90) T (cid:90) (cid:0) u + ω ¯ H [ q, ¯ q ] (cid:1) v ds dt , i.e. δ u ¯ L = u + ω ¯ H [ q, ¯ q ] in L ((0 , × (0 , T )) . Summarizing, we end up with the following strong formulation of the adjoint system(36) ρ ¯ q tt = (cid:0) σ ¯ q s − H ¯ q ⊥ ss (cid:1) s − (cid:0) G ¯ q ss + H ¯ q ⊥ s (cid:1) ss in (0 , × (0 , T )+ (cid:0) ¯ σq s − ¯ Hq ⊥ ss (cid:1) s − (cid:0) ¯ Gq ss + ¯ Hq ⊥ s (cid:1) ss ¯ q s · q s = 0 in (0 , × (0 , T )¯ q (0 , t ) = 0 t ∈ (0 , T )¯ q s (0 , t ) = 0 t ∈ (0 , T )¯ q ss (1 , t ) = 0 t ∈ (0 , T )¯ q sss (1 , t ) = τε (cid:0) ( q − q ∗ ) · q ⊥ s (cid:1) q ⊥ s (1 , t ) t ∈ (0 , T )¯ σ (1 , t ) = − τ ( q − q ∗ ) · q s (1 , t ) t ∈ (0 , T )¯ q ( s, T ) = − q t ( s, T ) s ∈ (0 , q t ( s, T ) = 0 s ∈ (0 , u should satisfy, for every v : [0 , × [0 , T ] → [ − , (cid:90) T (cid:90) (cid:0) u + ω ¯ H [ q, ¯ q ] (cid:1) ( v − u ) ds dt ≥ . It is easy to see that, in case of dissipation, the frictional forces − βq t and − γq sssst in the equations of motion(15), simply produce in (36) the analogous (but with opposite signs) terms β ¯ q t and γ ¯ q sssst respectively.4.2. An adjoint-based gradient descent method.
From a numerical point of view, solving the optimal-ity system obtained in the previous section is a complicate task, since we have to find an optimal 5-tuple( q ∗ , σ ∗ , ¯ q ∗ , ¯ σ ∗ , u ∗ ) satisfying (14), (36) and (37) in the whole space-time interval [0 , × [0 , T ]. After dis-cretization, this may result in a huge nonlinear system. Hence we adopt another (and simpler) approach,namely we use an iterative method that, at each iteration, splits the solution in three separate steps. Moreprecisely, given a guess for the control u , we first solve the system of motion (14) with u fixed, up to thefinal time. Then we solve the adjoint system (36) backward in time, again with u fixed and using the justcomputed ( q, σ ) as a datum. Finally, we plug q and ¯ q in the expression (35) of the gradient δ u ¯ L , and updatethe control u along the corresponding descent direction. The procedure is iterated until convergence on thecontrol. For the discretization of both forward and backward systems, we employ the Velocity Verlet schemepresented in Section 2.4, whereas the box constraints on the control are enforced, at each iteration, by theprojection on the interval [ − , [ − , ( · ). Summarizing, we build the following algorithm.Let us discuss how to build a suitable initial guess u (0) for the algorithm above. Once q ∗ is fixed, we Algorithm 2 Assign u (0) , a target point q ∗ , a penalty τ >
0, a tolerance tol >
0, a step size α > i = 0 repeat Compute the solution ( q ( i ) , σ ( i ) ) of (14) Compute the solution (¯ q ( i ) , ¯ σ ( i ) ) of (36) Compute δ u ¯ L ( q ( i ) , ¯ q ( i ) ) using (35) Update the control u ( i +1) = Π [ − , (cid:0) u ( i ) − αδ u ¯ L ( q ( i ) , ¯ q ( i ) ) (cid:1) and set i = i + 1 until (cid:13)(cid:13) u ( i ) − u ( i − (cid:13)(cid:13) ∞ < tol employ Algorithm 1 for the stationary control problem. By construction, this provides an equilibrium ˜ q ( s )minimizing the distance of the tentacle tip from q ∗ , and we can easily synthesize the corresponding optimalcontrol ˜ u ( s ) using the relation ˜ u ( s ) = ω ( s ) ˜ q s ( s ) × ˜ q ss ( s ). Then, we define u (0) extending ˜ u to a functionon [0 , × [0 , T ] which is constant in time, i.e. we set u (0) ( s, t ) := ˜ u ( s ) for all t ∈ [0 , T ], and we referto it as to the static optimal control . It is clear that, if we plug u (0) in the controlled dynamical system(14), the corresponding evolution will convergence to ˜ q only in presence of friction and for T → ∞ (see thesimulations in Section 2.4). Here the aim is to show that Algorithm 2 can do much better, producing a dynamic optimal control that changes both in space and time. To this end, we choose the parameters for thesimulations as in Section 2.4, with the exception of the friction coefficients, that we set as β ( s ) = 2 − s and γ ( s ) = 10 − (2 − s ). In this way, we accentuate the elastic behavior of the tentacle, and we can focus on thereduction of oscillations performed by the optimization. Moreover, we set T = 4, τ = 10 − , tol = 10 − and q ∗ = (0 . , − . N = 10 to have a reasonable computational time,and we use a fixed gradient descent step size α for a simple code implementation. In this setting, we foundthat α = 10 − produces an almost monotone decrease of the functional J in (29) up to convergence of thecontrol. Finally, for a qualitative analysis of the results, it is useful to compute the following discretizationsof the three terms appearing in J , J u ( n ) = ∆ s N (cid:88) k =0 ( u nk ) , J q ∗ ( n ) = 12 τ | q nN − q ∗ | , J v ( n ) = ∆ s N (cid:88) k =0 ρ k | v nk | , respectively the control energy, the target energy and the kinetic energy at time t n = n ∆ t .Figure 5 shows the resulting evolution at different times. For comparison, we represent in each frame theprofile associated to the dynamic optimal control (black line), the profile associated to the static optimalcontrol (gray line), and the equilibrium configuration (dashed line) corresponding to the target point q ∗ (small circle). Moreover, Figure 6 and Figure 7 show the behavior in time of the energies J v , J q ∗ and J u ,respectively for the static and the dynamic controls. ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 23
While the static control keeps the tentacle oscillating around the target point, waiting for the slow frictiondissipation, the optimized dynamics is really rich and deserves some comments. We observe that, at the verybeginning, the dynamic control strongly activates to compensate the high acceleration given by the targetpotential, then relaxes to reduce the control energy (0 ≤ t ≤ . . ≤ t ≤ . µ . This suggests thatthe tip stabilization can be faster as N → ∞ , and we can also expect a smoother behavior of the kineticenergy. Finally, we observe that, as t → T = 4, the dynamic control energy converges to the static one.This is not a mere consequence of the choice for the initial guess u (0) , since we know that the static controlis optimal if the kinetic energy vanishes. t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 1 . t = 1 . t = 1 . t = 1 . t = 2 . t = 2 . t = 2 . t = 3 . t = 3 . t = 4 . Figure 5.
Static optimal control vs dynamic optimal control
References [1] P. K. Agarwal and H. Wang. Approximation algorithms for curvature-constrained shortest paths.
SIAM Journal onComputing , 30(6):1739–1772, 2001.[2] A. M. Bruckstein, R. J. Holt, and A. N. Netravali. Discrete elastica. In
International Conference on Discrete Geometryfor Computer Imagery , pages 59–72. Springer, 1996.[3] A. Burchard and L. E. Thomas. On the cauchy problem for a dynamical euler’s elastica.
Communications on PartialDifferential Equations , 28(1-2):271–300, 2003.[4] E. J. Cockayne and G. W. C. Hall. Plane motion of a particle subject to curvature constraints.
SIAM Journal on Control ,13(1):197–220, 1975. (a) (b) (c)
Figure 6.
Time evolution of kinetic energy J v (a), target energy J q ∗ (b) and controlenergy J u (c) for the static optimal control (a) (b) (c) Figure 7.
Time evolution of kinetic energy J v (a), target energy J q ∗ (b) and controlenergy J u (c) for the dynamic optimal control [5] A. Della Corte, F. Dell’Isola, R. Esposito, and M. Pulvirenti. Equilibria of a clamped euler beam (elastica) with distributedload: Large deformations. Mathematical Models and Methods in Applied Sciences , 27(08):1391–1421, 2017.[6] Y. Engel, P. Szabo, and D. Volkinshtein. Learning to control an octopus arm with gaussian process temporal differencemethods. In
Advances in neural information processing systems , pages 347–354, 2006.[7] S. D. Fisher and J. W. Jerome. Stable and unstable elastica equilibrium and the problem of minimum curvature. In
Minimum Norm Extremals in Function Spaces , pages 90–106. Springer, 1975.[8] P. Fritzkowski and H. Kaminski. Dynamics of a rope as a rigid multibody system.
Journal of mechanics of materials andstructures , 3(6):1059–1075, 2008.[9] P. Fritzkowski and H. Kaminski. Non-linear dynamics of a hanging rope.
Latin American Journal of Solids and Structures ,10(1):81–90, 2013.[10] A. Fedotov, V. Patsko, and V. Turova. Reachable sets for simple models of car motion. In
Recent Advances in MobileRobotics . InTech, 2011.[11] M.R. Hestenes. Multiplier and gradient methods.
Journal of Optimization Theory and Applications , 4:303–320, 1969.[12] B. A. Jones and I. D. Walker. Kinematics for multisection continuum robots.
IEEE Transactions on Robotics , 22(1):43–55,2006.[13] R. Kang, A. Kazakidi, E. Guglielmino, D. T. Branson, D. P. Tsakiris, J. A. Ekaterinaris, and D. G. Caldwell. Dynamicmodel of a hyper-redundant, octopus-like manipulator for underwater applications. In
Intelligent Robots and Systems(IROS), 2011 IEEE/RSJ International Conference on , pages 4054–4059. IEEE, 2011.[14] A. Kazakidi, D. P. Tsakiris, D. Angelidis, F. Sotiropoulos, and J. A. Ekaterinaris. Cfd study of aquatic thrust generationby an octopus-like arm under intense prescribed deformations. 115:54–65, 2015.[15] R. Levien. The elastica: a mathematical history.
University of California, Berkeley, Technical Report No. UCB/EECS-2008-103 , 2008.[16] Anna Chiara Lai, Paola Loreti, et al. Robot’s finger and expansions in non-integer bases.
Networks and HeterogeneousMedia , 7(1):71–111, 2012.[17] Paola Loreti and Anna Chiara Lai. Robot’s hand and expansions in non-integer bases.
Discrete Mathematics & TheoreticalComputer Science , 16(1), June 2014.[18] Anna Chiara Lai, Paola Loreti, and Pierluigi Vellucci. A continuous fibonacci model for robotic octopus arm. In , pages 99–103. IEEE, 2016.
ODELING AND OPTIMAL CONTROL OF AN OCTOPUS TENTACLE 25 [19] Anna Chiara Lai, Paola Loreti, and Pierluigi Vellucci. A fibonacci control system with application to hyper-redundantmanipulators.
Mathematics of Control, Signals, and Systems , 28(2):15, 2016.[20] C. Laschi, B. Mazzolai, V. Mattoli, M. Cianchetti, and P. Dario. Design of a biomimetic robotic octopus arm.
Bioinspiration& biomimetics , 4(1):015006, 2009.[21] A. A. Markov. Some examples of the solution of a special kind of problem on greatest and least quantities.
Soobshch.Karkovsk. Mat. Obshch , 1:250–276, 1887.[22] S. Montgomery-Smith and W. Huang. A numerical method to model dynamic behavior of thin inextensible elastic rods inthree dimensions.
Journal of Computational and Nonlinear Dynamics , 9(1):011015, 2014.[23] M.J.D. Powell.
A method for nonlinear constraints in minimization problems . Academic Press, New York, NY, 1969.[24] S. C. Preston. The motion of whips and chains.
Journal of Differential Equations , 251(3):504–550, 2011.[25] D. Rus and M. T. Tolley. Design, fabrication and control of soft robots.
Nature , 521(7553):467–475, 2015.[26] German Sumbre, Yoram Gutfreund, Graziano Fiorito, Tamar Flash, and Binyamin Hochner. Control of octopus armextension by a peripheral motor program.
Science , 293(5536):1845–1848, 2001.[27] F. Tr¨oltzsch.
Optimal control of partial differential equations: theory, methods, and applications , volume 112. AmericanMathematical Soc., 2010.[28] D. Yong. Strings, chains, and ropes.
SIAM review , 48(4):771–781, 2006.[29] Y. Yekutieli, R. Sagiv-Zohar, R. Aharonov, Y. Engel, B. Hochner, and T. Flash. Dynamic model of the octopus arm. i.biomechanics of the octopus reaching movement.
Journal of neurophysiology , 94(2):1443–1458, 2005.[30] Y. Yekutieli, R. Sagiv-Zohar, B. Hochner, and T. Flash. Dynamic model of the octopus arm. ii. control of reachingmovements.
Journal of neurophysiology , 94(2):1459–1468, 2005.(S. Cacace)
Dipartimento di Matematica e Fisica, Universit`a degli Studi Roma Tre, Largo S. Murialdo, 1, 00154Roma, Italy, [email protected] (A. C. Lai)
Dipartimento di Scienze di Base e Applicate per l’Ingegneria, Sapienza Universit`a di Roma, Via A.Scarpa, 16, 00161 Roma, Italy, [email protected] (P. Loreti)(P. Loreti)