To examine the accuracy of the NMPC predictions, a 27 year old male subject was selected to perform reaching tasks. An Optotrak Certus motion capture system (Northern Digital Inc., NDI) was used to measure the arm trajectory at 30 Hz. In these experiments, the subject was seated with the arm elevated at the shoulder level. An active marker attached to the back of the hand (as shown in Figure 1A) has been used to capture center-out hand motion trajectories. The subject was asked to move his hand from an initial central position to one of eight final targets spread evenly on a circle of 20 cm radius at a self-selected convenient speed. The experiment was repeated 10 times for each target with 2 min rest intervals between each set. The subject also performed reaching to moving targets. He was instructed to reach to a target, which was relocated to another position midway through the movement. The subject was instructed to adjust his motion to reach to the moving target. The subject was also given 5 min to rest before performing the reaching moving targets experiment.

During the reaching trails, the electromyography (EMG) activity of seven muscles (anterior/middle/posterior deltoid, long/lateral triceps brachii, biceps brachii, and brachioradialis) were collected at 2000 Hz using a Trigno portable EMG system (Delsys Inc.). The EMGs were band-pass filtered (5–800 Hz cut-off), rectified, low-pass filtered (4 Hz cut-off), and normalized to maximum voluntary contractions (MVCs). We performed a Pearson correlation analysis to investigate the correlation between the EMGs and muscle activity predicted by NMPC and DO simulations. We resampled both the captured EMGs and simulation results with a sampling rate of 100 Hz, and performed the Pearson correlation analysis using the “Corr” command in MATLAB. The experiments have been approved by the Office of Research Ethics at the University of Waterloo and carried out with written informed consent from the subject. The subject gave written informed consent in accordance with the Declaration of Helsinki.

In this research, a planar arm model similar to the one developed by Ghannadi et al. (2015) was used. The model consisted of torso, upper arm, and forearm to simulate the hand motion. The torso was fixed and the shoulder and elbow were modeled using revolute joints. Six muscle groups including shoulder and elbow mono-articular flexors/extensors and two bi-articular flexors/extensors were used to actuate the arm as shown in Figure 1B. A modified Hill-type muscle model with muscle excitation-to-activation dynamics was used to simulate the skeletal muscle contraction dynamics (see Appendix A for details). The muscle parameters of the planar arm model (i.e., insertion and origin positions, maximum isometric force, fiber optimal length, slack length, and pennation angle) were tuned to represent the dynamics of the upper extremity in the experimental condition (reaching targets in a horizontal plane elevated at the shoulder level). These parameters were tuned through a series of optimizations so that the planar model provides the same joint torques as a high-fidelity three-dimensional upper extremity model (Ghannadi et al., 2015). Kinematic and dynamic parameters of the arm model and the Hill muscle model parameters used here can be found in Tables 1, 2, respectively.

The arm motion is controlled by complex commands descending from the CNS, which are the combination of the motion prediction (feedforward control) from an internal representation of the body and environment (or so-called internal model; Desmurget and Grafton, 2000), and the corrective command from the sensory organs to correct any errors due to uncertainty or unknown environment (feedback control). This complexity is captured here by a model-based NMPC with a receding horizon. The NMPC uses a control-oriented model (COM) representing the human's internal model to predict the optimal trajectory, and feedback information to correct the prediction errors. The NMPC predicts the optimal dynamics of the system (x, u) over a prediction horizon as shown in Figure 2A by minimizing the following cost function: (1)J=Ψ (t0+tph)+∫t0t0+tphψ (x(t),u(t)) dt (2)subject to: 0<u(t)<1 where Ψ is the cost evaluated at the end of prediction horizon, ψ is the cost evaluated during the prediction horizon, and t_ph is the length of prediction horizon. The state variables at the current time (t_o) are obtained from the current sensory information. The input (ū) is an optimal open-loop solution over the prediction horizon. If there are no external disturbances and no model uncertainty in the system, with infinitely long prediction horizon, the open-loop solution can be applied to the system for all time t > t_o. However, for the finite horizon case and in the presence of noise and uncertainty, the open-loop solution should only be applied until the next sampling time (t₀ + δ). At the new time step, the optimal solution is re-evaluated with the new initial conditions for the receding horizon and iteratively applied to the system. By incorporating the feedback information, the NMPC is converted from a completely open-loop controller to an optimal closed-loop controller. The NMPC can handle constraints on both the states and the inputs. In musculoskeletal models, the muscle activation command must be non-negative and less than one, and constraints on states can be added to avoid unphysiological movements.

