Effects of model inaccuracies on reaching movements with intermittent control
Igor Gindin aff001; Miri Benyamini aff001; Miriam Zacksenhouse aff001
Authors place of work:
Faculty of Mechanical Engineering, Technion Israel’s Institute of Technology, Haifa 32000, Israel
Published in the journal:
PLoS ONE 14(10)
Background and objectives
Human motor control (HMC) has been hypothesized to involve state estimation, prediction and feedback control to overcome noise, delays and disturbances. However, the nature of communication between these processes, and, in particular, whether it is continuous or intermittent, is still an open issue. Depending on the nature of communication, the resulting control is referred to as continuous control (CC) or intermittent control (IC). While standard HMC theories are based on CC, IC has been argued to be more viable since it reduces computational and communication burden and agrees better with some experimental results. However, to be a feasible model for HMC, IC has to cope well with inaccurately modeled plants, which are common in daily life, as when lifting lighter than expected loads. While IC may involve event-driven triggering, it is generally assumed that refractory mechanisms in HMC set a lower limit on the interval between triggers. Hence, we focus on periodic IC, which addresses this lower limit and also facilitates analysis.
Theoretical methods and results
Theoretical stability criteria are derived for CC and IC of inaccurately modeled linear time-invariant systems with and without delays. Considering a simple muscle-actuated hand model with inaccurately modeled load, both CC and IC remain stable over most of the investigated range, and may become unstable only when the actual load is much smaller than expected, usually smaller than the minimum set by the actual mass of the forearm and hand. Neither CC nor IC is consistently superior to the other in terms of the range of loads over which the system remains stable.
Numerical methods and results
Numerical simulations of time-delayed reaching movements are presented and analyzed to evaluate the effects of model inaccuracies when the control and observer gains are time-dependent, as is assumed to occur in HMC. Both IC and CC agree qualitatively with previously published experimental results with inaccurately modeled plants. Thus, our study suggests that IC copes well with inaccurately modeled plants and is indeed a viable model for HMC.
Control theory – Covariance – Eigenvalues – Noise reduction – Simulation and modeling – Gaussian noise
Stochastic optimal feedback control (OFC) has been strongly advocated as a framework for investigating human motor control (HMC) [1–5]. The inherent noise in HMC is handled by an optimal observer that estimates the state, while inherent delays are overcome by a predictor that predicts the current state from time-delayed state estimation [3, 6]. Standard HMC theories assume that the communication between these processes is continuous, though others argue that it is intermittent [5, 7–10]. Depending on whether the communication is continuous or intermitent, the resulting control is referred to as continuous control (CC) or intermittent control (IC), respectively. IC is advantageous since it reduces the computational and communication load [5, 11]. In the control literature, IC is known as sampled data systems (SDC) [11, 12], and is regaining interest for decentralized control. In the context of HMC, IC was suggested more than half a century ago to explain the psychological refractory period in tracking hand movements [13, 14]. However, the robustness of IC to model inaccuracies has not been addressed.
CC involves three main processes, as depicted in the left panel of Fig 1 (and further detailed in Continuous control): (i) Observer that combines the current noisy and delayed sensory measurement, y(t), with an internal model to estimate the time-delayed state x ^ ( t - τ ), where τ is the delay; (ii) Predictor that generates the predicted state, xp(t), given x ^ ( t - τ ); and (iii) Controller that derives the control signal u(t) from xp(t). IC requires two additional processes, as depicted in the right panel of Fig 1 (and further detailed in Intermittent control) [5, 12, 15]: (iv) Trigger that initiates information transfer between the observer and predictor at discrete times, tm (m = 1, 2, …), which are evenly spaced (clock-driven, periodic IC) or determined by specific events (event-driven, not considered here); and (v) System-Matched Hold (SMH) that generates the hold state xh(t), which approximates the predicted state in open-loop based on the internal model and the known control law. SMH facilitates IC since it accounts for the known control law between samples, rather than assuming that the control signal between samples is constant, as in the more common zero-order hold (ZOH).
Triggers are assumed to have a range of roles in HMC, including movement initiation (via the Go signal, for example) and transitions from one phase of locomotion to another [16, 17]. However, the main aspect of IC considered here is the use of SMH to approximate the predicted state in open-loop, in order to reduce the computational and communication burden. IC can explain important observations including: (i) highly variable timing of corrective sub-movements ; (ii) response to double stimuli with narrow pulses [5, 19], (iii) multi-peaked velocity profiles even in the presence of Gaussian noise  (though this can also be explained by CC with non-Gaussian process noise); and (iv) bursting activity in neural recording .
We focus on periodic IC, which involves clock-driven triggering. While HMC may involve event-driven triggering [14, 15, 22], it is generally assumed that refractory mechanisms set a lower limit on the sampling period [5, 13]. Hence, our analysis pertains to periodic IC at that lower limit.
To be a viable candidate for HMC, IC has to cope with inaccurately modeled plants, i.e., plants that differ from their internal model. Inaccurate models are common in daily activities, as when lifting lighter than expected loads or attempting to turn a locked steering wheel. Inaccurate models may also occur when cerebellar patients fail to update their internal models . Thus, our main objective was to evaluate and compare the effects of plant inaccuracies on IC and CC, and thus assess whether IC is a viable model for HMC. This was accomplished using both theoretical analysis and numerical simulations. Theoretical analysis was facilitated by considering time-invariant systems, which result from optimal control of time-invariant plants with respect to cost functions with time-invariant cost matrices. Numerical simulations were used to investigate movements with time-varying controller and observer gains, which are optimal when the cost function involve changing cost matrices, and in particular, when the cost matrices change at the desired reaching time. In either case, the controller and observer gains were determined using standard optimal control and estimation tools, as if the plant is accurately modeled.
Our second objective was to evaluate whether OFC, with either CC or IC, can provide an alternative model for the experimental results reported by Bhanpuri et al. , which involve reaching movements with inaccurately modeled plants. We demonstrate that both CC and IC reproduce well the main experimental results in , and in particular the observed overshoot and undershoot with heavier and lighter than expected mass, respectively.
Results include (i) theoretical criteria for determining the stability of CC or IC of inaccurately modeled linear time-invariant (LTI) plants with and without delays, under the assumption that the feedback and observer gains are also time-invariant, (ii) simulations of reaching movements that confirm those criteria, (iii) simulations of reaching movements with inaccurately modeled plants, subject to delays and time-varying feedback and observer gains, and (iv) overshooting and undershooting analysis of those simulations. Time-varying feedback and observer gains are optimal when the cost function involves time-dependent cost matrices, and in particular, when the cost matrices change at the desired reaching time.
A single-joint movement, such as flexing the elbow that was investigated in the dysmetria study in , is considered. For simplicity, the rotational movement is approximated by a translational movement, as in . The forearm and hand are modeled as a damped point-mass (with mass mPL and damping ratio γPL) to account for viscous damping at the elbow and the damping effect of any external device operated by the hand (e.g., exoskeleton , or joystick ). Stiffness is not included following , which concluded that inertia and damping accounted for much of the observed behavior. Actuation is generated by an over-damped second order muscle model suggested in [3, 25]. The internal model is the same as the actual plant model, except for the values of the mass and damping ratio, mIM and γIM, respectively, which may differ from the values mPL and γPL of the plant. The matrices defining the cost function, Eq (8), and the co-variances of the measurement and process noise were taken from  or motivated by . The model is detailed in Plant model.
Theoretical stability analysis
Theoretical stability analysis is conducted for LTI systems, i.e., LTI plants with time-invariant observer and feedback gains. Time-invariant feedback gains result from optmizing cost functions with time-invariant cost-matrices that penalize deviations from the target in the same way at all times (Eq (8)).
The standard equations for the resulting observer-predictor-controller (OPC) system  are briefly described below and detailed in Materials and Methods. Novel stability criteria are derived for periodic IC of inaccurately modeled plants. Lemma 1 focuses on the delay-free case, while Lemma 2 considers the effects of delays. The stability criteria are evaluated for translational reaching movements with the simple hand model mentioned above. The results are compared with corresponding stability criteria for CC, which are provided for completion in Stability of continuous LTI systems.
The dynamics of the LTI plant is described by the system matrix A ¯ and control matrix B ¯ (Eq (1)), which may differ from the system matrix A and control matrix B of the internal model (Eq (2)):
where x ∈ Rn is the state of the plant, u(t) ∈ Rm is the control signal, w(t) ∈ Rn is the process noise, and xIM ∈ Rn is the state of the internal model.
The compound effect of process and measurement delay is accounted for by introducing measurement delay τ:
where y(t) ∈ Rq is the measurement and v(t) ∈ Rq is the measurement noise. We also analyzed the effect of process delay and verified that the stability analysis is the same. Growing evidence suggests that HMC is subject to signal-dependent process and measurement noise [3, 4]. However, for simplicity, both are assumed to be white Gaussian noise with time-independent covariance matrix W and V, respectively (as in the examples in ).
As detailed in Continuous control, the observer generates the delayed estimated state, x ^ ( t - τ ), by combining the internal model and the delayed measurement y(t) according to the observer gain L(t) (Eq (9)). When τ > 0, a predictor is required to predict the current state xp(t) (Eq (10), left panel of Fig 1). The controller generates the control signal from either the estimated state (in delay-free systems, when τ = 0) or the predicted state (in time-delayed systems, when τ > 0), given the feedback gain K(t) (Eq (11)).
IC performs predictions only at discrete times tm (right panel of Fig 1). The predictor generates xp(tm), which provides the initial condition for the hold state, xh(t). Between samples, the hold state evolves continuously according to the SMH, Eq (13), captured by the feedback matrix:
Finally, the controller determines the control signal from the hold state (Eq 12). The evolution and reset of the hold state during IC of reaching movements with accurately and inaccurately modeled plants are illustrated in S1 File. The resulting control signal is also depicted and compared to the control signal generated by CC, to clarify the differences between the two control methods.
Periodic IC: Delay-free systems
The control of delay-free systems does not require a predictor, so the dynamics of the system depends only on the state of the plant and the state of the observer, captured by x o v ( t ) ≡ [ x ( t ) ′ x ^ ( t ) ′ ] ′. In periodic IC, the state of the observer is sampled at tm = mh, where h is the sampling period. Thus, stability depends on the discrete dynamics of xov[m]≡xov(tm), specified in Lemma 1.
Lemma 1. The discrete dynamics of delay-free LTI systems with periodic IC and sampling period h is given by xov[m + 1] = Apxov[m] + wov[m] where wov[m] is the discrete noise term,
and S is the solution of the Sylvester equation γS + Sβ + α = 0 with
The proof, detailed in S2 File, is based on solving the dynamics of xov(t), specified by Eq (14)), in continuous time, discretizing the solution at period h, and using the SMH (Eq (13)) to determine the hold state and resulting control signal between samples (Eq (12)). A similar approach was used in the networked system design , though here we rely on Sylvester equation to derive the discrete system dynamics.
From Lemma 1, the stability of delay-free IC systems with sampling period h can be determined from the eigenvalues of Ap specified in Eq (5). Since this is a discrete system, asymptotic stability is guaranteed when all the eigenvalues are inside the unit circle and is lost otherwise.
Lemma 1 was used to evaluate the stability of translational reaching movements with the muscle actuated hand model detailed in Plant model. This model includes 4 state variables (position, velocity, the force generated by the muscle and an internal state of the muscle), so Ap (Eq (5)) has 8 eigenvalues. The observer and feedback gain matrices L and K were optimized as if the plant is modeled accurately (see Materials and methods for more details). The effect of measurement noise was evaluated by considering both the nominal measurement noise (with the nominal covariance matrix V = V0 detailed in Plant model) and less noisy measurements with V = V0/10.
Considering movements of a forearm and hand of mass m = 1.5[Kg], we evaluated the eigenvalues of Aov for mIM = 2[Kg], to account also for the load. For completion, eigenvalues were evaluated for mPL ∈ (0.2, 4)[Kg], but for normal HMC only the range mPL ≥ 1.5[Kg] is relevant. Fig 2 depicts the maximum absolute eigenvalue as a function of mPL, for different sampling periods h (different blue lines, plotted with respect to the left y-axis), and two covariance matrices for the measurement noise (left panel: V = V0, right panel: V = V0/10). Non-smooth changes in the maximum absolute eigenvalue occur at the transition between different dominant eigenvalues.
For heavier than expected mass (mPL > mIM) the system remains stable independent of the sampling period h and measurement noise. For lighter than expected mass with V = V0 (left panel), stability is lost for mPL < 0.34[Kg] or mPL < 0.5[Kg] when h = 0.05[sec] or h = 0.1[sec], respectively, but is maintained throughout the investigated range when h = 0.2[sec]. With less noisy measurements (V = V0/10, right panel), stability is lost for mPL < 0.54[Kg] or mPL < 0.78[Kg] when h = 0.05[sec] or h = 0.1[sec], respectively, but is again maintained when h = 0.2[sec]. These results were verified in simulations of translational reaching movements to a target 0.2[m] away. Fig 3 depicts simulated trajectories with mIM = 2[Kg] and mPL = 0.4[Kg], when h = 0.1[sec] or h = 0.2[sec]. In agreement with the stability analysis in Fig 2, IC is stable when h = 0.2[sec] but is unstable for the shorter sampling interval h = 0.1[sec]. This holds with either the nominal measurement noise (V = V0, left panel) or less noisy measurements (V = V0/10, right panel).
Stability analysis of CC was performed for comparison based on the eigenvalues of Aov, Eq (15), which specifies the dynamics of xov(t) in continuous time, as detailed in Stability of continuous LTI systems. The stability of continuous systems is lost when the real part of any of the eigenvalues become positive. The maximum real part of the 8 eigenvalues of Aov is plotted on the corresponding panel of Fig 2 with respect to the right y-axis. Note that the two y-axes are aligned so the same y-level line marks the stability thresholds for both systems (unity on the left y-axis, and zero for the right y-axis). CC remains stable throughout the evaluated range when V = V0, but become unstable when V = V0/10 and mPL < 0.5[Kg], as verified by simulations in Fig 3.
We note that the range of stability becomes narrower when the measurements are less noisy, even though in that case the observer relies more on the measurements and less on the internal models, which may be inaccurate (i.e., the observe gains L are larger). The narrower stability range may be attributed to the over-reaction to the unexpected measurements.
Periodic IC: Time-delayed systems
The control of time-delayed systems involves a predictor that predicts the current state based on the estimation of the delayed state provided by the observer [27, 28]. In periodic IC of time-delayed systems, the predictor samples the observer at tm ≡ mh, receives x ^ [ m ] ≡ x ^ ( t m - τ ), and generates xp[m]≡xp(tm) (right panel of Fig 1). Considering the case where τ < h, the overall dynamics of the resulting OPC system depends on xtot[m] = [xov[m]′ xp[m − 1]′]′, as specified in Lemma 2.
Lemma 2. The discrete dynamics of LTI systems with delay τ under periodic IC with sampling period h > τ is given by xtot[m + 1] = Atotxtot[m] + wtot[m], where
Ao is defined by Eq (6), AF is defined by Eq (4), and Si, i = 1, 2, 3 are solutions to three Sylvester equations γiSi + Siβi + αi = 0 with
The proof, detailed in S2 File, is based on developing a difference equation for xp[m] that depends on x ^ [ m ] and a difference equation for xov[m] that depends on xp[m] and xp[m − 1]. The dynamics of xtot[m] is derived by combining those two difference equations.
From Lemma 2, the stability of time-delayed systems with periodic IC can be determined from the eigenvalues of Atot specified in Eq (7). Since this is a discrete system, asymptotic stability is guaranteed when all the eigenvalues are inside the unit circle, and is lost otherwise.
The left panel of Fig 4 depicts the maximum absolute eigenvalue of Atot when V = V0, h = 0.2[sec], and mIM = 2[Kg], as a function of mPL ∈ [0.2, 4][Kg], for different delays τ < h. A delay of τ = 0.1[sec] causes the system to become unstable when mPL < 0.5[Kg]. Interestingly, the system remains stable over all the investigated range with a longer delay of τ = 0.15[sec] (see also first row of Table 1). These results were verified in simulations, as depicted in the left panel of Fig 5 for reaching movements with mIM = 2[Kg] and mPL = 0.4[Kg]. In agreement with the stability analysis, IC is stable when τ = 0.15[sec] and unstable for shorter delay of τ = 0.1[sec].
Stability of time-delayed systems under CC was analyzed for comparison. The analysis was facilitated by converting the time-delayed continuous system to an equivalent time-delayed discrete-time system using the standard zero-order hold (ZOH), as detailed in Stability of continuous LTI systems. The equivalent discrete-time system is stable as long as all the eigenvalues of Aex, specified in Eq (16), are inside the unit circle. Note that an eigenvalue λd of a discrete-time system depends on the discretization time Δ, and evolves as λ d k at t = kΔ. For comparison with the eigenvalues of IC, consider the evolution of λd after t = h = 0.2[sec]. Here we used Δ = 10−3, so λd = 1.001, for example, would evolve to 1.001200 = 1.22.
The right panel of Fig 4 depicts the maximum absolute eigenvalue of Aex, indicating that time-delayed systems with τ = 0.001/0.1/0.15[sec] become unstable when mPL < 0.22/0.5/0.48[Kg], respectively (first row in Table 1). These results were verified in simulations, as depicted in the right panel of Fig 5 for reaching movements with mIM = 2[Kg] and mPL = 0.4[Kg]. In agreement with the stability analysis, CC is unstable with both τ = 0.1[sec] and τ = 0.15[sec].
The stability analysis in Fig 4 was repeated for different mIM and γIM, and the resulting stability thresholds are listed in Table 1. Neither CC nor IC is consistently superior to the other in terms of the range of loads over which the system remains stable. CC is superior when mIM = 6[Kg], while IC is superior when mIM = 4[Kg] and τ = 0.1[sec] or mIM = 2[Kg] and τ = 0.15[sec]. The effect of the delay is also inconsistent. For both mIM = 4[Kg] and mIM = 6[Kg], a system that is stable with τ = 0.15[sec] remains stable when the delay is shortened to τ = 0.1[sec], but this is not the case when mIM = 2[Kg], as indicated by Figs 4 and 5.
Reaching movements with time-varying gains
Human reaching is performed under time constraints that can be optimally satisfied with time-dependent feedback and observer gains [6, 29]. This hinders theoretical analysis, so we investigate the effects of plant inaccuracies using numerical simulations. The simulated task was to reach a target located 0.2[m] away from the initial position at T = 0.6[sec], with the hand model detailed in Plant model and measurement delay of τ = 0.15[sec]. Qualitative aspects of the responses are compared with previously published experimental results, both with ataxia patients and with healthy participants subjected to perturbations in the inertia and damping of the manipulated exoskeleton .
Reaching movements with inaccurately modeled plants
Figs 6 and 7 compare the trajectories generated by CC (solid lines) versus periodic IC (with h = 0.2[sec], dashed lines) when reaching with inaccurately modeled mass or damping, respectively. The parameters of the internal model were kept constant at the nominal values mIM = 2[Kg] and γIM = 8[Nsec/m]), while the parameters of the actual plant, mPL, and γPL, were smaller, the same, or larger than expected.
Both figures demonstrate that when the internal model is accurate, the responses generated by CC and by periodic IC are very similar, even though the sampling period, h = 0.2[sec], is third of the desired reaching time T = 0.6[sec]. This can be attributed to the accuracy of the SMH (Eq (13)) that IC uses to approximate the predicted state between samples. When the plant differs from the internal model, the responses generated by CC and IC deviate from each other, though the main characteristics remain the same. In particular, both controllers overshoot the target when the plant is heavier or the damping is smaller than expected. The effect of the mass may seem counter-intuitive, but agrees well with experimental results and other models . When the plant is lighter than expected, i.e., mPL = mIM − 0.5[Kg] (Fig 6), both CC and IC result in oscillatory trajectory, which are more pronounced with IC.
Overshooting and undershooting analysis
Overshooting or undershooting the target, known also as hypermetria and hypometria, respectively, are typical aspects of dysmetria. Following a study of motor dysmetria in ataxia patients during elbow movements, we quantified dysmetria by the difference between the target position and the position at the time of first correction . The latter was defined as the first time the velocity decreased below a threshold of 0.01[m/sec], or reached a local minimum after the peak velocity. The positions from which dysmetria was computed are marked by crosses (CC) or circles (IC) in Figs 6 and 7. It is evident that inaccurately modeled mass and damping have opposite effects on dysmetria according to both CC and periodic IC. In particular, while a lighter/heavier than expected plant results in undershoot/overshoot (Fig 6), a smaller/larger than expected damping ratio results in overshoot/undershoot (Fig 7).
Following , we evaluated the dependence of dysmetria on early velocity, which reflects the preplanned control before feedback driven corrections. Early velocity was defined as the velocity 0.15[sec] after movement onset, i.e., after the velocity exceeded 0.05[m/sec]. The analysis in Figs 6 and 7 was repeated for N = 9 different values of mPL ∈ (1.5, 2.5)[Kg] evenly spaced around MIM = 2[Kg] and separately for N = 9 different values of γPL ∈ (6, 10)[Nsec/m] evenly spaced around γIM = 8[Nsec/m]. Dysmetria and early velocity were extracted from each simulation and plotted against each other in Fig 8. Under both CC and periodic IC, dysmetria and early velocity are negatively correlated when reaching with inaccurately modeled mass (left panel, Fig 8), and positively correlated when reaching with inaccurately modeled damping ratio (right panel, Fig 8). Those correlations are in agreement with experimental results with healthy subjects , and the feedforward / feedback computational model suggested there.
Discrepancies between CC and IC in the accurate case (mPL = mIM and γPL = γIM) are due to the local velocity minima that emerges in the trajectories of IC at this case. Thus, the position of first correction, from which dysmetria was calculated, had been determined by different criteria: velocity decreasing below 0.01[m/sec] in CC and local velocity minima in IC. The relationship between dysmetria and early velocity is highly linear, as evident in the linear regression results listed on the graphs, except when mPL is much lower than mIM. Deviations from linearity in this range can be attributed to the early peak in the position that occurs when mPL is much lower than mIM, as is evident in Fig 6 for mPL = 1.5[Kg].
The negative correlation between dysmetria and early velocity is of particular interest since it was observed across a group of ataxia patients . Fig 9 demonstrates this negative correlation for two other sets of parameters. In both cases, the cost function was modified to reduce the penalty on deviations from the desired null velocity for t > T, with Q∞ = diag(1, 0.12, 0.022, 0) instead of the nominal Q∞ = diag(1, 0.22, 0.022, 0). The left and right panels depict the dysmetria analysis for mIM = 2[Kg] and mIM = 4[Kg], respectively. With the modified penalty, the correlation between dysmetria and early velocity under IC is more linear (left panel of Fig 8 versus the left panel of Fig 9). This can be attributed to the smaller initial peak in the position when mPL is much lower than mIM.
Materials and methods
Optimal control minimizes a cost function J that usually penalizes for deviations from the desired state and for the control effort. Stable linear systems converge to the origin, which is redefined to be the target state by change of variables (e.g., position is measured from the target position). Thus deviations from the origin are just the norm of the state vector, which can be weighted and computed in different ways. A standard cost function, which facilitates analysis, is the quadratic cost function, based on the weighted quadratic norm :
where T is the desired time to reach the target, also known as terminal time or horizon. When deviations from the target state are important at all times, the cost function in Eq (8) is extended to the infinite horizon case by taking the limit T → ∞.
The first term in the integral penalizes for deviations of the state from the origin weighted by a positive semidefinite matrix Q(t). The diagonal terms in this matrix are the weights of the square magnitude of the corresponding state variables, while other terms are usually set to zero. The second term in the integral penalizes for the square magnitude of the effort weighted by a positive definite matrix R(t). When the control signal is scalar, as in our case, R(t) is a positive scalar function of time. The term outside the integral penalizes for deviations from the origin at the (finite) terminal time, weighted by a positive semidefinite matrix S(T).
When T is finite, the cost function does not constrain the feedback gains for t > T. Nevertheless, those gains are important when the target is not reached at t = T, as is usually the case when controlling inaccurately modeled plants. Here we extend the gains at t = T to t > T, and assure that they do not vanish by proper selection of S(T). As detailed below, S(T) is selected so this is equivalent to optimizing the infinite horizon cost function, with Q(t > T) = Q∞ and Q(t < T) = Q(t) (see Plant model).
Optimizing Eq (8) involves two sub-processes: observer that estimates the state and controller that generates the control signal from the estimated state . As mentioned above, we simplify the analysis by assuming that the noise is signal independent white Gaussian noise, despite physiological evidence that the noise might be signal-dependent [3, 4]. Under this simplifying assumption, the two sub-processes can be optimized independent of each other. When there are delays in the system, as in HMC, a third intermediate sub-process, known as the predictor, is required to predict the current state from the delayed estimated state. The control signal in that case is derived from the predicted, rather than the estimated, state. The resulting three sub-processes are also known as the Observer-Predictor-Feedback (OPF) structure .
Observer combines the internal model (Eq (2)) and delayed measurement (Eq (3)) to generate the estimated state x ^ according to [6, 27, 29]:
where L(t) is the observer gain matrix. Under the assumption that the internal model is accurate, the optimal observer gain is the Kalman gain L(t) = P(t)C′ V−1(t) where P(t) is the solution to the Riccati differential equation P ˙ ( t ) = A P ( t ) + P ( t ) A ′ - P ( t ) C ′ V - 1 C P ( t ) + W [6, 29]. For infinite horizon, the Kalman gain is time-invariant L(V, W) = P∞C′ V−1 where P∞ is the solution of the algebraic Riccati equation 0 = AP∞ + P∞A′ − P∞C′ V−1CP∞ + W [6, 29].
Predictor predicts the current state, xp(t), given the estimated state, x ^ ( t - τ ), and the control signal u(σ) for σ ∈ [t − τ, t), based on the internal model (Eq (2)):
Controller generates the control signal u(t) given the predicted state:
where K(t) is the feedback gain matrix. For infinite horizon with time-invariant cost matrices (Q(t) = Q∞ and R(t) = R∞), the optimal feedback gain is time-invariant K = K ∞ ( Q ∞ , R ∞ ) ≡ R ∞ - 1 B ′ S ∞ where S∞ is the solution of the algebraic Riccati equation 0 = A ′ S ∞ + S ∞ A - S ∞ B R ∞ - 1 B ′ S ∞ + Q ∞ [6, 29]. For finite terminal time, the optimal feedback gain matrix is K(t) = R−1(t)B′ S(t) where S(t) is the solution to the Riccati differential equation: - S ˙ ( t ) = A ′ S ( t ) + S ( t ) A - S ( t ) B R ( t ) - 1 B ′ S ( t ) + Q ( t ) [6, 29] that is solved backward, starting with S(T). Due to plant inaccuracies, the plant may not reach the desired state at t = T, so feedback gains are relevant past the finite terminal time. For continuity, K(t > T) = K(T), and by setting S(T) = S∞(Q∞, R∞), K(T) = K∞(Q∞, R∞). This is equivalent to optimizing an infinite horizon cost funtion with Q(t > T) = Q∞ and Q(t < T) = Q(t).
IC performs predictions at discrete times tm, which can be evenly spaced (periodic IC, tm = mh, where h is the sampling period) or event driven (not considered in this work). At tm, the predictor receives x ^ ( t m - τ ) from the observer and generates xp(tm) according to Eq (10). The latter provides the initial condition for the hold state, xh(t) that determines the control signal:
Between samples, xh(t) evolves continuously according the feedback matrix (Eq (4)), defining the SMH :
Stability of continuous LTI systems
Stability is analyzed for LTI systems with time-invariant observer and controller gains, L and K. In that case, Eqs (1), (3) and (9) can be combined to describe the LTI dynamics of the overall state x o v ( t - τ ) = [ x ( t - τ ) ′ x ^ ( t - τ ) ′ ] ′:
where Ao and Bo are defined in Eq (6), and wov(t − τ) = [w(t − τ)′ Lv(t − τ)′]′ is the overall process noise.
In the delay-free case, the predicted state is the same as the estimated state. Thus, the control signal is determined from the estimated state given the state feedback gain K (Eq (11)), resulting in an autonomous system: x ˙ o v ( t ) = A o v x o v + w o v ( t ), where
Asymptotic stability is guaranteed as long as the real parts of all the eigenvalues of Aov are negative, and is lost otherwise.
For systems with delays, stability analysis is facilitated by converting the continuous system to an equivalent discrete-time system, using the standard zero-order hold . The discretization time Δ is selected so kτ = τ/Δ is an integer number. Thus, the overall dynamics of the equivalent discrete-time system depends on kτ samples of the predicted state. Defining the extended discrete state x e x ( k ) = [ x ( k - k τ ) ′ x ^ ( k - k τ ) ′ x p ( k - 1 ) ′ … x p ( k - k τ ) ′ ] ′, S3 File shows that with time-invariant observer and feedback gain matrices, Ld and Kd, respectively, xex(k + 1) = Aexxex(k) + wex(k), where
A ¯ d = e x p ( A ¯ Δ ), B ¯ d = A ¯ - 1 ( e x p ( A ¯ Δ ) - I ) B ¯, Ad = exp(AΔ) Bd = A−1(exp(AΔ) − I)B, and wex(k) is the discrete process noise.
Asymptotic stability is guaranteed as long as all the eigenvalues of Aex are within the unit circle, and is lost otherwise.
As mentioned above, the rotational movement of the hand around the elbow is approximated by a translational movement, as in . The translational position of the hand, ph, in response to the force, f, is governed by:
The force is generated by an over-damped second order muscle model with time constants μ1 = μ2 = μ = 0.04[sec]. The force generated by the muscle in response to the control signal u(t), corrupted by process noise, wg(t), is given by:
where g is an internal state of the muscle, and the variance of the process noise is σ g 2 = ( 0 . 46 [ N ] ) 2 .
Defining the hand state x ( t ) = [ p h ( t ) p ˙ h ( t ) f ( t ) g ( t ) ] ′, Eqs (17)–(19) can be expressed in the state-space representation of Eq (1) with:
The internal model is the same as the actual plant model, except for the values of the mass and damping ratio, mIM and γIM, respectively, which may differ from the values mPL and γPL of the plant.
The average mass of the human forearm and hand is m = 1.5[Kg] , so we consider mIM = 2–6[Kg] to account also for the mass of the load. The value of the linear damping ratio γhand was motivated by the angular damping ratio γθ ∈ (0.5–0.8)[Nmsec/rad] for elbow flexion/extension . Using standard conversion from joint to end-point damping , γ h a n d = γ θ / l a r m 2 where larm, is the length of the forearm. For larm ≊ 0.3[m], γhand ≊ (5.5 − 9) [Nsec/m]. Thus we use γIM = 8[Nsec/m] to account for viscous damping at the elbow and the effect of any external device operated by the hand.
Only the position and velocity of the hand are sensed, so C = [I2×2 02×2]. The nominal covariance matrix of the measurement noise is V 0 = d i a g ( [ σ p 2 σ v 2 ] ), where σ p 2 = ( 0 . 005 [ m ] ) 2 and σ v 2 = ( 0 . 04 [ m / s ] ) 2 are the variances of the (visual) measurement noise of the position and velocity, respectively. These values were estimated from , for mean position error of 0.09[m], and mean velocity of 0.28[m/sec].
The cost function, Eq (8), penalizes for control effort with R(t) = R∞ = 0.00001 and for deviations from the target with Q(t) = Q∞ = diag([1 0.22 0.022 0]) . While deviaions in the position are penalized most severely, this cost matrix also penalizes deviations from zero velocity and zero muscle force, which are required for staying at the target position. Finite reaching movements are imposed by setting Q(t < T) = 04×4, to allow for high velocity and force required to reach the target at T = 0.6[sec].
Intermittent control (IC) has been advocated as a more viable model for HMC [5, 19, 20] than continuous control (CC). Using system-matched hold (SMH), IC combines short-term open-loop control with intermittent closed-loop feedback to exploit the advantages of both. In particular, IC reduces the communication and computational burden associated with continuous prediction. For these reasons, IC is also gaining interest in the control literature, where it is known as sampled data control (SDC) . Furthermore, IC has been shown to provide better agreement with some HMC experiments as reviewed in the Introduction , and to increase robustness to perturbations during quiet and tiptoe standing [33, 34]. However, it is not clear how well IC maintains stability when controlling inaccurately modeled plants.
Our theoretical analysis, supported by numerical simulations, demonstrates that under some conditions, periodic IC may remain stable under a wider range of plant inaccuracies than CC, especially with less noisy measurements (right panels of Figs 2 and 3) or in the presence of delays (Figs 4 and 5). Interestingly, periodic IC with sample period of h = 0.2[sec] may remain stable under a wider range of inaccuracies with longer delays (e.g., with a delay of 0.15[sec] compared to a delay of 0.1[sec] or 0.001[sec]). While IC is not consistently superior to CC in terms of the range of loads over which the system remains stable, it is demonstrated to cope comparatively well with plant inaccuracies.
Simulations of reaching movements with time-varying gains, which are more relevant for HMC, demonstrate that periodic IC cope well with mismatches between the internal model and plant dynamics. Compared to CC, IC tends to generate more oscillatory responses and stronger overshoot/ undershoot. Interestingly, human experiments indicate that hypometric patients, whose internal model may not match the actual plant, also generate oscillatory reaching movements. Thus, both theoretical analysis and simulations suggest that IC is indeed a viable candidate for HMC.
Numerical simulations with different mass or damping inaccuracies agree well with previously reported experimental results with healthy participants where inertia and damping perturbations were introduced by an exoskeleton . In particular, both CC and IC result in the observed negative/ positive correlation between dysmetria and early velocity with inaccurately modeled mass/ damping, respectively. Damping inaccuracies result in a stronger effect on periodic IC than on CC, while mass inaccuracies result in a similar effect (left and right panels of Figs 8 and 9, respectively).
Experimental results with patients having cerebellar deficits were shown to be characterized by negatively correlated dysmetria and early velocity . Hence, it was hypothesized that cerebellar dysmetria is related to biased inertia in the internal model . A computational model based on feedforward/ feedback control was suggested to explain the experimental results. The feedforward control signal was derived from the desired trajectory based on the internal model of the plant, and then combined with a feedback control signal (Figure 4 in ). The desired trajectory was preplanned to generate a bell-shaped velocity profile (i.e., a minimum jerk profile ) with the desired movement amplitude and duration. In contrast, optimal control generates the trajectories online based on optimal state estimation and optimal control gains [3, 4, 6]. Optimal control gains are derived from a cost function that may account for both accuracy and control effort. The estimated state is derived from the internal model and sensory feedback, according to optimal Kalman gains that account for the processs and measurement noise. Here we demonstrate that OFC with either CC or IC provide an alternative model for the correlations between dysmetria and early velocity reported in .
Our analysis focused on periodic IC, though HMC may be dominated by aperiodic IC triggered by crossing prediction error thresholds [14, 15, 22]. The analysis of periodic IC is justified since refractory mechanisms are assumed to set a lower limit on the sampling period [5, 13]. Thus, our analysis considers the limit of small thresholds, when the lower limit on the sampling period is reached consistently. Comparing event- and clock-driven IC with the same number of updates per movement, event driven IC is expected to be more robust to plant inaccuracies since it would provide more updates when the deviations are larger, though this would require further investigation.
The assumptions that the measurement and process noise are Gaussian white noise (rather than signal dependent noise) would affect the observer and feedback gains. Those gains are also affected by other assumptions that are commonly made in HMC including: (i) weights in the cost function are time-independent, or at least vary only at the desired time to reach the target, (ii) damping is time-independent, contrary to evidence suggesting that it changes during the movement , and (iii) feedback gains are determined to optimize a quadratic cost function rather than other functions, especially those that lead to higher robustness. The effects of the observer and feedback gains were considered by analyzing cases with different measurement noise, which affect the observer gains, and simulating cases with different cost-functions, which affect the feedback gain.
Given the highly uncertain environment in which human operate, robust control may provide a more viable framework for HMC than optimal control. Optimal control is usually not robust to plant variations, since it is tailored to specific parameters. Here, we focused on optimal control due to its prevalence in theories of HMC. Nevertheless, robust controllers, based on H∞, as suggested for postural control , or robust satisficing, as suggested for two-alternative forced choice tasks , may be more relevant for HMC.
In summary, our main contribution is in demonstrating that IC is a viable model of HMC even when considering the effects of model inaccuracies. In addition, our study demonstrates that OFC with either CC or periodic IC reproduces the reported correlation between dysmetria and early velocity with inaccurately modeled mass or damping.
1. Todorov E, Jordan MI, Optimal feedback control as a theory of motor coordination, Nature neuroscience, 5(11), 1226–1235 (2002). doi: 10.1038/nn963 12404008
2. Todorov E, Optimality principles in sensorimotor control, Nature neuroscience, 7(9), 907–915 (2004). doi: 10.1038/nn1309 15332089
3. Todorov E, Stochastic optimal control and estimation methods adapted to the noise characteristics of the sensorimotor system, Neural computation, 17(5), 1084–1108, (2005). doi: 10.1162/0899766053491887 15829101
4. Shadmehr R, Krakauer JW, A computational neuroanatomy for motor control, Experimental Brain Research, 185(3), 359–381 (2008). doi: 10.1007/s00221-008-1280-5 18251019
5. Gawthrop P, Loram I, Lakie M, Gollee H, Intermittent control: a computational theory of human control, Biological cybernetics, 104(1), 31–51 (2011). doi: 10.1007/s00422-010-0416-4 21327829
6. Stengel RF, Optimal Control and Estimation. Courier Corporation (1994).
7. Ronco E, Arsan T, & Gawthrop P J., Open-loop intermittent feedback control: practical continuous-time GPC. IEE Proceedings-Control Theory and Applications, 146(5), 426–434, (1999). doi: 10.1049/ip-cta:19990504
8. Meyer DE, Abrams RA, Kornblum S, Wright CE, & Keith Smith JE, Optimality in human motor performance: ideal control of rapid aimed movements. Psychological review, 95(3), 340 (1988). doi: 10.1037/0033-295x.95.3.340 3406245
9. Neilson P. D., Neilson M. D., & O’dwyer N. J., Internal models and intermittency: A theoretical account of human tracking behavior. Biological Cybernetics, 58(2), 101–112, (1988). doi: 10.1007/bf00364156 3349110
10. Karniel A, Open questions in computational motor control, Journal of integrative neuroscience, 10(3), 385–411 (2011). doi: 10.1142/S0219635211002749 21960308
11. Garcia E, Antsaklis PJ, Montestruque LA, Model-based control of networked systems. Birkhuser (2014).
12. Montestruque LA, Antsaklis PJ, On the model-based control of networked systems, Automatica, 39(10), 1837–1843, (2003). doi: 10.1016/S0005-1098(03)00186-9
13. Vince MA, The Intermittency of control movements and the psychological refractory period, British Journal of Psychology. General Section, 38(3), 149–157, (1948). doi: 10.1111/j.2044-8295.1948.tb01150.x
14. Navas F, Stark L, Sampling or intermittency in hand control system dynamics, Biophysical Journal, 8(2), 252–302 (1968). doi: 10.1016/S0006-3495(68)86488-4 5639937
15. Gawthrop P J, Wang L. Event-driven intermittent control. International Journal of Control, 82(12), 2235–2248 (2009). doi: 10.1080/00207170902978115
16. McVea DA, & Pearson KG, Stepping of the forelegs over obstacles establishes long-lasting memories in cats. Current Biology 17, no. 16: R621–R623 (2007). doi: 10.1016/j.cub.2007.06.026 17714644
17. Hiebert GW, Whelan PJ, Prochazka A & Pearson KG, Contribution of hindlimb flexor muscle afferents to the timing of phase transitions in the cat step cycle, J. Neurophysiol. 75, 1126–1137 (1996). doi: 10.1152/jn.19126.96.36.1996 8867123
18. Doeringer JA, Hogan N, Intermittency in preplanned elbow movements persists in the absence of visual feedback, Journal of Neurophysiology, 80(4), 1787–1799 (1998). doi: 10.1152/jn.19188.8.131.527 9772239
19. van de Kamp C, Gawthrop PJ, Gollee H, & Loram ID, Refractoriness in sustained visuo-manual control: is the refractory duration intrinsic or does it depend on external system properties?. PLoS computational biology, 9(1), e1002843 (2013). doi: 10.1371/journal.pcbi.1002843 23300430
20. Fishbach A, Roy SA, Bastianen C, Miller LE, & Houk JC, Deciding when and how to correct a movement: discrete submovements as a decision making process. Experimental brain research, 177(1), 45–63 (2007). doi: 10.1007/s00221-006-0652-y 16944111
21. Groß J, Timmermann L, Kujala J, Dirks M, Schmitz F, Salmelin R, Schnitzler A, The neural basis of intermittent motor control in humans, Proceedings of the National Academy of Sciences, 99(4), 2299–2302 (2002). doi: 10.1073/pnas.032682099
22. Gollee H, Gawthrop P J, Lakie M, Loram ID. Visuo-manual tracking: does intermittent control with aperiodic sampling explain linear power and non-linear remnant without sensorimotor noise? The Journal of physiology, 595(21), 6751–6770 (2017). doi: 10.1113/JP274288 28833126
23. Bhanpuri NH, Okamura AM, Bastian AJ, Predicting and correcting ataxia using a model of cerebellar function. Brain, 137(7), 1931–1944 (2014). doi: 10.1093/brain/awu115 24812203
24. Benyamini M, Zacksenhouse M, Optimal feedback control successfully explains changes in neural modulations during experiments with brain-machine interfaces, Frontiers in systems neuroscience, 9 (2015). doi: 10.3389/fnsys.2015.00071 26042002
25. Winter D A, Biomechanics and motor control of human movement. John Wiley & Sons, (2009). doi: 10.1002/9780470549148
26. Saunders J. A., & Knill D. C., Visual feedback control of hand movements, Journal of Neuroscience, 24(13): 3223–3234, (2004). doi: 10.1523/JNEUROSCI.4319-03.2004 15056701
27. Kleinman D, Optimal control of linear systems with time-delay and observation noise. IEEE Transactions on Automatic Control, 14(5), 524–527 (1969). doi: 10.1109/TAC.1969.1099242
28. Kleinman DL, Baron S, Levison WH, An optimal control model of human response part I: Theory and validation, Automatica, 6(3), 357–369 (1970). doi: 10.1016/0005-1098(70)90051-8
29. Kwakernaak H, Sivan R, Linear optimal control systems, Vol. 1. New York: Wiley-interscience (1972).
30. NASA Handbook, NASA Human Integration Design Handbook (HIDH)-NASA (Vol. 3407). SP-2010 (2010).
32. Shadmehr R, & Wise SP, The computational neurobiology of reaching and pointing: a foundation for motor learning. MIT press (2005).
33. Asai Y, Tasaka Y, Nomura K, Nomura T, Casadio M, Morasso P, A model of postural control in quiet standing: robust compensation of delay-induced instability using intermittent activation of feedback control. PLoS One, 4(7), e6169 (2009). doi: 10.1371/journal.pone.0006169 19584944
34. Tanabe H, Fujii K, Suzuki Y, & Kouzaki M, Effect of intermittent feedback control on robustness of human-like postural control system. Scientific reports, 6, 22446 (2016). doi: 10.1038/srep22446 26931281
35. Flash T, Hogan N, The coordination of arm movements: an experimentally confirmed mathematical model, Journal of neuroscience, 5(7), 1688–1703 (1985). doi: 10.1523/JNEUROSCI.05-07-01688.1985 4020415
36. Yang Y, Pei L, Li H, An H∞, control model of human postural control in quiet upright standing, Control and Decision Conference (CCDC) IEEE.Chicago, 2483-2486 (2011).
37. Zacksenhouse M, Bogacz R, and Holmes P, Robust versus optimal strategies for two-alternative forced choice tasks, J Math. Psychology, 54:230–246 (2010). doi: 10.1016/j.jmp.2009.12.004