Comparison of linear and nonlinear feedback control of heart rate for treadmill running

ABSTRACT Heart rate can be used to define exercise intensity; feedback control systems for treadmills which automatically adjust speed to track arbitrary heart rate target profiles are therefore of interest. The aim of this study was to compare linear (L) and nonlinear (NL) controllers using quantitative performance measures. Sixteen healthy male subjects participated in the experimental L vs. NL comparison. The linear controller was calculated using a direct analytical design that employed an existing approximate plant model. The nonlinear controller had the same linear component, but it was augmented using static plant-nonlinearity compensation. At moderate-to-vigorous intensities, no significant differences were found between the linear and nonlinear controllers in mean RMS tracking error (2.34 vs. 2.25 bpm [L vs. NL], p=0.26) and average control signal power (51.7 vs.  m2/s2, p=0.16), but dispersion of the latter was substantially higher for NL (range 45.2 to 56.8 vs. 30.7 to  m2/s2, L vs. NL). At low speed, RMS tracking errors were similar, but average control signal power was substantially and significantly higher for NL (28.1 vs.  m2/s2 [L vs. NL], p<0.001). The performance outcomes for linear and nonlinear control were not significantly different for moderate-to-vigorous intensities, but NL control was overly sensitive at low running speed. Accurate, stable and robust overall performance was achieved for all 16 subjects with the linear controller. This points to disturbance rejection of very-low-frequency heart rate variability as the overriding challenge for design of heart rate controllers.


Introduction
Fitness training programmes often use heart rate to describe exercise intensity (Garber et al., 2011;Pescatello et al., 2014). A range of training protocols can be prescribed by combining intervals of varying duration and intensity (Weston, Taylor, Batterham, & Hopkins, 2014). During outdoor running, heart rate can be controlled manually by the runner observing a standard heart rate monitor or by feedback control using wearable sensor and mobile computing technology . For treadmill exercise, automatic feedback control of heart rate is particularly advantageous because the runner is then freed of the burden of having to continually adjust the treadmill speed to keep to the heart rate target. It is therefore important to develop heart rate control systems for treadmill exercise which are able to accurately track arbitrary heart rate target profiles by automatic and continuous adjustment of the plant control input (i.e. the treadmill speed).
It has recently been proposed that the principal design issue for feedback control of heart rate is disturbance CONTACT Kenneth J. Hunt kenneth.hunt@bfh.ch rejection of very-low-frequency heart rate variability, VLF-HRV (Hunt & Fankhauser, 2016;Sassi et al., 2015): accurate tracking of the heart rate target is important, but, at the same time, the control signal must not be excited too strongly in frequency bands where changes in the treadmill speed would be perceptible to and unacceptable for the runner. To address this challenge, a new approach was developed based on shaping of the plant input sensitivity function; in an experimental evaluation with 30 subjects, it was shown that this linear time-invariant approach was able to directly address the principal design challenge of VLF-HRV disturbance rejection and delivered robust and accurate tracking with a smooth, low-power control signal in all subjects (Hunt & Fankhauser, 2016). Moreover, controller calculation was based upon a single, approximate linear plant model that was not specific to any of the subjects tested, but was obtained previously in a separate system identification study (Hunt, Fankhauser, & Saengsuwan, 2015).
Other approaches have focussed not on disturbance rejection of VLF-HRV, but on the issues of parametric and structural plant uncertainty. By considering in detail the underlying physiological mechanisms which come into play during exercise, a nonlinear second-order statespace model of heart rate dynamics was proposed in which the control variable (the treadmill speed v) enters the plant quadratically, that is, in the form v 2 (Cheng, Savkin, Celler, Su, & Wang, 2008); controller synthesis was then based on piecewise linear approximation and a combination of linear-quadratic and H ∞ optimization; the method was experimentally tested with six healthy subjects during low-intensity walking.
An alternative control approach based on the v 2nonlinear model from Cheng et al. (2008) was developed and applied to both treadmill and cycle-ergometer exercise: see Scalzi, Tomei, and Verrelli (2012) and Paradiso, Pietrosanti, Scalzi, Tomei, and Verrelli (2013), respectively. This method countered the plant nonlinearity using the inverse nonlinearity, that is, the square-root function √ . , at the controller output. Although couched in the terminology of iterative learning control, the approach simply comprises, in its discrete-time implementation, a linear proportional-integral (PI) controller combined with static nonlinearity compensation. For the treadmill, this approach was tested with two healthy subjects, the PI gains having been selected empirically (Scalzi et al., 2012).
Both of the nonlinear control approaches reviewed above (Cheng et al., 2008;Scalzi et al., 2012) have interesting theoretical properties but, to date, there has been no experimental study which has performed a direct, headto-head comparison of linear and nonlinear strategies for heart rate control within a single subject cohort. Further, these previous studies have lacked objective, quantitative outcomes for analysis of controller performance upon which such a comparison could be based. The motivation for the present paper was the perceived need to directly address these limitations and, specifically, to determine whether nonlinear control brings any discernible benefits when compared to linear control.
To address the open questions identified above, the aim of the present work was to systematically compare linear (L) and nonlinear (NL) heart rate controllers using quantitative measures of root-mean-square tracking error and average control signal power. The NL controller used the static nonlinear compensation strategy of Scalzi et al. (2012), but with the improvement that the embedded PI controller parameters were obtained using a direct analytical design. The L controller was purposely constrained and parameterized to have a structure identical to the PI component of the NL controller; this strategy provided the means by which both the L and the PI/NL controller parameters could be calculated directly using an existing simple and approximate plant model (no plant identification was required for any of the subjects tested); this strategy also ensured that only a single factor was being compared in this study, namely the presence (NL) or absence (L) of the static nonlinear compensation function at the controller output.

