Computers and Electrical Engineering

A mobile robot acting in a human environment should follow social conventions, keeping safety distances and navigating at moderate speeds, in order to respect people in its surroundings and avoid obstacles in real-time. The problem is more complex in differential-drive wheeled robots, with trajectories constrained by nonholonomic and kinematics restrictions. It is an NP-hard problem widely studied in the literature, combining disciplines such as Psychology, Mathematics, Computer Science and Engineering. In this work, we propose a novel solution based on Membrane Computing, Social Force Model and Dynamic Window Approach Algorithm. The resulting model is able to compute, in logarithmic time, the best motion command for the robot, given its current state, considering the surrounding people and obstacles. The model is compatible with other membrane computing models for robotics and suitable for an implementation on parallel hardware. Finally, a visual simulator was implemented in ROS and C++ for validation and testing.


Introduction
The problem of motion planning in robotics can be described in the configuration space of a robot as follows: Given an initial robot configuration state and a goal configuration state, find a sequence of motion commands to move the robot towards the goal state while avoiding static and dynamic obstacles. In general, the motion planning problem is crucial in robotics and other areas such as verification, computational biology, and computer animation [1]. In particular, a critical case is given when a wheeled or legged robot is walking in a human environment populated with people, for example a pedestrian street, an office, a hospital or a marketplace. In this situation, the robot must not only avoid obstacles while trying to reach its goal, but also generate socially acceptable trajectories. An extensive survey on human-aware navigation can be found in [2]. Moreover, it implies several philosophical questions concerning the challenges of the application of these technologies to real-life scenarios. In the case a robot has to choose irretrievably between injuring one out of two persons, what should it do? And who would be the one condemned for such an act? [3]. On the other hand, several applications for these technologies can be found in the literature [4]. The problem is complex and involves multidisciplinary fields including Psychology, Mathematics, Computer Science and Engineering. Furthermore, even if the robot's optimal velocity vector could be computed for each step of time, there are several physical constraints to be considered, such as acceleration limits and nonholonomic restrictions.
In this paper, we present the first solution based on membrane computing to the motion planning problem in human environments. Membrane computing [5] is a branch of natural computing inspired in the structure and functions of the living cells. The computational devices in membrane computing are called membrane systems or P systems, and they have been successfully applied to a large variety of scientific areas: the modeling of microscopic and macroscopic biological systems [6,7]; and the study of neural models [8,9] incorporating fuzzy reasoning [10], among others. In particular, a variant of P systems called Enzymatic Numerical P systems (ENPS) has been used to simulate robot controllers [11,12] and to simulate robot navigation algorithms [13,14].
One of the most important features of any robot navigation system is the need to provide real-time solutions to complex situations. Therefore, most of the classical algorithms can only provide approximate solutions given a fixed time of computation. Several authors have studied how to accelerate such algorithms by using parallel hardware [15]. On the other hand, an intermediate computing paradigm can be used. In this sense, the membrane computing paradigm provides an inherently parallel computing framework with a large variety of simulators that can emulate computations over parallel hardware [16,17]. Thus, ENPS can be used to model solutions to robot navigation problems, and well-studied hardware simulators can apply the solutions to real case studies. The main advantage of this approximation is to use the computational power of membrane computing, which is inherently parallel, as well as widely-studied and robust simulation algorithms for parallel hardware. On the other hand, the main disadvantage is the complication of adapting such algorithms, requiring expert knowledge of both robotics and membrane computing.
The problem of motion planning can be divided in two sub-problems [18]: the global planning problem, in which only static obstacles are considered, and the local planning problem, in which the robot must follow a path while avoiding static and dynamic obstacles. If the dynamic obstacles are given by moving people, then it is a social local planning problem, since the robot should execute socially acceptable trajectories. In this paper, we provide a novel ENPS model for social local planning with a computational complexity of ( ( )) based on the Social Force Model (SFM) [19,20] and the Dynamic Window Approach Algorithm (DWA) [21]. The model is compatible with global planning ENPS models such as [14] and can be simulated in parallel hardware. A software simulator has been implemented for validation and testing purposes by using the Robot Operating System (ROS) 1 and C++. The simulator introduces in a virtual environment a number of human agents walking, following the SFM and a dual-wheeled nonholonomic robot, navigating according to the proposed ENPS model. Average velocities and distances to obstacles and other agents are measured for all the simulated agents to provide some results and discussion.
The rest of this work is structured as follows: Section 2 introduces the necessary preliminaries to understand the paper. Section 3 presents the proposed social local planning algorithm. Section 4 explains the corresponding ENPS model. Section 5 presents the simulation tool. Section 6 is devoted to analysis and discussion. Finally, Section 7 enumerates some conclusions and future work.

