Application of autonomous pathfinding system to kinematics and dynamics problems by implementing network constraints
AApplication of autonomous pathfinding system tokinematics and dynamics problems byimplementing network constraints
Kei-Ichi Ueda ∗ Abstract
A neural network system in an animal brain contains many modulesand generates adaptive behavior by integrating the outputs from the mod-ules. The mathematical modeling of such large systems to elucidate themechanism of rapidly finding solutions is vital to develop control meth-ods for robotics and distributed computation algorithms. In this article,we present a network model to solve kinematics and dynamics problemsfor robot arm manipulation. This model represents the solution as anattractor in the phase space and also finds a new solution automaticallywhen perturbations such as variations in the end position of the arm orobstacles occur. In the proposed model, the physical constraints, tar-get position, and the existence of obstacles are represented by networkconnections. Therefore, the theoretical framework of the model remainsalmost the same when the number of constraints increases. In addition,as the model is regarded as a distributed system, it can be applied towardthe development of parallel computation algorithms.
Over the past few decades, trajectory control for robot arm manipulation hasgarnered considerable attention in robot engineering (Featherstone, R. & Orin,D. (2000); Rodriguez, G., Jain, A., & Kreutz-Delgado, K. (1992)). A funda-mental kinematics problem in such a control system is to determine a feasibleposition for the arm joint such that the end point approaches the target loca-tion. This issue is formulated in terms of boundary value problem (BVP) intwo-dimensional space, which can be expressed as follows:BVP: Find x l ∈ Q l ( l = 2 , , . . . , L −
1) that satisfies the following conditions: | x l − x l +1 | = d l , ( l = 1 , . . . , L − x = ( x s , y s ) ∈ Q , x L = ( x g , y g ) ∈ Q L , where Q l ⊂ R is the feasible region for the joint l and L ≥ x l corresponds to the position of l th joint in the two-dimensionalspace. As BVP has multiple solutions in general, we need to solve the singu-lar equation of x l . In addition, if the obstacles are assumed to be placed inthe domain or Q l (cid:40) R , additional constraints should be implemented in theequation. Numerical algorithms have been proposed for kinematics problemsbased on iterative methods (Aristidou, A. & Lasenby, J. (2011); Unzueta, L.,Peinado, M., Boulic, R., & Suescun, ´A. (2008)) and neural network models(Tejomurtula, S. & Kak, S. (1999); K¨oKer, R. (2013); Toshani, H. & Farrokhi, ∗ Faculty of Science, Academic Assembly, University of Toyama, Toyama 930-8555, Japan a r X i v : . [ n li n . AO ] D ec M. (2014)). For dynamics problems, we need to consider additional constraintsto obtain a smooth arm motion. Several approaches have been proposed forsolving dynamics problems such as optimization methods (Wada, Y., Koike, Y.,Vatikiotis-Bateson, E., & Kawato, M. (1994); Poggio, T. & Girosi, F. (1990)),self-organizing maps (Kuperstein, M. (1988); Walter, J. A. & Schulten, K. I.(1993)), neural network models (Wada, Y. & Kawato, M. (1993); Narendra, K.S. & Parthasarathy, K. (1990); Glasius, R., Komoda, A., & Gielen, S. C. A. M.(1995)), and reservoir computation (Polydoros, A. S., & Nalpantidis, L. (2016)).As the number of components in the system increases, the formulation ofthe algorithms becomes complicated in general. In this study, we propose anew framework for modeling kinematics and dynamics problems. Here, thesolution of BVP is represented by a path in the network connecting nodes,which correspond to the boundary values. Further, the constraints in the armlength and the presence of obstacles are described by the addition and removalof the network links. Thus, the network construction procedure remains almostthe same as the number of system’s components and that of constraints increase.In fact, we need to attach or detach the network links according to the physicalconstraints and the position of obstacles.Autonomy is an important concept that should be considered while con-structing a robot system with adaptive behavior (Volpe, R., Nesnas, I., Estlin,T., Mutz, D., Petras, R., & Das, H. (2001)). When humans encounter un-familiar environment, they autonomously develop strategies and execute newactions. There is an increasing demand for the development of such an au-tonomous system, whose control algorithm is based on just the variables of thesystem. In addition, as the number of components in the system increases, thedistributed processing is required to decrease the computation time. Therefore,the development of effective autonomous and distributed systems has receivedconsiderable attention in industries. If the system is formulated in terms ofdifferential equations, the flexibility of the system against environmental vari-ation can be regarded as the switching of attractor in the phase space. Thus,the elucidation of the mathematical mechanism for the robustness of attractorswitching is vital to improve the performance of the system. Ueda et al. (Ueda,K. I., Yadome, M., & and Nishiura, Y. (2015)) proposed a network model toshow flexible attractor switching. This network model has been applied topathfinding problems and shows the following properties: (1) The model canspontaneously find one of the possible paths connecting two target points. (2)It begins to find another path when perturbations such as removal of pathsoccur. Using the above properties and implementing network constraints, weconstruct a solver that can autonomously find a solution of BVP and finds an-other possible solution when the existing solution becomes impractical due toperturbations.To apply the pathfinding model to BVP, we formulate a discretized versionof BVP, which is called DBVP. We define a two-dimensional lattice in (cid:101)
Ω and aset of the lattice points Γ as follows: (cid:101) Ω l := { x ∈ R | x = ( ξ i,j , η i,j ) , ( i, j ) ∈ Γ } , l ∈ { , . . . , L } ,ξ i,j := x min + ( i − · ( x max − x min ) / ( J x − ,η i,j := y min + ( j − · ( y max − y min ) / ( J y − , Γ := { ( i, j ) | i = 1 , . . . , J x , j = 1 , . . . , J y } . Because we assume that the solutions are attained at the lattice points, weformulate the DBVP as follows:DBVP: Find X l ∈ (cid:101) Q l ⊂ (cid:101) Ω l ( l = 2 , , . . . , L −
1) that satisfies the following condi-tions: | X l − X l +1 | ∈ [ d l − ∆ d l , d l + ∆ d l ] , ( l = 1 , . . . , L −
1) (1) X = ( x s , y s ) ∈ (cid:101) Q , X L = ( x g , y g ) ∈ (cid:101) Q L . In general, due to the discretization, we need to consider the margin ∆ d l as theconstrain for d l . The value of ∆ d l is determined by the geometrical constraintand can be reduced if J x and J y increase.Firstly, we apply the pathfinding model to DBVP. The boundary values forthe base position X and the end position X L are given as the start and targetpoint in the network. The boundary condition and the physical constraint forthe robot arm are described by the network. Secondly, we extend the DBVPmodel to the dynamics problem. As the network contains excitatory and in-hibitory connections between the nodes and integration system does not exist,the model represents a distributed system. Therefore, our study is potentiallyuseful for the development of the parallel computation algorithms to solve kine-matics and dynamics problems. We apply the model proposed in Ueda, K. I., Yadome, M., & and Nishiura, Y.(2015), which can find one of the possible paths connecting the start and targetpoints in hierarchical network consisting of nodes and directional excitatory andinhibitory links. The node dynamics is described by differential equations. Thesolution path is described by a stationary state of the model.
The network-construction procedure of the pathfinding system based on a hi-erarchical network is shown in Fig. 1. We assume that the start and targetpoints of the network are at the top and bottom layer, respectively. Excitatoryand inhibitory links are attached according to the following rules:(P1) The nodes corresponding to the point k in the network (Fig. 1(a)) areplaced at the P and N layers (Fig. 1(b)) . These nodes are called node k + and node k − , respectively.(P2) The excitatory links directed from node m + to k + and from m − to k − are attached (Fig. 1(b)) if a connection exists between point m and k (Fig. 1(a)). The existence of excitatory interaction directed from nodes m ± to k ± is represented as a m ± ,k ± , where a m ± ,k ± = 1 and a m ± ,k ± = 0indicate the presence and absence of such interactions, respectively.(P3) There are inhibitory links from nodes m + , l + , and l − to node k + and fromnodes m + , k + , and k − to node l + if there are excitatory links from node m + to both the nodes k + and l + . Similarly, there are inhibitory linksfrom nodes m − , l − and l + to node k − and from nodes m − , k − and k + tonode l − if there are excitatory links from node m − to the nodes k − and l − .According to the above procedure, the activated state, which is defined as ONstate, propagates from top to bottom in the P layer and from bottom to thetop in the N layer.We add excitatory links between P and N layers at the start and target nodesto form a loop. Thus, the solution path connecting the start and target nodesis represented by the nodes with ON state forming a loop network architecture.We use two different descriptions for the target point, which are discussed inSec. 5.(P4) The network has an excitatory link from node k − s to node k + s and fromnode k + g to node k − g , where k + s ( k − s ) and k + g ( k − g ) indicate the nodes at thestart and target positions in the P (N) layer, respectively. The boundarycondition for the start position is expressed in terms of the link connectionˆ a k − ,k + = (cid:40) k ± = k ± s a k + ,k − = (cid:40) k ± = k ± g , (2)ˇ a k = (cid:40) k = k − g . (3) The pathfinding models with boundary conditions given by equations (2) and(3) are referred to as Model I and Model II, respectively. According to the rules(P1) - (P4), Model I is described by˙ u k ± = f ( u k ± , v k ± ) + µ H (cid:32) (cid:88) m ± ∈ Λ ± a m ± ,k ± H ( u m ± − θ ) (cid:33) + µ ˆ a k ∓ ,k ± H ( u m ± − θ ) − µ H (cid:32) (cid:88) m ± ∈ Λ ± (cid:88) l ∈ Λ a m ± ,l ± a m ± ,k ± H ( u m ± − θ ) H ( u l ± + u l ∓ − θ ) (cid:33) + Aσ k ± ( t ) , ˙ v k ± = g ( u k ± , v k ± ) . (4)Similarly, Model II is described by˙ u k + = f ( u k + , v k + ) + µ H (cid:32) (cid:88) m + ∈ Λ ± a m + ,k + H ( u m + − θ ) (cid:33) + µ ˆ a k − ,k + H ( u m ± − θ ) − µ H (cid:32) (cid:88) m + ∈ Λ + (cid:88) l ∈ Λ a m + ,l + a m + ,k + H ( u m + − θ ) H ( u l + + u l − − θ ) (cid:33) + Aσ k + ( t ) , ˙ v k + = g ( u k + , v k + ) , ˙ u k − = f ( u k − , v k − ) + µ H (cid:32) (cid:88) m − ∈ Λ − a m − ,k − H ( u m − − θ ) (cid:33) + µ ˇ a k − − µ H (cid:32) (cid:88) m − ∈ Λ − (cid:88) l ∈ Λ a m − ,l − a m − ,k − H ( u m − − θ ) H ( u l − + u l + − θ ) (cid:33) + Aσ k − ( t ) , ˙ v k − = g ( u k − , v k − ) , (5)where k ± = 1 ± , ± , . . . , K ± , Λ ± := { ± , ± , . . . , K ± } , Λ := { , , . . . , K } , and τ ∈ [0 , T max ] is dimensionless time. The dot above u and v indicates theirderivative with respect to τ , and µ i and θ i ( i = 1 ,
2) are positive constants. Thefunctions f and g are described by the sigmoid FitzHugh–Nagumo equation(Rotstein, H. G., Kopell, N., Zhabotinsky, A. M., &. Epstein, I. R. (2003)): f ( u, v ) = − au + bu − cu + d − v , g ( u, v ) = (cid:15) [ p tanh (( u − q ) /r ) + s − v ].We set ( a, b, c, d, p, q, r, s, (cid:15) ) = (1 . , . , . , . , . , . , . , . , . σ k ± represents Gaussian noise with zero mean and unit variance,and A is the noise amplitude. We assume that the interaction function H isthe Heaviside step function and has a threshold θ such that the connectivityof the link switches dynamically: H ( u − θ ) = 1 for u > θ and H ( u − θ ) = 0otherwise. The third term on the right-hand side of u k ± -equation ensures thatthe system finds a single path from the multiple feasible solutions with the sameroute at the P and N layer. Based on this formulation, the node k + receives aninhibitory input when u m + > θ and u l + + u l − > θ . For each node, we defineON, OFF1, and OFF2 states depending on u k . The definition of the node statesand the process for finding the solutions are described in A and B. Here, we callthe OFF1 and OFF2 states as OFF state. To recapitulate, the system describedby equations (4) and (5) has the following properties: • The system robustly finds one of the possible paths connecting the startand target nodes if the solution exists. The associated nodes in P and Nlayers acquire ON state. • All the nodes acquire OFF state if no solution exists. This implies thatthe system can terminate the search process if no possible solution exists. • The system automatically starts the search process when the existing pathis damaged, and also terminates the search process when it finds a newpossible path.The representative time sequences of u k + for Model I and Model II are shownin Fig. 2. It is evident that the models successfully find solutions and exhibitflexible attractor switching when the target position is changed.Figure 1: (a) An example of hierarchical network. The start and target nodesare placed at the top and bottom layer of the hierarchy. (b) A network modelfor the network in (a). Arrows indicate the excitatory links. The direction ofthe excitatory link in the N layer is opposite to that in the P layer. Additionalexcitatory links are added at the start and target nodes.Figure 2: Numerical results obtained by Model I (Left) and II (Right). (a) Anexample of network structure, where only the excitatory links in the P layer areshown. The target position is considered at node 17 for τ < T max and at 20for τ > T max . Circles represent the node state when the solution converges tostationary state for τ < T max and τ > T max . Black and white circles indicatethe ON and OFF states, respectively. (b) Numerical solution of equation (4) forthe network. Only the time sequences of the nodes in P layer are shown. Thehorizontal and vertical directions indicate u k + and τ , respectively. It is clearthat the system successfully finds a possible solution when the target positionis varied.Figure 3: Schematic of excitatory link network. (a) According to (P2 (cid:48) ), exci-tatory links are connected from l to l + 1 in the P layer and from l + 1 to l in N layer. To form a loop structure of excitatory links, the links representingthe start and end positions are installed from N to P layer and from P to Nlayer, respectively. (b) Excitatory links are attached if the corresponding nodessatisfy the constraint in equation (1). The solution of DBVP is represented by a path connecting the start and targetnodes, which correspond to the base and end points, respectively. The nodes areplaced on a two-dimensional square lattice and the solution of DBVP, i.e., X l isrepresented by the position of nodes in ON state. The physical constraint andthe existence of the obstacles are described by link connection and disconnection.The network consists of L pairs of lattices and the l th pair is used to representthe position of l th joint. The lattice size is J x × J y , where J x and J y representthe grid size of the x - and y -coordinates, respectively. Thus, the resolutionof the approximation method is improved as J x and J y increase. Nodes arelocated at every lattice point. Therefore, the total number of nodes is 2 J x J y L .The minimum and maximum values of the x -coordinate ( y -coordinate) are x min ( y min ) and x max ( y max ), respectively. For notational convenience, the serialnumber of each node in given in each layer. To distinguish whether the nodebelongs to P or N layer, the k th node in P and N layer is called k + and k − node, respectively. We define a set of the serial number of nodes located in the l th lattice in P and N layer as Λ + ( l ) and Λ − ( l ), respectively, i.e.,Λ ± ( l ) = (cid:8) [1 + ( l − · J x J y ] ± , [2 + ( l − · J x J y ] ± , . . . , [ J x J y + ( l − · J x J y ] ± (cid:9) . The x - and y -value of the node k ± ( i, j, l ) ∈ Λ ± ( l ) are defined by ξ k ± = ξ k ± ( i,j,l ) = ( ξ k ± ( i,j,l ) , η k ± ( i,j,l ) ) ,ξ k ± ( i,j,l ) := x min + ( i − · ( x max − x min ) / ( J x − ,η k ± ( i,j,l ) := y min + ( j − · ( y max − y min ) / ( J y − . Model I is applied to DBVP by modifying (P2) and (P4) as follows:(P2 (cid:48) ) The links in P layer are attached if the two corresponding nodes satisfythe physical constraint in equation (1) (Fig. 3). This implies that a m + ,k + = | ξ m + − ξ k + | ∈ [ d l − ∆ d l , d l + ∆ d l ] , m + ∈ Λ + ( l ) ,k + ∈ Λ + ( l + 1) , l ∈ [1 , . . . , L − a k − ,m − ( k − ∈ Λ − ( l + 1) , m − ∈ Λ − ( l ))are determined according to (P2).(P4 (cid:48) ) The node numbers for the start (target) node in P and N layer are denotedas k + s ( k + g ) and k − s ( k − g ), respectively. According to (P4), the excitatorylinks are attached from the target node in P layer to the target node in Nlayer and from the start node in the N layer to the target node in P layer,i.e., ˆ a k − ,k + = (cid:40) k ± = k ± s , ˆ a k + ,k − = (cid:40) k ± = k ± g , where k ± s ∈ Λ ± (1) and k ± g ∈ Λ ± ( L ) are the boundary conditions. Due tothe inhibitory interaction, only a single pair of nodes acquire ON state forevery l ∈ { , . . . , L } when the system finds a solution.Due to these assumptions, the solutions of the model necessarily satisfy thephysical constraints and boundary conditions. We consider the following cases as perturbations: (1) variation of boundary value X L during computation, and (2) the existence of obstacles. In this section, forFigure 4: (a) Time sequence of the node state in P and N layers. The circlesindicate ON state. (b) The positions of the joints are displayed in the lattice(black: τ = 490 < T max /
2, gray: τ = T max = 1000). Here, L = 4 and l = 1 . / ( L −
1) = 0 . l th joint is displayed in the lattice. (a) Case 1. (b)Case 2. (c) Case 3. The gray regions in (a) and (b) indicate the forbidden region (cid:101) Q c for all the joints and for joints 4 and 5, respectively. The feasible region forjoint 4 is indicated by circles.simplicity, we consider that x min = y min = 0, x max = y max = 1, and J x = J y =21. The Euler–Maruyama method is used for time integration, where the timegrid is set as ∆ τ = 0 .
01. For obtaining the approximate solution, ∆ d l should bedetermined so that the union of the circles covers the entire region of the solutionspace. We consider ∆ d l = √ x (= √ y ) so that the model represented byequation (4) can robustly find a solution, where ∆ x = ∆ y = 1 /J x = 0 .
1. Theparameters are set as θ = 1 . θ = 3 . µ = 1 . µ = 9 .
0, and A = 1 . × − .As an initial state, the node at the start position in the P layer is considered tobe in ON state and the other nodes are in OFF state for all experiments. To confirm that the system can flexibly find a new solution, the boundary value X L is changed during the computation. The boundary values X and X L aregiven by X = (0 , , X L = (cid:40) (1 ,
0) for τ < T max / , (1 ,
1) otherwise , where T max = 1000. Initially, ON state propagates from layer 1 to L in the Player and then from layer L to 1 in the N layer (Fig. 4 (a1)). The systemsuccessfully finds one of the possible solutions before τ = T max /
2, i.e., only onenode is in ON state at every layer (Fig. 4 (a2)). The transient dynamics isobserved just after the position of X L is changed (Fig. 4 (a3)). The systemsuccessfully finds a new solution before τ = T max (Fig. 4 (a4)). Here, we consider the case in which obstacles are placed in the system. Theexistence of the obstacles is represented by the removal of excitatory links di-rected to the nodes located at the positions of obstacles. We assume that thelinks emanating from the nodes in the forbidden region are removed in the Player, and the links directed to the nodes in the forbidden region are removedin the N layer. This implies that a k + ,m + = a m − ,k − = 0 , if ( ξ k + , η k + ) ∈ (cid:101) Q cl ( l = 1 , · · · , L ) , where (cid:101) Q cl := (cid:101) Ω \ (cid:101) Q l .Other connections are determined according to (P2 (cid:48) ) and (P3). We examinenumerical results for the following three cases. For all the cases, we considerthat L = 6, l = 0 . T max = 500, X = (0 , X L = (1 . , . Case 1
The forbidden region is given by (cid:101) Q cl = { ( x, y ) | x ∈ [0 . , . , y ∈ [0 , . } , l ∈ { , . . . , L } Case 2
The forbidden region for the joints 6 and 7 is given by (cid:101) Q cl = (cid:40) { ( x, y ) | x ∈ [0 , , y ∈ [0 , . } , if l = 6 , ∅ , otherwise . Case 3
This case corresponds to the combination of the constrains in the above cases.Some joints are restricted to a specific position or to a specific region. Forexample, we consider the following constraint: (cid:101) Q cl = (cid:40) (cid:101) Ω \{ ( x, y ) | x = 0 . , y = 0 . } , if l = 4 ∅ , otherwise . This implies that (cid:101) Q = { ( x, y ) | x = 0 . , y = 0 . } . It may be noted that(0 . , . ∈ (cid:101) Q l . Figure 5 shows the numerical solution for these three cases. Itis clear that the model successfully finds one of the possible solutions satisfyingthe constraints. The variation in the position of obstacle is represented by the attachment andremoval of links. As established in earlier studies (Ueda, K. I., Yadome, M.,& and Nishiura, Y. (2015)), the proposed model can autonomously find a newsolution when the network structure varies during the computation. Owing tothis property, the system autonomously begins to find new solution when anobstacle destroys the existing solution. We assume that the obstacle motion isdescribed by (cid:101) Q cl ( τ ) = { ( x, y ) ∈ R | x ∈ [0 . , . ∩ (cid:101) Ω l , y ∈ ([0 , . ϕ τ ] ∪ [0 .
45 + ϕ τ , L y ]) ∩ (cid:101) Ω l } , X = (0 , ∈ (cid:101) Ω , X L = (1 . , . . ϕ τ ) ∈ (cid:101) Ω L ,ϕ τ = 0 . × (cid:106) τ (cid:107) , (6)where (cid:98) x (cid:99) = max { n ∈ Z | n ≤ x } . Figure 6 proves that the model autonomouslystarts to find another solution when the existing solution enters the forbiddenregion, and every joint corresponding to the new solution enters the feasibleregion.Figure 6: Transition of the solution of DBVP when the obstacles move accordingto equation (6). The gray regions indicate the forbidden region. Transientdynamics is initiated when τ = 10000 × ( m + 1) ( m = 0 , , . . . ). The modelfinds a new solution before t reaches 10000 × m + 5000 for every m . Here, L = 4and l = 1 . / ( L −
1) = 0 .
5. Solution for (a) m = 0, (b) m = 2, (c) m = 4, (d) m = 6, and (e) m = 8.1 Here, we measure the rate of increase in the number of steps when L is in-creased. The number of steps is defined as N step when the solution convergesto a stationary state or when a single pair of the nodes acquires ON state forevery pair in the layers. This implies that if the solution is found in time T (cid:48) ,then N step = T (cid:48) / ∆ τ . The positions ( x , y ) and ( x L , y L ) are fixed at (0 , i, j ) = (1 , , .
5) [( i, j ) = (21 , L , which is confirmed inFig. 7.Figure 7: Average number of steps required to find a solution across 20 trialsas a function of L . We extend the model for the motion problem to solve the orbit problem. It isnot guaranteed that the system represented by equation (4) exhibits a smoothmotion when the target position is given because only physical and boundaryconditions are employed as constraints. Therefore, to apply our algorithm forgenerating smooth trajectories of robot arm, additional constraint with respectto the continuity of feasible regions is required, which is called time constraint .We consider the following DBVP with time constraint (DBVPT).DBVPT: Find X t,l ∈ (cid:101) Q l ( t = 1 , . . . , T ; l = 1 , . . . , L ) that satisfies the followingconditions: | X t,l − X t,l +1 | ∈ [ d l − ∆ d l , d l + ∆ d l ] , ( t = 1 , . . . , T ; l = 1 , . . . , L − , (7) | X t,l − X t +1 ,l | ∈ [0 , ∆ b ] , ( t = 1 , . . . , T − l = 1 , . . . , L ) , (8) X t, = ( x s ( t ) , y s ( t )) ∈ (cid:101) Q , X t,L = ( x g ( t ) , y g ( t )) ∈ (cid:101) Q L , (9)where equations (7), (8), and (9) represent the physical constraint, time con-straint, and boundary condition, respectively. ∆ b represents the maximumspeed of arm motion between t and t + 1. Small ∆ b generates smooth armmotion. Schematics of matrix and network structure are shown in Fig. 8. The2P and N layers at the l th column and t th row of the matrix, which represent theposition of l th joint at time t , are denoted by ( t, l ) + and ( t, l ) − , respectively.The indexes of the variables are changed as follows: k ± ( i, j, l, t ) = ( i + ( j − J x + ( l − J x J y + ( t − J x J y L ) ± , ξ k ± ( i,j,l,t ) := ( ξ k ± ( i,j,l,t ) , η k ± ( i,j,l,t ) ) ,ξ k ± ( i,j,l,t ) := x min + ( i − x max − x min ) / ( J x − ,η k ± ( i,j,l,t ) := y min + ( j − y max − y min ) / ( J y − , Λ ± ( t, l ) = { (1 + ( l − J x J y + ( t − J x J y L ) ± , · · · , ( J x J y + ( l − J x J y + ( t − J x J y L ) ± } , Λ( t, l ) = { (1 + ( l − J x J y + ( t − J x J y L ) , · · · , ( J x J y + ( l − J x J y + ( t − J x J y L ) } , (cid:101) Ω := { x ∈ R | x = ξ k ( i,j, · , · ) , ( i, j ) ∈ Γ } , The node number at matrix ( i, j ) in ( t, l ) + and ( t, l ) − is denoted by k + ( i, j, l, t )and k − ( i, j, l, t ), respectively. The feasible region restricted by the time con-straint is defined as B + ( m + ( i, j, l, t )) = (cid:8) k + ∈ Λ + ( t, l ) | | ξ k + − ξ m + | ∈ [0 , ∆ b ] , m + ∈ Λ + ( t + 1 , l ) (cid:9) ,B − ( m − ( i, j, l, t )) = (cid:8) k − ∈ Λ − ( t, l ) | | ξ k − − ξ m − | ∈ [0 , ∆ b ] , m − ∈ Λ − ( t − , l ) (cid:9) . We define respectively a set of node positions for time t that satisfy the boundaryconditions and physical constraint as (cid:101) S ± B,P ( t ) ⊂ Λ ± ( t, × · · · × Λ ± ( t, L ) =: (cid:101) Λ ± ( t ) . A set of plausible solutions satisfying the time constraint is defined as follows: (cid:101) B + ( k + t ) := B + ( k + t, ) × · · · × B + ( k + t,L ) , k + t = ( k + t, , . . . , k + t,L ) ∈ Λ + ( t, × · · · × Λ + ( t, L ) . The network connection procedures (P2 (cid:48) ) and (P4 (cid:48) ) are modified as (P2 (cid:48)(cid:48) )and (P4 (cid:48)(cid:48) ), and a procedure (P5 (cid:48)(cid:48) ) for the time constraint is added.(P2 (cid:48)(cid:48) ) The links in P layer are attached if the corresponding two nodes satisfythe condition in equation (1), i.e., a m + ,k + = | ξ m + − ξ k + | ∈ [ d l − ∆ d l , d l + ∆ d l ] , m + ∈ Λ + ( t, l ) ,k + ∈ Λ + ( t, l + 1)0 otherwise ,l ∈ , . . . , L − , t ∈ , . . . , T. The connections in N layer, a k − ,m − ( k − ∈ Λ − ( t, l ) , m − ∈ Λ − ( t + 1 , l )) aredetermined according to (P2).(P4 (cid:48)(cid:48) ) We denote the node number for the start (target) node in P and N layerat time t as k + s ( t ) ( k + g ( t )) and k − s ( t ) ( k − g ( t )), respectively. According to(P4 (cid:48) ), we attach excitatory links from the target node in P layer to thetarget node in N layer and from the start node in N layer to the targetnode in P layer, i.e.,ˆ a k − s ( t ) ,k + s ( t ) = 1 , ˆ a k + g ( t ) ,k − g ( t ) = 1 , t = 1 , . . . , T (P5 (cid:48)(cid:48) ) Excitatory links are attached from node k + t,l ∈ Λ + ( t, l ) to k + t − ,l ∈ Λ + ( t − , l ) if k + t − ,l ∈ B + ( { k + t,l } ). Similarly, excitatory links are attached fromnode k − t ∈ Λ − ( t, l ) to k − t +1 ,l ∈ Λ − ( t + 1 , l ) if k t +1 ,l ∈ B − ( { k t,l } ). Thepresence and absence of the connection is represented by b m + ,k + = (cid:40) k + ∈ B + ( { m + } ) , ,b m − ,k − = (cid:40) k − ∈ B − ( { m − } ) , ,i = 1 , . . . , J x ; j = 1 , . . . , J y . T, l ) ± to (1 , l ) ± ( l ∈ [1 , . . . , L ]). Thisimplies that the system can quickly find a solution under the condition (C).(C) For any t ∈ [2 , . . . , T ] and element k + t ∈ (cid:101) S +B,P ( t ), there exists k + t − ∈ (cid:101) S +B,P ( t − ∩ (cid:101) B + ( k + t ) (Fig. 9(b)(c)).We use Model I for t = 1 and Model II for t = 2 , . . . , T . We assume thatthe node k receives an excitatory signal when a m (cid:48) ,k b m (cid:48)(cid:48) ,k = 1, and both thenodes m (cid:48) and m (cid:48)(cid:48) are in the ON state, i.e., the node can be in the ON state ifit satisfies the physical and the time constraints. The model can be expressedas follows:˙ u k + = f ( u k + , v k + ) + µ G (cid:32) (cid:88) m + ∈ Λ + a m + ,k + H ( u m + − θ ) (cid:33) × G (cid:32) (cid:88) m + ∈ Λ + b m + ,k + H ( u m + − θ ) (cid:33) + µ ˆ a k − ,k + H ( u k − − θ ) − µ G (cid:32) (cid:88) m + ∈ Λ + (cid:88) l ∈ Λ ( a m + ,l + a m + ,k + + b m + ,l + b m + ,k + ) H ( u m + − θ ) H ( u l + + u l − − θ ) (cid:33) + Aσ k + ( t ) , ˙ v k + = g ( u k + , v k + ) , ˙ u k − = f ( u k − , v k − ) + µ G (cid:32) (cid:88) m − ∈ Λ − a m − ,k − H ( u m − − θ ) (cid:33) × G (cid:32) (cid:88) m − ∈ Λ − b m − ,k − H ( u m − − θ ) (cid:33) + µ (cid:101) H ( u ˆ k − θ ; l − ) − µ G (cid:32) (cid:88) m − ∈ Λ − (cid:88) l ∈ Λ ( a m − ,l − a m − ,k − + b m − ,l − a m − ,k − ) H ( u m − − θ ) H ( u l − + u l + − θ ) (cid:33) + Aσ k − ( t ) , ˙ v k − = g ( u k − , v k − ) , (10)where k ± = k ± ( i, j, l, t ) ∈ Λ ± := ∪ l,t Λ ± ( t, l ) and (cid:101) H ( x ˆ k ; l ) = (cid:40) H ( x ˆ k ) (if l = L − )0 (otherwise) , ˆ k = ˆ k ( i, j, l, t ) = (cid:40) k + g (1) (if l = L − and t = 1) (Model I) k − g ( t −
1) (if l = L − and t ≥
2) (Model II) . The initial data is taken such that the node at the target position ξ k − g (1) isin ON state. The ON state propagates according to the following sequence:(i) The ON state propagates from (1 , L ) − to (1 , − layer and from (1 , L ) − to ( T, L ) − layer (Fig. 8(c)).(ii) The ON state propagates from (1 , l ) − layer to ( T, l ) − layer and from( t, L ) − to ( t, − layer (Fig. 8(c)).(iii) The ON state propagates from ( T, l ) + layer to (1 , l ) + layer and from ( t, + to ( t, L ) + layer (Fig. 8(d)).(iv) After the ON state reaches ( T, L ) + , the nodes k ± T ∈ (cid:101) S ± B,P ( T ) are selectedduring the process (iii) (Fig. 8(d)).(v) Due to (C), there exists k ± T − ∈ (cid:101) S ± B,P ( T − ∩ (cid:101) B + ( k + T ).(vi) The process (v) successively occurs for t = T − , . . . , L = 4 and T = 4,where the start and target points are given by( x k ± s ( t ) , y k ± s ( t ) ) = (0 ,
0) ( t = 1 , . . . , T )and ( x k ± g (1) , y k ± g (1) ) = (1 , , ( x k ± g (2) , y k ± g (2) ) = (0 . , . , ( x k ± g (3) , y k ± g (3) ) = (0 . , . , ( x k ± g (4) , y k ± g (4) ) = (0 . , . . We consider that θ = 1 . θ = 3 . µ = 1 . µ = 9 . A = 1 . × − ,∆ b = 0 .
25, and J x = J y = 11. Model I can be easily applied to the kinematics problem under the case thatthe joint number increases or decreases during computation. Such situationsoccur when a tool is being handled or an additional joint is attached. Forexample, when the joint number increases from L to L + 1, we add L + 1 th layer as well as excitatory and inhibitory links between L th and L + 1 th layersaccording to (P2 (cid:48) ), (P3), and (P4 (cid:48) ). Further, we remove the existing connectionat the target point in L th layer and attach new excitatory connection at thenew target point in L + 1 th layer. In numerical simulation, we add a segmentwith length l = 2 . / ( L −
1) = 0 . τ = 500. L is changed from 4 to 5, and a ij is set according to (P2 (cid:48) ) and (P3) during computation. The new target positionin L + 1 th layer is ( x, y ) = (1 . , . J x and J y . However, the number of steps increases exponentially as J x and J y increase.Another potential way is to successively decrease the area of the search regionas the model finds a solution, but J x and J y should be constant. Here, weset x min , x max , y min , y max to x l min (0) , x l max (0) , y l min (0) , y l max (0) and replace themby x l min (1) , x l max (1) , y l min (1) , y l max (1), respectively, such that | x l max (1) − x l min (1) | and | y l max (1) − y l min (1) | become smaller than | x l max (0) − x l min (0) | and | y l max (0) − y l min (0) | as the model finds a solution for a given x l min (0) , x l max (0) , y l min (0) , y l max (0).This algorithm is expressed as follows:Step 0 Set n ←
0, 0 < r < x l min ( n ) , x l max ( n ) , y l min ( n ), and y l max ( n ).Step 1 Find a solution by using Model I. Let the solution be X c = ( x c , y c ) and n ← n + 1.Step 2 x l min ( n ) , x l max ( n ) , y l min ( n ) , y l max ( n ) are given by∆ d l ← r ∆ d l ,x l min ( n ) = x c − ( J x − d l / , x l max ( n ) = x c + ( J x − d l / ,y l min ( n ) = y c − ( J y − d l / , y l max ( n ) = y c + ( J y − d l / . We stop the computation if | X l − X l +1 | become smaller than the expectedprecision for all l = 1 , . . . , L −
1. Otherwise, we return to step 1.The robustness of the calculation can be enhanced if r approaches 1, but theiteration time between steps 1 and 2 increases. In future, we hope to derive theoptimal value of r for which the model robustly finds a solution at every stepwith minimum iteration steps.5Figure 8: Schematic of the network model. (a) The arrows indicate the excita-tory links in the P layers. Gray circles indicate the time constraint B − ( k , ).(b) Three-dimensional network structure for the excitatory links of (a). (c) Thegray circles indicate that only the nodes in the N layer are in ON state. (Step1) Due to the initial data, the node at the target point in (1 , L ) − layer acquiresON state. (Step 2) The node at (1 , L − − and (2 , L ) − layer acquires ON state.(Step 3) As both (1 , L − − and (2 , L ) − layer are in the ON state, the node at(2 , L − − also acquires ON state. In addition, the node at (1 , L − − acquiresON state. (Step 4) As both (1 , L − − and (2 , L − − layer are in ON state,the node at (2 , L − − acquires ON state. (d) The black circles indicate thatnodes in both P and N layers are in ON state. (Step 1) The node at the startpoint of ( t, + layer (1 ≤ t ≤ T −
1) can be in ON state when both the nodesat the start in ( t, − and ( t + 1 , + acquire ON state, but ( T, + acquires ONstate when ( T, − is in ON state. Thus, the node at the start point ( T, + acquires ON state faster than the other nodes in P layer. (Step 2) Similar tothe case of N layer, a possible pair in ( T, + and ( T − , + acquires ON state.(Step 3) Subsequently, a possible pair in ( T, + and ( T − , + acquires ONstate. (Step 4) As both ( T, + and ( T − , + layer are in ON state, the nodeat ( T − , + acquires ON state.6Figure 9: (a) Concept of the model represented by equation (10). (b) Blacklines indicate the solutions that satisfy the physical and time constraints. Timeconstraint satisfies the condition (C). In fact, point P is in B + ( m + T,l − ), andpoint Q is also in B + ( T − l, m + T,l ). (Right) Time constraint does not satisfythe condition (C). In fact, point P is in B + ( m + T,l − ), but point Q is not in B + ( m + T,l ).7Figure 10: Solution of equation (10) displayed on x – y plane.In our numerical experiments, we mainly considered the constraints for theposition of the joints. We can formulate the model such that every segment doesnot enter the forbidden regions. Figure 12 shows an example of the solution ofDBVP when the connections are determined as follows: a k + ,m + = (cid:40) | ξ k + − ξ m + | ∈ [ d l − ∆ d l , d l + ∆ d l ] , and the line segment connecting ξ k + and ξ m + does not enter (cid:101) Q l , where k + ∈ Λ + ( l ), m + ∈ Λ + ( l + 1), l = 1 , . . . , L − X l for (a) τ = 400 and (b) τ = 1000. Anew 4 th link is attached at τ = 500. Acknowledgements
This work was supported by JSPS KAKENHI grant numbers 18H04940 and17K05361.
References
Aristidou, A. & Lasenby, J. (2011). FABRIK: A fast, iterative solver for theinverse kinematics problem. Graphical Models, 73(5), 243-260.8Figure 12: A possible solution when a ij is set such that each segment does nottouch the forbidden region. The gray region indicates the forbidden region. Weconsider that X = (0 , X L = (1 , . L = 6, and l = 2 . / ( L −
1) = 0 . k when the state of node m isgiven. (b) Examples of the possible stationary states at the branching point ofthe excitatory links. Stationary state of the node k + when the states of nodes m + , l + , and l − are given. (c) Solution finding process of the proposed model.The arrows indicate the excitatory links. ON waves initially propagate along theexcitatory links. Inhibitory interaction occurs after an ON wave in the N layerreaches the node 1 − (right-most panel). The left and right halves of the circleindicate the nodes in the P layer and the N layer, respectively. (d) Recoveryprocess when the target position is changed from node 6 ± to 7 ± . When thenode 5 − acquires OFF1 state, post-inhibitory rebound occurs at the node 2 ± ,and nodes 3 + , 4 + , 4 − , 3 − and 2 − successively acquire ON state.9Featherstone, R. & Orin, D. (2000, April). Robot dynamics: equationsand algorithms. In Proceedings 2000 ICRA. Millennium Conference. IEEEInternational Conference on Robotics and Automation. Symposia Proceed-ings (Cat. No. 00CH37065) (Vol. 1, pp. 826-834). IEEE.Glasius, R., Komoda, A., & Gielen, S. C. A. M. (1995). Neural net-work dynamics for path planning and obstacle avoidance. Neural Networks,8(1), 125-133.K¨oKer, R. (2013). A genetic algorithm approach to a neural-network-based inverse kinematics solution of robotic manipulators based on errorminimization. Information Sciences, 222, 528-543.Kuperstein, M. (1988). Neural model of adaptive hand-eye coordina-tion for single postures. Science, 239(4845), 1308-1311.Narendra, K. S. & Parthasarathy, K. (1990). Identification and controlof dynamical systems using neural networks. IEEE Transactions on NeuralNetworks, 1(1), 4-27.Poggio, T. & Girosi, F. (1990). Regularization algorithms for learningthat are equivalent to multilayer networks. Science, 247(4945), 978-982.Polydoros, A. S., & Nalpantidis, L. (2016, October). A reservoir com-puting approach for learning forward dynamics of industrial manipulators. In2016 IEEE/RSJ International Conference on Intelligent Robots and Systems(IROS) (pp. 612-618). IEEE.Rodriguez, G., Jain, A., & Kreutz-Delgado, K. (1992). Spatial opera-tor algebra for multibody system dynamics. Journal of the AstronauticalSciences, 40(1), 27-50.Rotstein, H. G., Kopell, N., Zhabotinsky, A. M., &. Epstein, I. R.(2003). A canard mechanism in systems of globally coupled oscillators. SIAMJournal on Applied Mathematics, 63(6), 1998-2019.Tejomurtula, S. & Kak, S. (1999). Inverse kinematics in robotics usingneural networks. Information Sciences, 116(2-4), 147-164.Toshani, H. & Farrokhi, M. (2014). Real-time inverse kinematics of re-dundant manipulators using neural networks and quadratic programming:A Lyapunov-based approach. Robotics and Autonomous Systems, 62(6),766-781.Ueda, K. I., Yadome, M., & and Nishiura, Y. (2015). Multistate net-work model for the pathfinding problem with a self-recovery property. NeuralNetworks, 62, 32-38.Unzueta, L., Peinado, M., Boulic, R., & Suescun, ´A. (2008). Full-bodyperformance animation with sequential inverse kinematics. Graphical Mod-els, 70(5), 87-104.Wada, Y., Koike, Y., Vatikiotis-Bateson, E., & Kawato, M. (1994). Acomputational model for cursive handwriting based on the minimizationprinciple. In Advances in Neural Information Processing Systems (pp.727-734).0Wada, Y. & Kawato, M. (1993). A neural network model for arm tra-jectory formation using forward and inverse dynamics models. NeuralNetworks, 6(7), 919-932.Walter, J. A. & Schulten, K. I. (1993). Implementation of self-organizingneural networks for visuo-motor control of an industrial robot. IEEE Trans-actions on Neural Networks, 4(1), 86-96.Volpe, R., Nesnas, I., Estlin, T., Mutz, D., Petras, R., & Das, H.(2001, March). The CLARAty architecture for robotic autonomy. In 2001IEEE Aerospace Conference Proceedings (Cat. No. 01TH8542) (Vol. 1, pp.1/121-1/132). IEEE. A Single node dynamics and definition of nodestates
The definition of node states and the procedure for determining the parametersin equation (4) are based on an earlier study (Ueda, K. I., Yadome, M., &and Nishiura, Y. (2015)), The dynamics of an isolated node is described by thefollowing equation: ˙ u = f ( u, v ) + µ I e − µ I i , ˙ v = g ( u, v ) . (11)where f ( u, v ) = − au + bu − cu + d − v , g ( u, v ) = (cid:15) [ p tanh (( u − q ) /r ) + s − v ].We set ( a, b, c, d, p, q, r, s, (cid:15) ) = (1 . , . , . , . , . , . , . , . , . u -nullcline ( f = 0) and v -nullcline ( g = 0) be-comes sufficiently small at the peak point of u -nullcline ( u = u p in Fig. 14(b)).The terms I e and I i correspond to the excitatory and inhibitory inputs, re-spectively, and they are either 0 or 1. In equation (4), due to the networkconstruction, each node receives one of the following signals: ( I e , I i ) = (0 , , , I e , I i ) ≡ (0 , , ,
0) acquire OFF1, OFF2, and ON states, respectively. The pa-rameters θ , θ , µ , and µ are set such that the system exhibits PIR and H ( u m ± − θ ) H ( u l ± + u l ∓ − θ ) = 1 when the nodes m ± , l ± , and l ∓ are in ONstate. B Pathfinding system
B.1 State transition
The fundamental role of the excitatory link is to propagate ON state and that ofthe inhibitory link is to select a solution at the branching point of the excitatorylinks. Suppose that the node k receives an excitatory link from node m . Thenode k acquires ON state if the node m is in ON state and acquires OFF1 statefor the remaining cases (Fig. 13(a)). According to (P3), inhibitory links areinstalled at the branching point of the excitatory links. Typical cases are shownin Fig. 13(b). It is evident that the node k + acquires OFF2 state when thenodes m + , l + , and l − are in ON state. For the remaining cases (ii)–(vi), thenode k + does not receive inhibitory signals. The search and selection processis schematically shown in Fig. 13(a). We assume that node 1 + is in ON stateand the others are in OFF1 state. ON state propagates along 1 + → + → + and 1 + → + → + , and then along 5 − → − → − in the N layer. Inhibition1Figure 14: (a) Time sequence of the solution for the single-node model rep-resented by equation (12), where I e and I i are varied according to equation(13). (b) The red line indicates the trajectory of the solution shown in (a). Thesolid black lines denoted by (A), (B), and (C) correspond to u -nullcline when( I e , I i ) = (1 , , , v -nullcline, which is independent of I e and I i .occurs when the nodes in P and N layers at point 4 ± acquire ON state and thatat point 2 ± acquire OFF2 state. B.2 Postinhibitory rebound
PIR is crucial for the self-recovery property of the model. According to ournetwork construction procedure, PIR occurs at the branching point of the ex-citatory links. Thus, the OFF2 state is observed when the corresponding nodereceives both excitatory and inhibitory signals. We show this fundamental be-havior by using a simple model with external forces corresponding to excitatoryand inhibitory signals. ˙ u = f ( u, v ) + µ I e , ˙ v = g ( u, v ) + µ I i , (12)where I e and I i are defined as follows: I e = I i = (cid:40) t ∈ [2000 , , , (13)For τ ∈ [0 , τ ∈ [2000 , B.3 Recovery process
Figure 13(b) shows the recovery process of the model (4) when the target posi-tion is changed. After the target position is changed, nodes 7 − , 6 − , 5 − , 1 − , and1 + in the N layer successively acquire OFF1 state. When the node 5 − acquiresOFF1 state, the nodes 2 + and 2 − acquire ON state due to PIR, and ON statepropagates along 2 + → + → + → − → − and 2 − → − → ++