Subjects and test protocol
Sixteen healthy male subjects participated in the experimental evaluation and comparison of the linear (L) and nonlinear (NL) heart rate control strategies (subject details - Table 1). Each subject was tested using the L and NL controllers according to a systematic protocol which was reviewed and approved by the Ethics Committee of the Swiss Canton of Bern (Ref. ; all subject gave written, informed consent prior to participation. The L and NL tests were carried out on separate days and the study design was counterbalanced by randomizing the order of presentation of each test condition, that is, L then NL vs. NL then L, for each subject. Each measurement was preceded by a 10-minute warm up at moderate intensity followed by a 10-minute rest. The formal measurement for each controller lasted 45 min and proceeded according to a standardized target heart rate profile (see Figure 4). The profile was selected to allow investigation of both steady-state regulation and dynamic tracking performance: for the first 35 min, target heart rate was varied by ±10 beats per minute (bpm) around a mid-level. The final 10 min were designed as a low-intensity cool-down period where target heart rate was set to 15 bpm below the lowest level used during the 0-35 min stage.
The mid-level target heart rate during the first 35 min, HR mid , was set individually for each subject at the border between the moderate and vigorous exercise intensity regimes: this boundary is defined as 76.5% of the age-related maximal heart rate prediction, thus HR mid = 0.765 · HR max where HR max 220 − age [bpm] (Pescatello et al., 2014;Shargal et al., 2015).
Controllers were implemented in a PC using the Matlab/Simulink Real Time Workshop (The Mathworks, Inc., USA) connected to a treadmill (model Venus, h/p/cosmos Sports and Medical GmbH, Germany) using a serial communication protocol. Heart rate was measured using a chest belt (model T34, Polar Electro Oy, Finland) and receivers integrated in the treadmill.

Primary outcome measures and statistical analysis
Comparison of the performance of the L and NL controllers was based on objective, quantitative outcome measures: root-mean-square tracking error (RMSE) for the heart rate HR; and average control signal power, defined formally as the average power of changes in the control signal v (P v ), where v is the treadmill speed reference. Thus, where i are the discrete indices of the evaluation data. HR nom is the target nominal heart rate response which was obtained by simulating the nominal linear closedloop transfer function. Both outcomes were calculated for two phases of each test: a main evaluation interval 420 ≤ t ≤ 1800 s (7-30 min), and a low-speed evaluation interval 2100 ≤ t ≤ 2700 s (35-45 min) (cf. Figures 4 and 7).
Prior to hypothesis testing, normality of the data sets was checked using the Kolmogorov-Smirnov test with Lilliefors correction. Paired two-sided t-tests were applied for normal data, and Wilcoxon signed rank tests otherwise, to determine whether there were any significant differences in the means of the two primary outcomes, RMSE and P v , between the two test conditions L and NL; the null hypothesis was that no differences existed. The significance level was set to α = 0.05; statistical analysis was performed using the Matlab Statistics and Machine Learning Toolbox (The Mathworks, Inc.).

Linear and nonlinear control structures
The linear control structure, L (Figure 1(a)), comprises a conventional single-degree-of-freedom feedback loop with a linear, time-invariant, discrete-time compensator C(z −1 ) and a nominal plant P.
The nonlinear control structure, NL (Figure 1(b)), has the same linear function C, but it is augmented with a static nonlinear (square-root) compensation function √ . in order to directly compensate the plant nonlinearity, that is, the putative influence of the control signal v (treadmill speed reference) within the plant in the quadratic form v 2 . To ensure that the gain of the overall nonlinear compensator at a chosen operating point v o is identical to the gain of the linear function C, a static gain f must be introduced whose value is the inverse of the small-signal gain of the square-root function at the operating point v o . Since v = √ x, where x is an intermediate signal introduced for convenience, the small-signal gain of the nonlinear compensation function is dv It can readily be shown that the nonlinear structure NL (Figure 1(b)) is identical to a linear PI compensator augmented with the static square-root compensation function (see Figure 1(c)), as proposed in the original references for this approach (Paradiso et al., 2013;Scalzi et al., 2012). To see this, consider the transfer function C PI from e to x in Figure 1 (4) The rhs of Equation (4) has the same form as the constrained transfer function C(z −1 ) (see Equation (6)), that is, the numerator is a polynomial of degree 1 in z −1 and the denominator is 1 − z −1 . To make the two structures equivalent, it is necessary to set C PI = fC (cf. the two structures in Figure 1(b) and 1(c)), viz.
These transfer functions are identical when k y + μ = fs 0 and −k y = fs 1 . Should it be so desired, these relations allow the PI gains k y and μ to be calculated, given the parameters s 0 and s 1 : k y = −fs 1 , μ = f (s 0 + s 1 ). However, in the present work, the L and NL control structures were implemented directly according to Figure 1(a) and 1(b).

Feedback design
The parameters of the linear compensator C were obtained by following an algebraic approach to obtain a closed-form analytical solution. The algebraic problem was set up and constrained so as to obtain a transfer function for C with the same structure as the PI compensator forming part of the nonlinear method in Figure 1(c). Thus, as noted above, the sought-after compensator is required to have a numerator polynomial S(z −1 ) of degree 1 in z −1 and a denominator R(z −1 ) equal to 1 − z −1 , whereby the is a linear transfer function and P is the plant. The principal signals in each structure are the heart rate reference HR * , the actual heart rate HR, the treadmill speed reference v (plant control input), a disturbance d which mainly represents heart rate variability, and the error e. (a) Linear control structure, L. (b) Nonlinear control structure, NL. The linear transfer function C is the same as in L, part (a). f is a static gain whose value is the inverse of the small-signal gain of the square-root function at the operating x is an intermediate signal. (c) Equivalent nonlinear control structure: linear PI controller C PI with static nonlinearity compensation (Scalzi et al., 2012). k y is a proportional gain and μ is the integrator gain. For a proof of the equivalence of the structures in parts (b) and (c), see Section 2.3, Equations (4)-(5).
latter term gives integral action, that is where s 0 and s 1 are real coefficients to be determined. The control design approach is purposely based upon a linear first-order plant model with steady-state gain k and time constant τ , expressed in continuous (P c (s)) or discrete where the double arrow denotes discretization with sample period T s . The discrete model parameters, expressed in terms of k and τ , are Here, the nominal plant model was taken from a previous identification study as the average of 48 individual models (Hunt et al., 2015). The methodology employed for parameter estimation and model validation is fully detailed in that reference. Briefly, the 48 individual models were obtained empirically from 24 subjects who were each tested at two different mean speed levels representing moderate and vigorous exercise intensities. Individual models of heart rate response were first-order, linear, time-invariant transfer functions of the form given by Equation (7). Model parameters k and τ were obtained using least-squares optimization; model validation was carried out using a normalized root-mean-square error (NRMSE) between the model and measured outputs ('model fit') and also by the absolute RMSE. The average model thus obtained was that is, k = 24.2 bpm/(m/s), τ = 57.6 s (cf. P c in Equation (7)). Thus, in the present work, no system identification was carried out and the nominal plant model was not specific to any of the 16 subjects tested. It can readily be proven that, for the first-order plant model P d in Equation (7), the controller transfer func- (6) is the correct structure to allow full and arbitrary placement of two closed-loop poles. A complete proof of this result can be found elsewhere ; in brief, the foregoing problem formulation results in the characteristic equation The closed-loop poles are the roots of the characteristic polynomial , which, from Equation (10), has degree 2: It can be seen that there must exist a unique solution s 0 , s 1 to the algebraic problem defined by Equations (10) and (11) because the degree of both sides is the same and this degree, 2, is equal to the number of unknowns. The unique solution of this pole assignment problem is found by equating coefficients of like degree on both sides of Equation (10) to give Here, the characteristic polynomial was determined from the discrete-time equivalent of a second-order transfer function with a desired 10-90% rise time t r and critical damping ζ = 1, thus giving A complete proof of these results can again be found in .

Controller calculation
Calculation of the controller parameters s 0 and s 1 in C(z −1 ), Equation (6), is straightforward: (1) The sample period and desired rise time were selected as T s = 5 s and t r = 236 s (see below for the rationale behind these choices).
(2) The nominal plant model was calculated according to Equations (7)-(9) to be that is, b 0 = 2.0121 and a 1 = −0.9169. (3) φ 1 and φ 2 were calculated from Equation (13)  This procedure gave the controller transfer function (6), which was implemented as part of the L and NL control strategies (Figure 1(a) and 1(b), respectively). For the NL controller, the speed operating point was selected as v o = 2.5 m/s and the required scaling factor f was calculated from Equation (3) to be f = 2v o = 5. The choice v o = 2.5 m/s was made because the observed mean speed for 30 subjects in a previous study of similar experimental design wasv = 2.50 m/s (Hunt & Fankhauser, 2016).
(As noted following Equation (5), the computed parameter values for s 0 and s 1 in Equation (15) and for f would, incidentally, lead to the following PI controller parameters in the equivalent structure of Figure 1(c): k y = −fs 1 = 0.1222, μ = f (s 0 + s 1 ) = 0.0117. This procedure of designing the linear part of the controller using a unique, closed-form analytical solution to determine the PI parameters is advantageous, because an empirical, trial-and-error tuning procedure for k y and μ is avoided.) The rationale for the choice of sample period as T s = 5 s was that it is recommended for digital control to have ∼ 4 to 10 samples over the rise time of the plant (Åström & Wittenmark, 2011). The nominal plant used here had τ = 57.6 s, Equation (9), but it was previously observed that this value could be as low as τ =27.3 s during moderateto-vigorous intensity treadmill exercise (Hunt et al., 2015).

Frequency-domain analysis
The input sensitivity function U o links the HRV disturbance signal d to the control signal, that is, to the treadmill speed reference signal v ( Figure 1); it is defined by Substituting for C and P d from Equations (6) and (7), and employing Equations (10) and (11), This shows that U o is causal (but not strictly causal) or, equivalently, that U o is proper in z (but not strictly proper). As a consequence, |U o | does not roll off to zero at high frequency, but tends instead to a finite value (cf. Figure 2, thick blue line). By way of contrast, Figure 2 also shows the input sensitivity function gain |U sp | (thick red line) corresponding to a strictly proper compensator employed elsewhere (Hunt & Fankhauser, 2016), which does roll off to zero at high frequency; that compensator was demonstrated in Hunt and Fankhauser (2016) to give accurate, stable and robust overall control performance. This led to the rationale for the choice of t r = 236 s above: t r was optimized to make the bandwidth of U o the same as the bandwidth of U sp , namely frequency 0.0046 Hz, whence these two functions intersect at this frequency and at the −3 dB gain (Figure 2).  Hunt and Fankhauser (2016); U o was designed to have the same bandwidth as U sp , at frequency 0.0046 Hz, which is why these two functions intersect at 0.0046 Hz and −3 dB.
To obtain the input sensitivity function U NL for the nonlinear controller structure (NL, Figure 1(b)), the function C has to be scaled by the product of the factor f = 2v o and the small signal gain of the square-root function, 1/2v, that is, by v o /v, which gives U NL is seen to be identical to U o at the operating point v = v o , but diverges from U o at off-nominal speeds: see Figure 2, which shows |U NL | over the speed range 1 to 4 m/s, together with |U o |; clearly, the gain |U NL | becomes much larger at lower speeds due to the effect of the smallsignal gain of the square-root function, 1/2v, on U NL as described by Equation (18).

Main evaluation interval
There were no significant differences in the primary outcomes RMSE and P v for the linear (L) and nonlinear (NL) controllers over the main evaluation interval 420 ≤ t ≤ 1800 s: the mean RMS tracking errors were, respectively, 2.34 and 2.25 bpm (p = 0.26, Table 2); the means of the average control signal powers were 51.7 and 60.8 × 10 −4 m 2 /s 2 (p = 0.16, Table 2). Dispersion of RMSE values for the L and NL cases was similar (Table 2, where similar standard deviations can be observed for RMSE, and Figure 3(a)), but the dispersion of P v values was substantially and strikingly higher for the nonlinear controller: the ranges for P v were 45.2 to 56.8 vs. 30.7 to 108.7 × 10 −4 m 2 /s 2 , L vs. NL (Figure 3(b)); the standard deviations of P v were 3.3 vs. 24.1 × 10 −4 m 2 /s 2 (Table 2 and Figure 3(b)).
The non-significant differences in the primary outcomes are signified by p > 0.05 in both cases, and also by the inclusion of the value 0 in the 95% confidence intervals for the mean differences (Table 2, Figure 3(a) and 3(b)).
The overall mean running speed during the main evaluation interval, averaged across all 16 subjects, was 2.40 m/s for L and 2.42 m/s for NL. Since, for the nonlinear controller structure NL, the overall compensator parameters are dependent upon running speed v, comparative test results are presented for both the L and NL controllers for the subjects who had the lowest, closestto-nominal (v o = 2.5 m/s), and highest individual mean speedv during the evaluation period ( Figure 4). These data highlight the strong dependence in the NL controller of average control signal power, P v , on speed: for the NL controller, and for the three subjects highlighted, P v was on a very wide range, being 30.7, 56.3 and 108.7 ×    Table 2); the green lines link the sample pairs from each subject; the red horizontal bars depict mean values. D is the difference between the paired samples: D = L − NL. MD is the mean difference (red horizontal bar), with its 95% confidence interval (CI) in blue. Inclusion of the value 0 within the 95% CIs signifies non-significant differences between the means; this conforms with p > 0.05 in both cases (Table 2).
10 −4 m 2 /s 2 for the highest, middle and lowest mean speeds (Figure 4(f), 4(d) and 4(b), respectively). For the L controller, on the other hand, P v was very similar for all three subjects: 54.9, 50.9 and 56.8 × 10 −4 m 2 /s 2 , highest to lowest speeds, Figure 4(e), 4(c) and 4(a), respectively. For the nonlinear controller NL, the dependence of P v on speed v was further investigated using correlation analysis. There was a very strong negative linear correlation between v and P v , with correlation coefficient r = −0.92 (p = 6.2 × 10 −7 ) (cf. Figure 5). Since the smallsignal gain of the nonlinear compensation element in the NL structure is equal to 1/(2v) (Section 2.3), and since P v involves the square of changes in the control signal, Equation (2), it would be expected that P v ∝ 1/v 2 and that a nonlinear fit of this form to the observed P v − v relationship would give a higher fit: such a nonlinear fit gave R = 0.96 ( Figure 5; here, R is the coefficient of determination). For the NL controller, there was very weak or no correlation between v and RMSE; for the L controller there were very weak or no correlations between v and P v and between v and RMSE; (r < 0.31 for these three cases).