Preliminaries
This Section provides the reader with the basic concepts and notations used throughout this paper, from the description of the problem addressed in this study to the methods involved in the proposed solution.

Global and local planning
The problem of motion planning in robotics can be partially solved by a software architecture in three levels [18]: global planner, local planner and PID controller. Firstly, the global planner is the module in charge of solving the motion planning problem without considering dynamic obstacles (usually surrounding people). The output of the global planner is a trajectory or plan representing the path from the starting point to the goal area. This process is usually carried out by using a pre-computed map of static obstacles, i.e., walls, furniture, and other fixed obstacles that have been previously tracked. The computation is done before the robot starts to move by applying algorithms such as variants of RRT [22,23], and RRT* [24], among others. After global planning has been performed, the local planner module generates motion commands to follow the trajectory in a safe manner, by considering the information given by the robot sensors in real time. The Dynamic Window Approach for Obstacle Avoidance Algorithm [21] and the Pure Pursuit Algorithm [25] are well-known algorithms for local planning. Finally, a PID controller module [26] manages the power of the motors to fit each motion command and maintain a constant velocity until the next command is processed.

Differential wheeled robots
Depending on the motion constraints imposed, two types of mobile robots can be distinguished: holonomic and nonholonomic. On the one hand, holonomic robots are those that can move in any direction in the configuration space. On the other hand, nonholonomic robots are systems where the velocities (magnitude and/or direction) and other derivatives of the position are constrained. In this paper, we consider the type of differential wheeled robots in which each wheel is driven independently; in particular, we focus on two-wheeled robots. A differential wheeled robot can be characterized in a 2D space by using a tuple ( , , ), known as its 2D pose, where , are the Cartesian coordinates of its center and is its yaw angle, i.e., the orientation of its head. A differential wheeled robot is a nonholonomic robot, since the possible velocities are constrained for each value of ; e.g., it cannot move laterally along its axle. It should be noted that in this paper we also use poses in order to characterize the position and orientation of people surrounding the robot.  The possible trajectories for a differential wheeled robot can be expressed as circular trajectories passing through the center of the robot, as is shown in Fig. 1. Each trajectory is centered in a point called Instantaneous Center of Curvature (ICC). If the ICC is located at the center of the robot, the corresponding motion produces a contrary rotation of the wheels, which is an in-place rotation. On the other hand, if the ICC is located at the infinity along the robot perpendicular axis, then the wheels move at the same velocity, resulting in a straight-line motion.
The motion commands used for differential wheeled robots, i.e., the commands generated by the local planner, are pairs of linear and angular velocities ( , ), where is the desired linear velocity of the robot (i.e., the magnitude of its instant velocity vector) and is the desired angular velocity of the robot over the ICC. Therefore, the PID controller computes the individual linear velocities and for each wheel in order to fit and by applying the kinematic equations for differential drive [27]. Finally, the PID powers the motors in an open loop to maintain a constant velocity until the next motion command is applied.

