The robotic arm we used in this research comprises nine servo actuators (7 × Dynamixel's XM540-W270, 2 × Dynamixel's XM430-W350). All actuators are capable of 40 N radial load and have an embedded Cortex-M3 microcontroller. The M3 is coupled with contactless 12 bit absolute encoders, allowing the retrieval of the actuator's position, velocity, and trajectory as feedback for position estimation. The XM540 actuators are used for arm movements and have a stall torque of 10.6 Nm (at 12 v input). The XM430 actuators are used to manipulate the EE (grasping, rotating) and have a stall torque of 4.1 Nm (at 12 v input). Actuators are manufactured by ROBOTIS (Korea). Actuators do not have torque sensors and were therefore actuated by current specifications. The relation between the driven current and the generated rotational velocity is not linear as it has to account for friction. In this work, the current-speed association was estimated as described below. Arm chassis is based on 3D printed grippers (allowing function-tailored customization), ridged and lightweight T-slot extruded aluminum arms, and aluminum brackets. The chassis is connected to the actuators with industrial-grade slewing bearings, and it was assembled by Interbotix (Downers Grove, Illinois). Motion control was evaluated on the Nvidia Xavier chip (Jetson AGX Xavier) and then realized on Intel's Loihi chip. Communication with the daisy-chained servos was based on TTL half-duplex asynchronous serial communication, handled by Dynamixel's U2D2 control board. Overall, the arm design provides 6 DOF, 82 cm reach, 1.64 m span, 1 mm accuracy, and 750 g payload.

To simulate the robot described above, we used the Multi-Joint dynamics with Contact (MuJoCo) physics simulation framework. The robotic arm and joint's accurate dynamic were specified using CAD-derived mechanical description and inertia and mass matrices. CAD and dynamic specifications were provided by Trossen Robotics (USA). The simulation was developed using Nengo, a Python package for building, testing, and deploying NEF-based neural networks.

FK transform a robot's configuration to the cartesian coordinates of its EE. Here, it was implemented using transformation matrices, which characterize the relative transformation (rotation, translation) from each joint to the next. For our five joints robot, FK would take the form of T = T₀₁T₁₂T₂₃T₃₄T₄₅T₅₆, where T_ij is the transformation matrix in homogenous coordinates from the reference frame at joint i to reference frame at joint j (indices 0 and 6 refer to the world's and EE coordinate reference frames, respectively). Initializing T with the appropriate set of rotations and translations and multiplying it by a zero vector [0, 0, 0, 1]^T (in homogenous coordinate) will result in our FK model T_x(q):

Where s_x = sin(q_x), c_x = cos(q_x), and q_x is actuator x angle of rotation. T_x(q) returns the EE position in the world's coordinate system, where the origin is at the robot's base. The numerical coefficients were derived by calculating the transformations while taking into account the robot's geometry, retrieved from the robot's CAD file.

IK refers to the transformation in which a robot's configuration is computed from its EE's desired location. Generally, IK cannot be analytically solved, and it is, therefore, usually numerically optimized. Here we used the Jacobian inverse for IK. We calculate the Jacobian J of T, which relates the change of the EE position x to the change of joint angles q: J(q)=δxδq. The Jacobian relates a change in robot configuration q. to a change in EE position ẋ with:

Equation (2) allows us to specify a target in task space—that is, the cartesian space centered on the EE's origin—rather than the space that can be directly controlled, the configuration space. Note that the Jacobian has to be recalculated along the trajectory. With our robotic system, the calculated Jacobian has the shape of (3, 5), where 3 is the number of space dimensions (task space) and 5 is the number of joints (configuration space). To compute IK, we need to invert Equation (2). Since the Jacobian is not necessarily invertible, a common practice is to use its pseudo-inverse form J⁺ allowing to compute q.=J+ẋ. Therefore, given an error in space coordinates x_d as the difference between the EE current position x_c and its target position x_y, the appropriate change in joint space can be computed using:

Where d (q) is the change in joint angles, for which the robot's EE will get closer to its target. In each iteration, this equation is re-evaluated until x_d is within some accuracy threshold. Once the joint configuration for a given target point is concluded, control signals, which achieve it, can be calculated using PID control. Further details are provided in Lynch (2017).

PID control is used universally in applications requiring accurate control (Ang et al., 2005). Given a target position, a PID controller will continuously reduce an error signal by providing the robot's actuator with the appropriate control signal to come closer to its target. To do so, the PID controller generates a signal u (t), which is proportional to the value of the error signal e (t) (accounting for the current value of the error), to the error integrated value over time (accounting for the past values of the error), and to the error derivative (accounting for the projected value of the error), using:

Where K_p, K_i, and K_d are the proportional, integral, and derivative gain coefficients, respectively.

In this work, we've implemented IK and PID on Intel's neuromorphic research chip Loihi (Davies et al., 2018). NEF was compiled on the board using the nengo_loihi library (version 0.19) (Lin et al., 2018). The nengo_loihi library was designed to execute Nengo models on Loihi boards. It contains a Loihi emulator backend for rapid model development and a hardware backend for running models on the board itself. Nengo Loihi's hardware backend uses Intel's NxSDK API to interact with the host and configure the board. Each Loihi chip is comprised of 128 neuron-cores; each simulates 1,024 neurons and has 4,096 ports. Each chip also has × 86 cores, which are used for spike routing and monitoring. Communication was established via an SSH channel between our local computer and a virtual machine installed on Intel's neuromorphic research cloud.

IK was implemented with NEF using Nengo and tested on our robotic arm. Our model schematic is shown in Figure 1A. In terms of joint angles, the configuration of the robot was introduced through a node into a neuron ensemble denoted Current q. Ensemble Current q is fully connected to another neuron ensemble, denoted Target q. These synaptic weights are modulated to minimize the value decoded from neuron ensemble Error using PES optimization. Ensemble Distance to target encodes the EE distance from its designated target by subtracting the current position of the EE (calculated using Equation 1, applied on the value decoded from Target q) with its desired position (given as an input through node xyz Target). Both Target_q and the distance to target ensembles are connected to a 5D ensemble, allowing for non-linear computation between the five robot's joint states. The difference between the current and the desired robot's joint configuration is calculated through Equation (3), which is implemented through the connection to ensemble Error. When the Error decoded value is 0, ensemble Target q encodes the desired robot configuration. The Error ensemble is connected to an inhibition signaling node. Upon actuation (initiated once a sufficiently accurate result is achieved), the error signal is zeroed; thus, Target q is stabilized at its current state. Here, we performed IK on our robot's 5 joints, starting from an initial configuration where all joints were zeroed (Figure 1B). As learning progresses, the new robot configuration is calculated (Figure 1C), while the error is continually minimized (Figure 1D). Raster plots of ensembles q, target q, and error are shown in Figure 1E. While the spiking pattern, which represents the initial joint configuration, is constant, the target's spiking pattern changes as the error spiking pattern becomes more amorphous, indicating convergence to zero.

We further analyzed the model by modulating neurons' encoders and learning rates. Each neuron's intercept defines the part of the representation space in which the neuron responds by firing, and it is reflected on the neuron tuning curve (firing rate as a function of input). Note that the intercept is the input value for which the neuron initiates spiking at a high rate. Distributing intercepts uniformly between −1 and 1 makes sense for 1D ensembles for which they create a uniform spanning of that space. This is not the case for high dimensional ensembles. In our implementation, we used high dimensional ensembles to represent the DOF of our robotic system. With uniformly distributed intercepts, the resulted tuning curves are uniformly distributed, resulting in a non-efficient spanning of the representation space. As a result, the system does not converge to its target, as is evident from the error's non-decreasing value (Figure 2A). Changing the intercept distribution to follow a triangular pattern modulates the neurons' tuning curved distribution such that the representation space is adequality spanned. As a result, the system converges to its target, as is evident from the decreasing error (Figure 2B) (Gosmann and Eliasmith, 2016). This modification of the intercept distribution is crucial for accurate representation in 5D space and it is briefly described in DeWolf et al. (2020) and Gosmann and Eliasmith (2016). Our model relies on PES-based optimization, and it is therefore constrained to a prespecified learning rate. We tested our model with three different learning rates, and as expected, error flattening is slower as we increase the learning rate (Figure 2C). Suppose we permit our system to keep optimizing. In that case, the ensembles' fluctuating encoded values induce continuously changing results, where changes in one joint's angle are compensated with changes in other joints (convergence is driven toward a zero gradient potential field). Therefore, we are inhibiting learning once some accuracy threshold is reached, thus holding the computed weights constant and providing a stable target configuration (Figure 2D). Once the desired configuration is calculated, the robot should be actuated accordingly by using, for example, a PID-controller.