Low-speed evaluation interval
The apparent tendency for the control signal (treadmill speed reference v) to become very sensitive with the NL controller when speed is low can also be seen during the final 10-minute low-speed evaluation interval (2100 ≤ t ≤ 2700 s) of each test described above: the level of activity in the control signal increases substantially as mean speed reduces (Figure 4(f) to 4(d) to 4(b)). The behaviour of the control signal with the L controller, in contrast, was similar during the final 10-minute periods for the three cases shown (Figure 4(e), 4(c) and 4(a)) and for all other subjects.
There was no significant difference in RMSE for the linear (L) and nonlinear (NL) controllers over the formal low-speed evaluation interval: the mean RMS tracking errors were, respectively, 2.29 and 2.30 bpm (p = 0.94, Table 3, Figure 6(a)). However, there was a very substantial and significant difference in P v between the L and NL cases during the low-speed phase: the means of the average control signal powers were 28.1 and 138.7 × 10 −4 m 2 /s 2 , respectively (p = 0.00024, Table 3,  The nonlinear fit, which is of the form ∝ 1/v 2 , is substantially higher, with R = 0.96. (r is the correlation coefficient; R 2 is the coefficient of determination.) Figure 6(b)). This low-speed outcome evaluation considered n = 13 subjects (Table 3); three subjects were excluded from this part of the analysis because, during the low-speed phase, their heart rate remained above the target heart rate even when the treadmill speed reference v was set by the controller to its lower limit, thus rendering the system open-loop. Dispersion of P v values was again substantially and strikingly higher for the nonlinear controller: the ranges for P v during the low-speed evaluation phase were 18.2 to 45.9 vs. 21.4 to 469.1 × 10 −4 m 2 /s 2 , L vs. NL (Figure 6(b)); the standard deviations of P v were 8.7 vs. 142.3 × 10 −4 m 2 /s 2 (Table 3 and Figure 6(b)).
The test results for the subject with the highest value of P v = 469.1 × 10 −4 m 2 /s 2 clearly illustrate that the NL controller can be hypersensitive at low speed (Figure 7): the overall performance of the linear controller was good throughout the whole test, including the lowspeed phase (Figure 7(a)), but the control signal with the nonlinear controller was unacceptable at low speed (Figure 7(b)).
The low-speed sensitivity of the NL controller is in concordance with the theoretical considerations of Section 2.6, where inflation of the input sensitivity function gain |U NL | at low speed was clearly seen (cf. Equation (18) and Figure 2). Table 3. Low-speed evaluation interval: primary outcome measures for paired comparisons and p-values for comparison of means (see also Figure 6).  Table 3); the green lines link the sample pairs from each subject; the red horizontal bars depict mean values. D is the difference between the paired samples: D = L − NL. MD is the mean difference (red horizontal bar), with its 95% confidence interval (CI) in blue. Inclusion of the value 0 within the 95% CI for RMSE, (a), signifies a non-significant difference between the means (p > 0.05, Table 3); a significant difference for P v is signified by 0 lying outwith the 95% CI, (b), (p < 0.05, Table 3).