In this article, the optimal dynamics over the prediction horizon were calculated using the GPOPS-II optimal control package that utilizes an orthogonal collocation method (Patterson and Rao, 2014). This method is a direct (simultaneous) optimization method in which both states (x) and inputs (u) are parameterized using a series of connected Legendre polynomials and become part of a Nonlinear Programing (NLP) problem. Here, the arm model described in the previous section was used as the COM of the NMPC and the simulation model for planar reaching and pointing tasks. The dynamic equation of the arm model can be described by: (3)x˙=f(x(t),u(t)),    x(0)=x0 where x ϵ ℝ^{10 × 1} are the arm model state variables consisting of shoulder angle and angular velocity plus elbow angle and angular velocity, and muscle activation states, and x₀ is the vector of the initial states. The muscle excitation inputs u ϵ ℝ^{6 × 1} represent the ratio of excited motor units to the maximum number of motor units in that muscle.

In this research, for the particular case of a goal-directed reaching task, the terminal cost of the NMPC cost function (Ψ) was removed, and the integral part (ψ) is computed from the summation of two terms: (i) a choice of specific physiological cost function, and (ii) a trajectory tracking error. Therefore, the NMPC cost function shown in (1) is converted to: (4)J=∫t0t0+tph(p(ζ(t)-ζdes)2+q GM(u (t))) dt  where p and q are cost function weightings, and ζ and ζ_des are the hand position and its desired final value (in the Cartesian coordinate system), respectively. The simulated hand position ζ varies on the prediction horizon, while the desired final value ζ_des was kept constant. The hand position is calculated from.

where θ₁ and θ₂ are shoulder and elbow angles, and L₁ and L₂ are upper arm and forearm lengths, respectively. The physiological cost G^M(u(t)) is defined as:

The term G^M in the cost function represents the neural excitation effort to perform the reaching tasks.

GPOPS-II finds the optimal dynamics of each given horizon by minimizing Equation (4) while satisfying inequality constraints related to muscle excitation (2) and equations of motion (3) using the Sparse Nonlinear Optimization (SNOPT) solver (Gill et al., 2005). An hp-adaptive mesh refinement method (Liu et al., 2015) has been used within GPOPS-II to refine the individual interval widths and the polynomial degree to reach a final optimal solution. Then, the first five-percent (e.g., 50 for 1000 ms prediction horizon) of optimal activations are applied to the muscles, and the arm motion is simulated. The new position and orientation of the arm are measured and sent back to the NMPC as initial conditions of the next iteration. In this research, we have assumed that the sensory organs can measure the exact joint angles and angular velocities. Uncertainty can be added to the measurements to account for the noise within the sensory organs, and to simulate the variability in the movement repetition. However, this has not been included in the scope of this work. The optimal muscle activations are shifted and considered as the initial guess of the next iteration.

In addition to the NMPC simulations, a dynamic optimization (DO) using GPOPS-II was performed to simulate the same task. Unlike the NMPC simulations, which continue until the position tracking error passes a certain threshold, the final simulation time and the final position of the hand are explicitly specified in the DO simulations. The DO cost function is:

where t_f is the final simulation time, and x_f is the state vector corresponding to the target position. The same physiological cost function as in the NMPC simulations (Equation 6) was used in DO to compute the optimal hand position trajectory and muscle activations.