A PID controller integrates three error modulations to provide the desired actuation, such that the system would approach a target position. These three signals are described in Equation (4) and are schematically presented in Figure 3A. Our PID controller holds a model of engine actuation. Here, we modeled the actuators using a basic speed-torque (implemented with a driving current) model, which corresponds to our physical actuators. In our model, the actuator experiences static friction and it responds exponentially fast once its gears become active, saturating at some maximum speed. As we stop driving the engine, it is losing momentum due to dynamic friction. When actuation is reversed (current in reversed direction), the position is changed accordingly. The actuation model is shown in Figure 3B. In this work, we implemented the PID with spiking neurons using NEF. The model schematic is shown in Figure 3C. The robot's current configuration is introduced through node Current q. We subtract a feedback signal y (t) from it to compute an error signal. This error signal is propagated to the output ensemble through three paths: 1. Proportional path in which the error is proportionally transformed through a gain factor k_p, producing signal e_p(t); 2. Integration path in which the error is integrated using a neuromorphic integrator (see Methods for further details). The result is scaled by a gain factor k_i, producing signal e_i (t); and 3. A derivative path, implemented by connecting the error ensemble to a 2D derivative ensemble. To implement derivation, the error is propagated through two synapses: one with a short time scale (τ) and the other with a longer one. The two values are subtracted and scaled by a gain factor k_d, producing signal e_d (t). These error signals e_p(t), e_i(t), e_d(t) are summed in the output ensemble, delivered to the engine as u (t), which also feedback as y (t). A running example of this model is shown in Figure 3D. Given a target angle position for the engine (normalized to 1), the error is quickly reduced as the engine's location approaches its target. Raster plots for y (t) and u (t) are shown in Figure 3E. These stable spiking dynamics of the control signals reflect fast convergence to target.

To further analyze our neuromorphic PID control, we examined it as a P (proportional path was enabled), a PI (proportional and integrative paths were enabled), and a PD (proportional and derivative paths were enabled) controller. By implementing all three models, we demonstrated the classic PID characteristics in a neuromorphic implementation. Particularly, the P controller was shown to fall short of reaching the target, the PI controller reached the target with inefficient dynamics, and the PD controller had an improved reaching dynamic, but it failed to reach the target accurately (Figure 4A). We further examined our model by changing the number of neurons. Allocating 250 neurons per ensemble per dimension produced accurate results. Reduced number of allocated neurons dramatically affected performance and stability (Figure 4B). This result is compatible with the noise characteristics of NEF-based representation in which the decoders-induced static noise is proportional to 1N2, where N is the number of neurons (Eliasmith and Anderson, 2003). Synaptic time constants also constrain neuromorphic implementations. Reducing these time-constants inhibits the integration dynamic (Equation 8), as demonstrated in Figure 4C.

Control was evaluated on a physical 6 DOF robotic arm (described in the Methods). To demonstrate robot performance, we utilized the MuJoCo physical simulator, in which we accurately described the dynamic and mechanical characteristics of our physical arm. Here, we used it to demonstrate the integration of neuromorphic IK and PID control in a physical setting. While IK was used to derive the robot configuration from the desired EE location in space, PID was used to actuate the robot, generating an EE trajectory toward the target. We created a uniformly distributed 10,000 target points in a 2 × 2 × 1 meters volume (Figure 5A). We tested each point for reachability using IK, constructing a 3D reachability map, where a black point designates a reachable point with an accuracy of at least 1 mm (Figure 5B). Arm base is located at the origin (0, 0, 0). For demonstration, we randomly chose two points and used PID control to generate robot motion. The selected points, the final arm configuration, and the generated EE trajectories are shown in Figure 5C. Trajectories are linear (minimal path) as expected. Distance to target curve is shown in Figure 5D, demonstrating fast convergence to target.

We implemented both IK and PID control on Intel's Loihi chip. When implementing IK with different learning rates, the same error convergence pattern appears in both simulation and on the board (Figures 6A,B). However, superimposed results showed that the Loihi could converge better, as its error reduced faster than in the simulated model (Figure 6C). We found it to be consistent for various learning rates (Figure 6D). When implementing PID control on the Loihi, we found it hard to implement the derivative pathway, as it is currently cannot support high synaptic time constants. Defining a long time-constant is essential, as the generated control signal fluctuates, affecting our capacity to derive the signal's rate of change accurately. However, for a time-constant of 10 ms, and a configuration of k_p= −1, k_i = −0.1, k_d= 0.35, the Loihi was able to converge faster than simulation to the desired target, pointing out its embedded learning accelerator (Davies et al., 2018).