Discussion
The aim of this study was to systematically compare linear (L) and nonlinear (NL) heart rate controllers using quantitative measures of root-mean-square tracking error and average control signal power. The primary findings are, first, that there were no significant differences in the primary performance outcome measures (RMS tracking error and average control signal power) between the linear and nonlinear controllers at moderate-to-vigorous exercise intensities (i.e. over the main evaluation interval). Second, the nonlinear control strategy was seen to be highly sensitive to running speed, with higher, and often unacceptable, average control signal power at low speed. This behaviour can be readily understood by considering the small-signal gain of the square-root function, which is given by dv/dx = 1/(2v), Section 2.3. Thus, the small-signal gain of the NL strategy tends to infinity as treadmill speed tends to zero, since lim v→0 dv/dx = ∞. This results, in turn, in the inflation of |U NL | as described above.
The secondary finding of the study is that accurate, stable and robust performance was achieved for all 16 subjects with the linear controller: RMS tracking error was low; the average control signal power was sufficiently low to give an acceptably smooth control signal (treadmill speed), and it had low dispersion across all subjects; and robustness was demonstrated since a single approximate model was used that was not specific to any of the subjects tested. These findings support a previous observation that a well-designed, model-based linear controller can give satisfactory overall heart rate control performance (Hunt & Fankhauser, 2016). The nonlinear controller was also accurate and robust, but with the tendency described above for the control signal to be very sensitive at low speed.
The fact that no significant differences in the primary outcomes were seen at moderate-to-vigorous intensities, even with the inclusion of n = 16 subjects in the comparison sample, suggests that any differences (effect sizes) which do actually exist are small and have little relevance from a practical standpoint: taking the sample effect sizes and dispersions observed here as estimates for a putative population, statistical power calculations reveal that sample sizes of n = 93 for RMS tracking error and n = 60 for average control signal power would be required to demonstrate differences between the L and NL controllers (with an assumed significance level α = 0.05 and statistical power of 1 − β = 80%).
The lack of significant differences in outcomes between L and NL at moderate-to-vigorous intensities, together with the adequacy of a simple approximate linear model and the low-speed sensitivity of the nonlinear controller, provides further support to the concept put forth previously, (Hunt & Fankhauser, 2016), that achieving an appropriate response to VLF-HRV is the overriding challenge in the design of heart rate controllers, whereas issues of parametric and structural plant uncertainty are secondary. One direct way to design for disturbance rejection of VLF-HRV is to use input-sensitivity shaping (Hunt & Fankhauser, 2016), while the static nonlinearity compensation strategy at the core of the NL approach tested here does not directly consider closed-loop frequency responses.
Although the performance of the linear controller L was satisfactory, this can be further improved by a simple extension of the structure of the controller transfer function C. Here, this transfer function was purposely constrained to take the form of Equation (6), C(z −1 ) = (s 0 + s 1 z −1 )/(1 − z −1 ), in order that the L controller and the NL controller, as documented in Scalzi et al. (2012), would be identical at the chosen operating point v o . Technically, this amounted to constraining C(z −1 ) to be merely causal, and not strictly causal (strict causality requires a pure delay term z −1 as a factor in the numerator of C). Equivalently, C(z −1 ) written in terms of z, that is, C(z −1 ) = (s 0 z + s 1 )/(z − 1), was proper in z, but not strictly proper. This has an important practical consequence: the gain of a strictly proper compensator will roll off to zero at high frequency, as lim ω→∞ |C| = 0, thus protecting the loop from high-frequency disturbances, whereas the highfrequency gain of a proper compensator will tend to a finite, non-zero, value (Franklin et al., 2008). This property carries over to the input sensitivity function U o , wherefore the gain |U o (e jωT s )| does not roll off to zero for the compensator C employed here (Figure 2). Because of this, the control signal with L remains somewhat sensitive to highfrequency HRV disturbances, as was observed through the mean value for average control signal power for L, P v = 51.7 × 10 −4 m 2 /s 2 (Table 2). This value can be contrasted with that obtained in a previous 'unconstrained' study, (Hunt & Fankhauser, 2016), using a strictly proper compensator which gave an input sensitivity function with high-frequency rolloff to 0 (see |U sp | in Figure 2): the mean value across 30 subjects was P v = 16.0 × 10 −4 m 2 /s 2 and the control signal was qualitatively perceived to be substantially smoother than in the present work; the mean RMS tracking error was, however, higher -2.96 vs. 2.34 bpm, previous study vs. controller L herethus demonstrating the inherent tradeoff in feedback design between tracking accuracy and control signal power.
This study, and the previous, related work (Hunt et al., 2015;Hunt & Fankhauser, 2016), focussed on empirical testing with healthy male subjects running at moderate-to-vigorous intensities. Further work should be carried out to determine whether the present findings and observations carry over to other populations and at low (walking) and ultra-high intensities. Formalized exercise training programmes are applied systematically in cardiac rehabilitation; the current recommendations for exercise intensity in such patients are also given in terms of heart rate (Mezzani et al., 2013;Pescatello et al., 2014). Since the dynamic response of heart rate in this type of pathology is likely to be more complex and variable than in healthy individuals, (Cheng et al., 2008), careful investigation needs to be carried out to determine whether nonlinear control approaches bring demonstrable benefits in such cases. With nonlinear strategies, caution will need to be exercised at the low speeds employed in patient groups, because of the potential control signal sensitivity highlighted above.
Based on the results and interpretations of this study, a number of recommendations can be made regarding the conduct of future studies which investigate new heart rate control strategies. The performance of feedback control systems for heart rate is dominated by broadspectrum heart rate variability, the very-low-frequency band being the most perceptible to the runner. Moreover, the HRV is itself highly variable over time within a given subject, and between subjects. Because of this degree of variability, it is recommended that: • any new control strategy should be systematically tested against the performance of a well-designed linear controller (n.b., a PI controller whose gains are set using an empirical, trial-and-error procedure does not fall into this category); • controller evaluation should employ objective, quantitative performance outcomes for tracking accuracy and control signal intensity (e.g. RMS tracking error, RMSE, Equation (1), and average control signal power, P v , Equation (2), as used here); • evaluations and comparisons should give priority to experimental investigations with human test subjects, rather than focussing on simulation results (simulations are 'doomed to succeed' in this application, whereas the unpredictability and variability of HRV in practice introduces a large degree of uncertainty); • comparative evaluations should use a common subject cohort; • comparative investigations should employ carefully designed protocols and formal statistical analyses with sample sizes (n = number of subjects tested with each controller) sufficient to allow meaningful conclusions to be drawn; experience from the present work and from related reports indicates that sample sizes in the range of 10-20 are appropriate for detection of practically relevant effect sizes.

Conclusions
There were no significant differences in the primary performance outcome measures between the linear and nonlinear controllers at moderate-to-vigorous intensities, but the nonlinear control strategy was overly sensitive at low running speed. Accurate, stable and robust performance was achieved for all 16 subjects with the linear controller designed using a single, approximate plant model. These results provide further evidence that disturbance rejection of VLF-HRV is the overriding challenge in the design of heart rate controllers, whereas issues of parametric and structural plant uncertainty are apparently secondary. Future studies which investigate new heart rate control approaches should include comparative data from well-designed linear controllers, they should employ quantitative performance outcomes, and they should use experimental designs with a common subject cohort of appropriate sample size.

Disclosure statement
No potential conflict of interest was reported by the authors.