The dynamic window approach for obstacle avoidance algorithm
Given the pose of a differential wheeled robot, its current instant velocity vector and the map of surrounding obstacles, there are three types of possible circular trajectories: 1. The set of circular trajectories that cannot be executed in a time due to the robot velocity vector and its hardware acceleration limits. For example: A robot moving forward cannot begin to move backwards in a very short period of time.
We consider these trajectories as non-reachable. 2. The set of reachable trajectories that can provoke a collision with surrounding obstacles, i.e., when there are one or more obstacles in the trajectory curve and the robot cannot stop in a safe manner within time due to its velocity vector and hardware acceleration limits. We consider these trajectories as non-admissible. 3. The remaining circular trajectories are called safe trajectories. They can be evaluated according to some fitness function related to the robot's goal.
In [21], the problem of computing and selecting safe circular trajectories is addressed and the authors propose the Dynamic Window Approach for Obstacle Avoidance Algorithm (DWA), which is shown in Algorithm 1.
The algorithm uses a set of predefined circular trajectories , a.k.a. motion commands. It should be noted that a motion command ( , ) ≡ (0, 0) to stop the robot must be included in the set. Then, the algorithm assigns a fitness value for each safe circular trajectory and selects a trajectory with the minimum fitness value. The fitness function is user-defined; in [21], the authors use a weighted sum of features, such as the distance to obstacles in the trajectory, the velocity of the robot and the distance to the goal position.
Furthermore, in [21], the authors give the following definition for the subset of reachable trajectories in time : where ( , ) are the current robot linear and angular velocities and (̇,̇) are the robot linear and angular acceleration limits.
Set the fitness value ← ( , , , , ) end if end for Execute a circular trajectory ( , ) ∈  with minimum fitness value.
Additionally, the authors provide the following definition for the subset of admissible trajectories: is the Euclidean distance to the closest obstacle in the circular trajectory.

The Social Force Model
The Social Force Model (SFM) was introduced in [19] and widely extended in [28,29], among other works. The SFM studies the motion of pedestrians as particles acting under the influence of forces in the environment. These forces can be attractive or repulsive, being produced by people's intentions, static obstacles (walls, furniture, etc.) and dynamic obstacles (people moving in the environment). Attractive forces usually represent people's intentions to reach specific goals. On the other hand, repulsive forces model comfort distances between people in movement, as well as safety distances to obstacles during the execution of trajectories.
In this paper, we consider a Social Force Model based on [28,29], which can be summarized as follows: The instant velocity ( ) and position ( ) at time for a pedestrian (the subject hereafter) trying to reach a desired goal as socially as possible are given by: ( ) = ( − ) + ⋅ ( ) and ( ) = ( − ) + ⋅ ( ), where ( ) is the instant force for the subject at time . Considering unitary mass, such a force is equal to the acceleration. Following the SFM, ( ) is: where 1 , 2 and 3 are the weights for the corresponding forces, and: is the subject's attractive force to at time , where is a parameter called the comfort subject's velocity, is a parameter known as relaxation time and is the unitary vector from ( − ) to . • ( ) = ∑ ( ) is the subject's repulsion force at time produced by the surrounding pedestrians, where is the number of pedestrians and ( ) is the repulsion force from pedestrian at time . The latter is defined as , where: , , ′ and are parameters; is the distance between the subject and pedestrian at time − ; is the interaction direction given by