In this section, a goal-directed reaching task is simulated and the effect of prediction horizon length variation on the hand trajectory and muscle activation is studied. Here, the hand is initially at rest and in a natural position (θ1=44o and θ2=58o) and moves toward a target 20 cm to the left of its initial position as shown in Figure 2B. This task was simulated using NMPC with 0.2, 0.3, 0.4, 0.5, and 0.8 s prediction horizons, and DO with a fixed time duration. In the NMPC simulations, the cost function weightings (p = 20 and q = 1) were kept the same. The DO final simulation time was chosen to be 1.5 s, in accordance with the experimental reach duration (1.431 ± 0.176 s in 10 repetitions).

Figure 2B demonstrates the hand trajectories for the aforementioned prediction horizons (t_ph). Despite the slight differences between the trajectories, the solutions with 0.4, 0.5, and 0.8 s prediction horizons closely correlate with the experimental results. The Pearson correlation factors between NMPC simulations and the experimental trajectories are 0.9351, 0.9758, 0.9947, 0.9991, 0.9994 respectively for 0.2, 0.3, 0.4, 0.5, 0.8 s. The Pearson factor between the DO and experiment is 0.867.

As shown in Figure 3, the elbow and shoulder angle variations with 0.2 s prediction horizon are less than those with longer prediction horizons. This signifies the importance of prediction horizon; with short prediction horizons, the controller is more cautious and takes longer to reach a desired position. Simultaneously, this results in smaller muscle activations for shorter horizons since the transient time to reach the final position is longer (see Figure 4). The reaching error at 1.5 s of the simulation is about 23% for 0.2 s prediction horizon, and reduces to 0.36% for 0.8 s prediction horizon. As expected for this motion, the shoulder mono-articular flexor (muscle 1), elbow mono-articular flexor (muscle 3) and bi-articular flexor (muscle 5) are activated at the beginning to accelerate the body; then, the antagonistic muscles are activated to reach a full stop at the desired position.

In the DO, similar to NMPC, flexor muscles are active at the beginning of the motion to accelerate the hand toward the target, then the extensor muscles are activated to stop the hand movement (a bang-bang control strategy). Since DO has to stop at the specified final time (1.5 s) the extensor muscle activities are larger than NMPC predictions at the decelerating phase of motion. On the other hand, the flexor muscle activities at the accelerating phase of the simulation are larger for NMPC than DO because the trajectory error at the beginning of the motion is large and exponentially reducing when it gets closer to target.

As shown in Figure 4, the muscle activations predicted by NMPC and DO simulations can capture the general trends of the experimental measurements. It can be observed that as the prediction horizon increases, the NMPC predicts larger muscle activities at the beginning of the motion. A Pearson correlation analysis was performed between the flexor muscle activities predicted by the NMPC and DO simulations and the EMGs from experimental measurements; the correlation coefficients are presented in Table 3. The correlation coefficient for extensor muscles are not reported since the EMG activity of these muscles were minimal in the experiments. As shown in Table 3, the correlation coefficient of flexor muscles reduces when the prediction horizon increases in the NMPC simulations, while the DO predictions correlate better with the experiments.

As suggested by Morasso (1981), subjects tend to move in straight lines with bell-shaped tangential-velocity profiles when reaching for a target. Figure 2B shows that the NMPC with long-enough prediction horizon can realistically predict the hand trajectory, while Figure 5 shows the differences in the velocity profiles of the NMPC and DO predictions. In these simulations, the NMPC tends to accelerate the hand more quickly than the DO; in DO, the optimization knows the final time and can distribute the acceleration over a longer time. On the contrary, in the receding horizon NMPC, the controller can only predict the motion as far as the prediction horizon. This results in the fast acceleration at the beginning the motion due to large tracking errors and slow deceleration at the end of motion due to small tracking errors.