Enzymatic numerical P systems
Membrane systems or, simply, P systems are computing models studied in Membrane Computing. They have an associated structure (given by a rooted tree or a directed graph), whose nodes (unit processors) are called compartments (membranes, cells or neurons). Compartments basically work with two types of elements: (a) multisets of symbols of a given (working) alphabet (called objects) and (b) some kinds of rules, generally inspired by nature (chemical reactions, neurons, dynamics of real ecosystems, etc.). Starting from an initial situation, the multisets of objects evolve according to fixed semantics for the rules. Membrane systems basically work with natural numbers through the multiplicity of objects in different compartments. Numerical P systems (NPS) are a special type of membrane systems where the concept of multisets of objects is replaced by numerical variables (real numbers), while the evolution rules are, in this case, programs. Numerical variables and programs are associated with compartments in a way that, from an initial situation (initial values associated with variables), they evolve dynamically, according to fixed semantics for programs.
There are membrane systems where the application of the evolution rules is controlled by certain elements, such as P systems with promoters and inhibitors. Inspired by this fact, Enzymatic Numerical P systems (ENPS) were introduced in [11] as an extension of numerical P systems in which some variables, named enzymes, control the execution of the programs. This computing model is being applied in different areas, especially in the simulation of control mechanisms of mobile and autonomous robots.
is an alphabet of labels, containing symbols; 2. is a rooted tree with nodes bijectively labeled by elements from ; 3.
ℎ , ℎ ∈ , is a finite set of numerical variables ,ℎ associated with the membrane labeled by ℎ, ℎ , ℎ ∈ , is a finite set of programs associated with the membrane labeled by ℎ, where each program ∈ ℎ has the following form } using the set of enzymes ℎ = { 1,ℎ , … , ℎ ,ℎ }, in such a way that it disables the program when the output is . the right hand side ,1 | ,ℎ , … , , | , is an expression, with: ,1 , … , , being natural numbers, and ,1 , … , , numerical variables from membrane ℎ, the parent membrane of ℎ and all children of ℎ according to .
An enzymatic numerical P system, = ( , , {( ℎ , ℎ , ℎ , ℎ (0)) | ℎ ∈ }), of degree ≥ 1, can be viewed as a set of membranes, labeled by elements of , arranged in a hierarchical structure given by a rooted tree, whose root is called the skin membrane, thus the membrane labeled by ℎ has: (a) a finite set ℎ of numerical variables ,ℎ (their initial values are defined by the set ℎ (0)), (b) a subset ℎ of variables from ℎ , called enzymes, and (c) a finite set ℎ of programs, where each program has the following form A configuration or instantaneous description of an ENPS at time ∈ N is a tuple ( 1 ( ), … , ( )), where ℎ ( ), for ℎ ∈ , provides the values of all numerical variables ,ℎ from the set ℎ at time (the value of ,ℎ at time ∈ is denoted by ,ℎ ( )). The initial configuration of is the tuple ( 1 (0), … , (0)).
is . When such a program is executed to a membrane labeled by ℎ from configuration  : (a) the value of function at time is computed: . Then the value , ( + 1) of the variable , at time + 1 will be the sum of total contributions that , receives from the neighboring membranes at time . It is worth pointing out that the enzymes ,ℎ with 1 ≤ ≤ ℎ , which enables the execution of programs, are variables from ℎ . Therefore, their values can change due to the contribution they receive from other programs and compartments.
Let us assume that there exists a global clock marking the time. In our semantics, all the programs that can be executed are selected in each time unit. The concept of transition step from a configuration  to  +1 and the concept of computation for ENPS are defined in [11].

A social local planning algorithm for differential wheeled robots
In this section, we propose a social local planning algorithm for differential wheeled robots, using the SFM, presented in Section 2.4, to compute the robot's optimal instant velocity vector for each time step.
In an actual environment, not all the velocity vectors can be executed by a robot due to acceleration limits and collision threats, as explained in Section 2.3. Moreover, some velocity vectors cannot be executed due to nonholonomic constraints, as explained in Section 2.2. Therefore, the problem is to find a circular trajectory, i.e., a motion command, approaching the optimal instant velocity vector. In our approximation, we apply the DWA algorithm to solve it. Thus, for each time step, we use a fixed set of motion commands and filter admissible and reachable commands with respect to acceleration limits and surrounding static/dynamic obstacles. Finally, a fitness function is computed for each filtered command, and one with the best fitness value is executed by the PID controller. The fitness function evaluates the difference between the corresponding circular trajectory and the optimal instant velocity vector, as explained below in this section. The computation is repeated in a loop until the robot's goal is reached or an error is produced, as shown in Algorithm 2. where 1 and 2 are parameters. More precisely: 1 is the weight for the error in the robot's position and 2 is the weight for the error in the robot's yaw angle. Such values can be set respectively to 1.0 and 0.1 by default and optimized or tuned by simulation or experimentation. Finally, it should be noted that must be low enough for a safe navigation in real time, e.g., 0.01 or 0.02 s.
The set of initial values is as follows: where the word input means a value obtained by the robot's on-board sensors and the word user means an user-defined value. The set of Enzymes is 1 ≡ { , 0 } and the set of programs 1 is as follows: In the first step of the computation, programs 1 and 2 are executed to generate the robot's force to the goal position. Programs 3, ∶ 1 ≤ ≤ and 4, ∶ 1 ≤ ≤ are executed to generate, respectively, the Euclidean distances from the robot to each pedestrian position and from the robot to each obstacle position.
In the second step of computation, programs 5, , 6, ∶ 1 ≤ ≤ generate the unitary vectors from the robot's position to each pedestrian position. Programs 7, , 8, ∶ 1 ≤ ≤ produce the corresponding unitary vectors from the robot's position to each obstacle position.

Simulator
A simulation tool was implemented in this work for validation and testing purposes. The Robot Operating System (ROS) framework was selected to develop the simulator. ROS is a flexible framework for writing robot software. It contains a wide variety of tools and libraries, which simplifies the task of creating robust software. ROS can be used together with Python and/or C++. In this work, the C++ programming language was selected. We adapted the BSD-licensed libraries LightSFM 3 and Pedlab 4 to simulate the human agents. The RVIZ tool 5 was used to generate a real-time 3D visualization. The simulator is GNU GPLv3 licensed and can be downloaded from https://github.com/RGNC/enps_dwa.   The software simulates an indoor environment with obstacles containing a configurable number of human agents walking according to a social random walk, i.e., each agent selects a nearby goal and walks approaching it by following the SFM explained in Section 2.4, when the goal is reached, another goal is selected. A dual-wheeled non-holonomic robot agent is also included, navigating according to the ENPS model proposed in Section 4. The robot's goals are randomly selected in the same way.
The ENPS model was simulated ad-hoc, i.e., the model definition was developed inside the simulator's source code. More precisely, the file simulator.hpp implements the behavior of the model.
A 3D visualization of the virtual world is generated in real-time by using RVIZ. A screenshot of the visualization is shown in Fig. 2. The visualization includes the environment, the agents, the force vectors, the simulated LIDAR scan, the people detections, the possible robot's motion commands, i.e., circular trajectories and the selected robot's trajectory at each time-instant.
The robot model is based on the TERESA robot, whose parameters are: radius of 0.3 m, linear velocity from 0 m/s to 0.6 m/s, angular velocity from −0.8 rad/s to 0.8 rad/s, people detection range (360 • ) of 2.5 m and obstacle detection range (360 • ) of 10 m. A more thorough explanation can be found in [30].
Regarding the computational complexity of the simulator, it is worth noting that it is greater than ( ( )), since the implementation is purely sequential for the sake of simplicity and according to the purposes of validation and testing.

Analysis
We run an experiment of 160 min by using the software presented in Section 5. Twenty human agents and 1 robot agent were included in the simulation. The robot and SFM parameters were set to the common parameters of the TERESA robot. A set of 33 possible angular velocities from −0.8 rad∕s to 0.8 rad∕s was combined with a set of 13 possible linear velocities from 0 m∕s to 0.6 m∕s, producing a whole set of 429 possible motion commands for the robot. The initial velocity of each human agent is randomly set by using a Gaussian distribution with = 1.2 m∕s and = 0.001.
For each agent and for each second of simulation, the software records in a file the instant velocity and the Euclidean distances to the corresponding nearest obstacle and nearest agent. Despite the natural differences between people and robot velocities shown Fig. 4. The distances from the robot and from each human to their corresponding nearest obstacle and nearest agent were measured at each second of simulation. Such distances can be used as metric to compare the social behavior of the agents. Despite the differences in the velocity PDFs, a high similarity is observed in the average distances.
in Fig. 3, we observed a similar social behavior when considering the average distances to obstacles and other agents during the simulation, as shown in Fig. 4.

Conclusions and future work
This work presents the first model based on membrane computing for the social navigation of robots in human environments. Such a model is based on the ENPS framework, the Social Force Model and the Dynamic Window Approach Algorithm. It is carefully designed for differential-drive wheeled robots where nonholonomic and kinematic constraints are taken into account. The computational complexity of the model is ( ( )) and it is compatible with other ENPS models for the navigation of robots. The main advantage of using the ENPS framework is the natural way to expose parallelism in the algorithm and, therefore, enabling the possibility of using parallel hardware simulators. A simulation tool was implemented in this work for validation and testing purposes. The simulation was conducted in a 3D environment and metrics about velocities and distances were recorded. We run experiments with the software in order to compare the social behavior of the simulated human agents and the simulated robot. Despite the differences related to non-holonomic constraints and velocity limits, we observed a similar social behavior in the simulation with respect to average distances to nearby obstacles and agents. The main research line for future work is the simulation of the model by using parallel architectures in order to approach the ( ( )) computational complexity, as well as the integration of this model in the navigation workflow of actual robots.

Declaration of competing interest
No author associated with this paper has disclosed any potential or pertinent conflicts which may be perceived to have impending conflict with this work. For full disclosure statements refer to https://doi.org/10.1016/j.compeleceng.2021.107408.