Reaching to eight different directions was also simulated using NMPC with 0.8 s prediction horizon. The target positions of the hand are located on a circle centered at the initial position of the hand with a radius of 20 cm (see Figure 1). In the NMPC simulations, a smooth 5th-order polynomial with zero initial and final velocities and accelerations is used as the desired straight-line hand trajectory. These trajectories begin at the initial position of the hand and end at the target positions. As shown in Figure 6, the NMPC is able to follow the desired trajectories, which are qualitatively correlated with the measured trajectories in experiments.

We have assumed that the CNS plans a trajectory to reach a target and constantly monitors the deviations from this trajectory and the target position. In this section, we study the case where the target position is suddenly relocated. Here, the hand is initially at rest at point O (Figure 7A) and moves toward the target at point A. Then 1 s later, the target position suddenly moves to the point B. This protocol was achieved in the lab by manually moving the target from its initial point A to B when the subject reached half the way to A. In this simulation, the time delay related to the visual cognition of this change [about 150 ms, (Jeannerod, 2006)] has not been considered.

Figure 7 depicts that the NMPC controller can track the location of the target and correct the hand trajectory to reach the new target. As shown in Figure 7A, it seems that when the target moves, the subject over-compensates by moving the hand to the right, while the NMPC finds a trajectory that minimizes both position error and control effort. It is not possible to simulate this scenario with DO; it is one of its disadvantages compared to NMPC, which is able to make online adjustments to the hand trajectory.

The focus of this paper has been on the proof-of-concept of the NMPC as a possible model for CNS control of human movement. The current implementation of this approach is computationally expensive and is not real-time. For instance, with a prediction horizon of 0.5 s, the NMPC takes 0.45 ± 0.24 s to find the optimal dynamics at each time step and re-plan the movement. These simulations were performed on a computer with an Intel Core™ i7-4790 processor and CPU 3.60 GHZ and RAM 16 GB. However, online NMPC methods such as the Continuation/GMRES method (Ohtsuka, 2004), advanced-step NMPC (Zavala and Biegler, 2009), and explicit MPC (Kouramas et al., 2013) can be used to achieve real-time performance. As an example, Mehrabi et al. (2016) developed a Newton/GMRES NMPC controller to control the functional electrical stimulation of knee extension. This controller by discretising the system dynamics and employing a fast online optimization method (GMRES), significantly reduced the computational time and allowed the real-time implementation of the NMPC. Furthermore, the COM can be further simplified using muscle synergy theory, in which the CNS coordinates human body movements by bundling individual muscles into groups. This allows a low-dimensional control input that significantly reduces of the size of the control problem (Sharif Razavian et al., 2015).

In this research, the first use of NMPC to simulate human motor control in reaching movements was presented. It was shown that NMPC can replicate certain properties of the human motor control system (i.e., path-planning, prediction, and target tracking), and can be used to realistically simulate reaching movements. Due to its feedback nature, it can correct the tracking errors for static targets or can follow a moving target seamlessly. The NMPC prediction horizon can represent the time horizon for which the CNS minimizes a physiological cost function. It should be noted that the NMPC conclusions from this research are specific to the cost function used here; stronger conclusions can only be made if more diverse cost functions are investigated. Nonetheless, this method opens up new opportunities to study challenging problems such as predictive forward dynamic simulation of biomechanics and biomechatronic systems.

As a possible future research direction, an online NMPC can be used to represent a user/patient in an assistive devices controller to facilitate the shared control between the device and user. This shared control allows the device to perform some tasks independently of the user by sensing information about the environment (Millán et al., 2010). By predicting the motion of the user and adjusting the trajectory online, the NMPC can reduce the cognitive workload imposed on the user, who does not need to consider low-level executions in the presence of external disturbances or obstacles (Tucker et al., 2015). The variability in limb movement is another known characteristic of reaching movements. This characteristic can be incorporated in the NMPC assistive device by accounting for noisy sensory information and sending noisy motor commands to musculotendon units.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.