By Ahmed Ramadan A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering Doctor of Philosophy 201 8 ABSTRACT By Ahmed Ramadan Low back pain (LBP) and neck pain (NP) are two major musculoskeletal disorder s that cause disability . Changes in motor con trol are associated with LBP and NP. These changes can be revealed by designing motor control tasks to test healthy and pain groups. The very first research We evaluated the within - and between - visit reliability of a novel seated balance test that utilized a torque - controlled robotic seat . This test assesses trunk motor control under three conditions: eyes open (EO), eyes closed (EC) , and eyes closed with vibration to the lumbar muscles (VIB). Intraclass correlation coefficients and c oefficients of multiple correlation respectively. Therefore, our novel seated balance test is reliable. Once a test is deemed reliable, further investigations of motor control can be con ducted. First, we devised an approach that infers the control intent (motor control strategy) by solving inverse Model Predictive Control (iMPC) problem s using input/output data . Model Predictive Control (MPC) allows us to model physiological constraints i n motor control. For a subject who performed the seated balance test, t minimizing the squared sum of the upper - body and lower - body angles. Motor control utilizes several physiolog ical pathways: visual, vestibular, proprioceptive, and intrinsic. Parametric models of these pathways provide insight into the relative contributions of various pathways. However, estimating model parameters is challenging because (1) there are many parame ters to be estimated , (2) the information in experimental data is limited since subjects do not tolerate long testing periods before fatiguing, and (3) non - invasive external sensors are not available to measure all signals of interest. Therefore, we devise d a systematic method of selecting sensitive parameter subsets to be estimated , while fixing the remaining parameters to values obtained from preliminary estimation. Our method , based on the Fisher Information Matrix (FIM) , was applied on a rotational head position tracking test. O ur method (1) reduced model complexity by only requiring five out of twelve parameters to be estimated, (2) significantly reduced parameter variability ( 95% confidence intervals ) by up to 89% of the original confidence interval, ( 3) maintained goodness of fit (v ariance a ccounted f or ) at 82%, (4) reduced computation time, where our FIM method was 164 times faster than the existing nonlinear Least Absolute Shrinkage and Selection Operator (LASSO) method , and (5) selected similar sens itive parameters to the LASSO method, where three out of five selected sensitive parameters were shared by FIM and LASSO methods. We investigated the feasibility of incorporating parameter reliability and model diversity in our FIM method to narrow down th e solution space from sensitive parameters to key parameters of motor control. First, we incorporated parameter reliability as a third selection criterion besides parameter variability and goodness of fit. Then, we incorporated model diversity as a fourth optional criterion. We applied this approach to head position tracking in axial rotation and flexion/extension. selected three key parameters out of twelve. The key parameters were associated with at least two physiological pathways out of four modeled pathways which is a measure of model div ersity. The average goodness of fit was >69%. Copyright by AHMED RAMADAN 2018 v ACKNOWLEDGEMENTS I would like to acknowledge my Ph.D. committee; my advisor : Prof. Jongeun Choi; PIs : Prof. Jacek Cholewicki and Prof . N. Peter Reeves; co - PI : Prof. Clark Radc liffe, and co - investigator : Prof . John Popovich Jr. who has been a true friend. My sincere thanks extend to my small family: Aya, Yaseen, and Hamza; and my big family overseas : Sania, Ramadan, Aya, and Mohamed for their emotional and physical support. vi T ABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ....................... viii LIST OF FIGURES ................................ ................................ ................................ ........................ x CHAPTER 1 INTRODU CTION ................................ ................................ ................................ .... 1 CHAPTER 2 RELIABILITY OF ASSESSING POSTURAL CONTROL DURING SEATED BALANCING USING A PHYSICAL HUMAN - ROBOT INTERACTION ................................ 5 Abs tract ................................ ................................ ................................ ................................ ....... 5 2.1 Introduction ................................ ................................ ................................ ...................... 6 2.2 Methods ................................ ................................ ................................ ............................ 7 2.2.1 Subjects ................................ ................................ ................................ ..................... 7 2.2.2 Data Collection ................................ ................................ ................................ ......... 8 2.2.3 Data Analysis ................................ ................................ ................................ .......... 13 2.3 Results ................................ ................................ ................................ ............................ 18 2.4 Discussion ................................ ................................ ................................ ...................... 20 CHAPTER 3 INFERRING INTENT OF TRUNK MOTOR CONTROL IN SEATED BALANCE USING INVERSE MODEL PREDICTIVE CONTROL ................................ ............................. 23 Abstract ................................ ................................ ................................ ................................ ..... 23 3.1 Introduction ................................ ................................ ................................ .................... 24 3.2 Model Predictive Control (MPC) ................................ ................................ ................... 26 3.3 Inverse MPC (iMPC) ................................ ................................ ................................ ..... 27 3.3.1 Estimating the Starting Point ................................ ................................ ............. 30 3.4 iMPC ver sus iLQR Example ................................ ................................ .......................... 31 3.5 The Seated Balance Experiment ................................ ................................ .................... 38 3.5.1 The pHRI ................................ ................................ ................................ ................ 38 3.5.2 Subject ................................ ................................ ................................ ..................... 41 3.5.3 Control Intent ................................ ................................ ................................ .......... 42 3.6 Discussion ................................ ................................ ................................ ...................... 45 3.7 Conclusion ................................ ................................ ................................ ...................... 47 CHAPTER 4 SELECTING SENSITIVE PARAMETER SUBSETS IN DYNAMICAL MODELS WITH APPLICATION TO BIOMECHANICAL SYSTEM IDENTIFICATION ...................... 48 Abstract ................................ ................................ ................................ ................................ ..... 48 4.1 Introduction ................................ ................................ ................................ .................... 49 4.2 Methods ................................ ................................ ................................ .......................... 51 4.2.1 Subjects ................................ ................................ ................................ ................... 51 4.2.2 Data Collection ................................ ................................ ................................ ....... 52 4.2.3 Parametric Model ................................ ................................ ................................ .... 52 4.2 .4 Parameter Estimation ................................ ................................ .............................. 55 vii 4.2.5 Sensitive Parameter Selection ................................ ................................ ................. 57 4.3 Results ................................ ................................ ................................ ............................ 58 4.4 Discussion ................................ ................................ ................................ ...................... 64 4.5 Conclusion ................................ ................................ ................................ ...................... 66 CHAPTER 5 FEASIBILITY OF INCORPORATING TEST - RETEST RELIABILITY AND MODEL DIVERSI TY IN IDENTIFICATION OF KEY NEUROMUSCULAR PATHWAYS DURING HEAD POSITION TRACKING ................................ ................................ .................. 67 Abstract ................................ ................................ ................................ ................................ ..... 67 5.1 Introduction ................................ ................................ ................................ .................... 68 5.2 Methods ................................ ................................ ................................ .......................... 69 5.2.1 Subjects ................................ ................................ ................................ ................... 69 5.2.2 Data Collection ................................ ................................ ................................ ....... 70 5.2.3 Data Analysis ................................ ................................ ................................ .......... 70 5.3 Results ................................ ................................ ................................ ............................ 80 5.4 Discussion ................................ ................................ ................................ ...................... 81 APPENDICES ................................ ................................ ................................ .............................. 84 Appendix A Fisher Information Matrix (FIM) ................................ ................................ ......... 85 Appendix B Least Absolute Shrinkage and S election Operator (LASSO) ............................... 87 Appendix C FIM and LASSO Methods for Sensitive Parameter Selection ............................. 88 BIBLIOGRAPHY ................................ ................................ ................................ ......................... 94 viii LIST OF TABLES Table 2.1: Characteristics of the subjects presented as mean standard deviation. ..................... 7 Table 2.2: Be tween - visit reliability ( ) of seated balance test performance measures in time domain ( ). is given as mean ( , ). ................................ ....................... 19 Table 2.3: Between - visit reliability ( ) of seated balance test performance measures in frequency domain ( ). is given as mean ( , ). ................................ ................ 19 Table 2.4: Within - and between - visit reliability expressed as of the empirical tr ansfer function estimates (ETFE) from perturbation torque to lower and upper body angles ( and , respectively) in frequency domain. values are given as mean standard deviation. .......... 19 Table 2.5: across visits and conditions. Values are given as mean standard deviation. 20 Table 2.6: across visits and conditions. Values are given as mean standard deviation. . 20 Table 3.1: The proposed iMPC algorithm ................................ ................................ ................... 31 Table 3.2: Physiological parameters of the model in Figure 3.6 . They are divided in two groups separated by a double line. The upper group include fixed parameters. The lower one includes estimated parameters to fit the experimental data. ................................ ................................ ........ 44 Table 4.1: characteristics presented as mean standard deviation. ................................ ................................ ................................ ................................ ...... 51 Table 4.2: The neurophysiological parameters of the head position tracking task along with their physiological description. The lower and upper bounds were obtained from preliminary estimation ................................ ................................ ................................ 53 Table 4.3: Parameter estimates of the complex model with 12 estimated parameters (Nominal) a nd the reduced models with five sensitive parameters selected by FIM and LASSO methods ( and , respectively). Parameter 95% confidence interval width was estimated for the complex and reduced models. Goodness of fit was measured by Variance Accounted For ( ). All values are given as mean standard deviation across subjects. ................................ ......... 60 Table 5.1 : Demographic characteristics of the subjects given as (mean ± standard deviation). . 71 Table 5.2 : Physiological parameter s used in the neck motor control model ( Figure 5.2 ) and their bounds. Bounds were based on preliminary estimates ( Table 5.3 ). has model parameters that may be estimated (candidate key parameters), while has fixed parameters. .................... 72 ix Table 5.3 : The estimates of preliminary (Nominal) and selected key model paramete rs and the corresponding Variance Accounted for ( ) in Yaw. Values are given as mean ± standard deviation across subjects (and visits for Nominal). Between - visit reliability was measured by ICC(A, k) with 95% confidence interval (CI). ................................ ................................ ............. 76 Table 5.4 : The estimates of preliminary (Nominal) and selected key model parameters and the corresponding Variance Accounted for ( ) in Pitch. Values are given as mean ± standard deviation across subjects ( and visits for Nominal). Between - visit reliability was measured by ICC(A, k) with 95% confidence interval (CI). ................................ ................................ ............. 77 Table 5.5 : The sensitive subsets (solution space) of the FIM method ( Figure 5 .3 IntraClass Correlations (ICCs) are given in parentheses. The selected key model parameters for estimation are highlighted. ................................ ................................ ................................ ............ 7 9 Table C . 1 : The algorithm of sensitive parameter selection using the LASSO method. ............... 92 x LIST OF FIGURES Figure 2.1: The seated balance test utilizes a physical human - robot interaction (pHRI). (A) The u pright equilibrium position. (B) A tilted equilibrium position. ................................ .................... 8 Figure 2.2: (A) The two degree - of - freedom physical model of the seated balance test utilizing a pHRI. (B) Block diagram of the pHRI system. The input is a torque signal perturbation, whereas the output is the displacement vector (lower and upper body angles). - seat/subject plant dynamics, - seat stiffness, - seat damping, - total seat torque. ................... 10 Figure 2.3: The time series input - output data of a single trial. (A) The perturbation torque is a pseudorandom ternary sequence of 5 Nm amplitude. (B, C) The lower and upper body angles , respectively. ................................ ................................ ................................ .................. 11 Figure 2. 4: The power spectral density of (A) the perturbation , (B) the error , and (C) the error of a trial, respectively. ................................ ................................ ............... 15 Figure 2.5: The frequency domain data represented in ampli tude and phase of the ETFE for the same trial/subject/condition across three visits. The curves are plotted in the passband region of . (A, B) The ETFE from to . (C, D) The ETFE from to . .............................. 16 Figure 2.6: The frequency domain data represented in amplitude and phase of the ETFE for the same trial/condition/visit across three subjects. The curves are plotted in the passband region of . (A, B) The ETFE from to . (C, D) The ETFE from to . .............................. 17 Figure 3.1: (A) Block diagram of a generic closed - loop system utilizing MPC in the feedback. (B) Schematic layout of the two - level iMPC approach. ................................ ................................ ..... 28 Figure 3.2: (A) The system states ( ) using a full - state, feedback gain versus the reference states ( ). (B) The control ( ) using a full - state, feedback gain versus the reference control ( ). ................................ ................................ ................................ ............................... 33 Figure 3.3: (A) The iMPC states ( ) versus the reference states ( ). (B) The iMPC control ( ) versus the reference control ( ). ................................ ................................ ................ 36 Figure 3.4: (A) The iLQR states ( ) versus the reference states ( ). (B) The iLQR control ( ) versus the reference control ( ). The iLQR here is the modified gradient - descent algorithm in (3. 16 ). ................................ ................................ ................................ ....................... 37 Figure 3.5: (A) The pHRI for our seated - balance experiment in the upright equilibrium position. (B) The pHRI model of our seated - balance experiment. The variables with subscript are for the robot model and the those with subscript are for the human model. ............................. 40 xi Figure 3.6: Block diagram of the closed - loop system of trunk motor control during seated balancing. The experimental output measurements a re the time series of . ...................... 41 Figure 3.7: MPC trajectories (dashed curves) versus experimental trajectories (continuous curves). ................................ ................................ ................................ ................................ .......... 42 Figure 4.1: (A) Ex perimental setup for head position tracking. The reference marker is r(t) and the head position marker is y(t). (B) The block diagram of the head position tracking model adopted from [22], [89]. The transfer functions come from [22], [23] after accounting for the change in output from acceleration to position. ................................ ................................ ............ 54 Figure 4.2: The best - fit model output (during a trial) when estimating the FIM - and LASSO - selected sensitive parameters. The experimen .............. 55 Figure 4.3: An outline of model - based bootstrap [103] using 1000 replicates to estimate 95% confidence interval (CI) of the parameter vector . ................................ ................................ ...... 57 Figure 4.4: Results of applying the FIM method (Appendix C) with to the experimental data . The optimal number of sensitive parameters per subset is . The number of subsets (each with parameters) is . ................................ ................................ ..... 59 Figure 4.5: Comparing the 95% confidence interval (CI) between estimating ALL parameters, FIM - selected parameters ( ), and LASSO - selected parameters ( - estimate and the error bars are the mean CI (across subjects) for (A) , (B) , and (C) that were selected as sensitive parameters by both FIM and LASSO methods. ................................ .. 63 Figure 5.1 : The experimental setup of head position tracking tasks in (A) axial rotation (Yaw), and (B) flexion - extension (Pitch). ................................ ................................ ................................ . 71 Figure 5.2 : The block diagram of neuromuscular control for the head position tracking tasks. The contributing pathways are visual, vestibular, proprioceptive, and intrinsic. ................................ 72 Figure 5.3 : The solution space (the sensitive parameter subsets) of the FIM method for (B) Yaw, and (B) Pitch. The solution space was discretized at 10 values of (the upper bound of parameter variability). is the maximum number of sensitive parameters per subset. is the number of sensitive subsets w ith parameters and parameter variability . ................................ ................................ ................................ ................................ 78 1 CHAPTER 1 I NTRODUCTION Low back pain (LBP) and neck pain (NP) are two common musculoskeletal disorder s that are ranked first and fourth , respectively, among causes for ye ars lived with disability [1] . Between 1990 and 2010, years lived with disability due to LBP and NP hiked by 42.6% and 41.0%, respectively [1] . Considering the rise in pain, it is no surprise that patients find conservative treatment methods, including pain medication and physical therapy, unsatisfactory [2] . Therefore, patients seek alternative treatment f or their pain. One - third of adults used some form of complementary and alternative medicine (CAM) [3] . Osteopathic manual therapy (OMT) is a CAM therapy that is being chosen by an increasing number of LBP and NP patients. OMT is suggested to be effective in relieving musculoskeletal pain [4] [6] . Subjective tests such as manual range of motion estimates by therapists, and survey statements provided by patients were used in most studies to support these findi ngs. Therefore, these research efforts are not sufficient to accurately determine the effectiveness of the therapy, characterize its mechanisms, or optimize treatment [7] . Therefore, the overall goal of this work is to develop objective tools for the ass essment of motor control of the lumbar and cervical spine systems. Changes in motor control are associated with LBP [8], [9] and NP [10] [12] . LBP patients exhibit degraded performance in postural control tasks, which can be quantified by the movement of the center of pressure (CoP) while st anding on a force plate [13] [15] or balancing on an unstable seat 2 [16] . NP patients showed diminished accuracy compared to healthy individuals in terms of the capacity to relocate the head in space [10] , which is believed to be linked to a reduction in somatosensory resolution [10], [17] . Researchers design motor control tasks to potentially reveal differences between healthy and p ain task (i.e., experimental data is reproducible). Reliability is often considered for input/output performance measures such as the root mean square error (R MSE) and the empirical transfer function estimates (ETFE) [18 ], [19] . After verifying reliability, the physiological pathways of motor control during the task can be investigated. Neuromuscular motor control has several physiological pathways that generate a neural drive to activate muscles during motor control tas ks. These pathways are visual, vestibular, proprioceptive, and intrinsic. The visual pathway responds to the position of the head in space during disturbance rejection tasks, whereas vision and cognition responds to the tracking error on a display during t racking tasks. Vestibular pathway responds to head acceleration during movements. Proprioceptive pathway responds to muscle fiber length and velocity changes (sensed by muscle spindles) and to the generated muscle force (sensed by the Golgi tendon organ). Intrinsic pathway has a passive spring/damper - like action. Parametric models for these physiological pathways were developed, in the literature, using the principles of control systems [20] [24] . Using system identification techniques, these models can be identified (e.g., model parameters were estimated) to gain insigh ts into key changes of motor control before and after OMT. Furthermore, when the system is performing poorly, this approach provides insight into possible mechanisms of impairment. 3 However, some difficulties are associated with these models, which limit th eir application. Difficulties occur mainly because these models can produce a large number of parameters that need to be estimated [25] . At the same time, the information in experimental data is limited because subjects do not toler ate a long data collection process before fatiguing, and non - invasive external sensors are not available to measure all signals of interest. To address these difficulties, this thesis makes the following contributions: In chapter 2, w e investigated the non parametric reliability of a novel seated balance test that utilized a torque - perturbed seat. T he seated balance test produced reliable performance measures . When comparing with other studies, there is some evidence to suggest that the physical human - robot interaction ( pHRI ) with a deterministic input (perturbation torque) helps improve reliability. In chapter 3, w e devised an inverse Model Predictive Control (iMPC) algorithm to infer the control intent of subjects while imposing physiological constraints (m aximum human torque). The novelty of this work stems from two aspects. First, there has been no previous research on iMPC in the literature to our knowledge. We proposed an iMPC algorithm to solve for the MPC cost function weights using experimentally coll ected output trajectories. Second, we imposed physiological constraints, which cannot be imposed with LTI controllers and their inverse optimal problems. In chapter 4, w e devised a systematic method of selecting sensitive parameter subsets to be estimated while fixing the remaining parameters to preliminary values . Our method was based on Fisher information matrix (FIM ) and was compared against the n onlinear Least Absolute Shrinkage and Selection Operator (LASSO) . Three out of five selected parameters were identical between both methods. T he 95% confidence intervals (CI) of the 4 three common parameters were reduced significantly versus the original CI . Goodness of fit was maintained at 82%. T he average computation time of our FIM method was less than 1% of th at in LASSO method In chapter 5, w e incorporated parameter reliability and model diversity in our FIM method to select key parameters affecting neck motor control . Including parameter reliability in system identification is the novel contribution in Chapte r 5. Using variability, goodness of fit, reliability, and diversity selection criteria on data of forty subjects, we selected one subset of key model parameters. This work has been facilitated by the MSU Center for Orthopedic Research (MSUCOR). MSUCOR is a state - of - the - art biomechanics research facility that merges the gap between clinical and basic research. This work was made possible in part by grant number U19 AT006057 from the National Center for Complementary and Integrative Health (NCCIH) at the Nati onal Institute of Health. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of NCCIH. 5 CHAPTER 2 RELIABILITY OF ASSES SING POSTURAL CONTRO L DURING SEATED BALANC ING USING A PHYSICAL HUMAN - ROBOT INTERACTION T his chapter is published at the Journal of Biomechanics [26] . Abstract This study evaluated the within - and between - visit reliability of a seated balance test for quantifying trunk motor control using input - output data. Thirty healthy subjects performed a seated balance test under t hree conditions: eyes open (EO), eyes closed (EC) , and eyes closed with vibration to the lumbar muscles (VIB). Each subject performed three trials of each condition on three different visits. The seated balance test utilized a torque - controlled robotic sea t, which together with a sitting subject resulted in a physical human - robot interaction (pHRI) (two degrees - of - freedom with upper and lower body rotations). Subjects balanced the pHRI by controlling trunk rotation in response to pseudorandom torque perturb ations applied to the seat in the coronal plane. Performance error was expressed as the root mean square ( ) of deviations from the upright position in the time domain and as the mean bandpass signal energy ( ) in the frequency domain. Intra - class correlation coefficients ( ) quantified the between - visit reliability of both and . The empirical transfer function estimates (ETFE) from the perturbation input to each of the two rotational outputs were calculated. Coefficients of multiple correlation ( ) quantified the within - and between - visit reliability of the averaged ETFE. of and 6 for all conditions were . The mean within - and between - visit were all for the lower body rotation and for the upper body rotation. Therefore, our seated balance test consisting of pHRI to assess coronal plane trunk motor control is reliable. 2.1 Introduction Several sensory systems are involved in maintaining upright postural control and balance: visual, vestibular, and pro prioception [20], [21], [27] . To determine the contribution of di fferent sensory feedback pathways to trunk postural control, experiments are performed with sensory information either removed (e.g., closed eyes to remove visual information) [16], [20] or corru pted (e.g., added vibration to reduce the quality of proprioceptive information) [28] [30] . Identifying the contributions of selected feedback pathways to trunk postural control is important for understanding the mechanisms of control and possible impairments. To that end, several studies have investigated the role of the trunk muscles in main taining postural control using sitting balancing tests [24], [31] [38] . Sitting posture has been selected, because the lower body is fixed so that ankle [39] , knee, or hip [40] stabilization strategies are not possible; the mo vement is limited to the trunk to maintain balance. For the seated balance test to be useful, it must reliably assess trunk control. Three studies addressed the reliability of seated balance tests, and the reported results range from poor to excellent [41] [43] . These inconsistencies could stem from the fact that unperturbed postural sway was used as a measure of performance. In such a case, the postural sway (output signal) is transient behavior contaminated with neuromuscular noise, which is not repeatable. Thus, data reproducibility and test reliability may vary from one study to another. To over come this problem, we developed a physical human - robot interaction (pHRI) system [31] that consists of a compliant 7 unstable seat and, at the same time, provides torque perturbation as a known input to the system [24], [44] [46] . This pHRI also allows for the precise regulation of test difficulty by controlling rotational seat stiffness. The purpose of this study was to assess th e within - and between - visit reliability of a seated balance test under three conditions: eyes open (EO), eyes closed (EC), and eyes closed with lumbar muscle vibration (VIB). The reliability was evaluated using both time and frequency domain analysis, whic h capture the control attributes such as accuracy and speed. This study uses nonparametric analysis of input - output data, therefore, the contributions of feedback pathways (e.g. , reflexive , and intrinsic) cannot be distinguished in this framework. Because the pHRI involved with this test was designed with a known input perturbation, we hypothesized that the performance measures in our seated balance test would be reliable. 2.2 Methods 2.2.1 Subjects Thirty subjects participated in this study ( Table 2 . 1 ). They were in self - reported good general health with no history of back pain lasting longer than three days or any neurological condition that could affect motor control. Subjects were instructed to wear their corrective lenses if their eyesi ght was impaired. Prior to testing, all subjects signed an informed consent. Michigan State Table 2 . 1 : Characteris tics of the subjects presented as mean standard deviation. Females (n=19) Males (n=11) Height (m) Weight (kg) Age (years) 8 2.2.2 Data Collection A pHRI was used in a seated balance test to assess trunk motor control in the coronal plane ( Figure 2 . 1 ). During the test, the robotic platform applied torque perturbations to the subjects. They were instructed to balance on the platform, with their arms crossed on the chest, using trunk motions in the coronal plane. (A) (B) Figure 2 . 1 : The seated balance test utilizes a physical human - robot interaction (pHRI). (A) The upright equilibrium position. (B) A tilted equilibrium position. The pHRI consisted of a two degrees - of - freedom system (a one degree - o f - freedom compliant - robot platform and a sitting subject) ( Figure 2 . 2 A). The first degree - of - freedom was the rotation and denoted by (lower - body angle). The second degree - of - 9 (trunk and above) rotation assumed to be about the fourth lumbar ( ) vertebra and denoted by (upper - body angle). A torque r (C062C - 13 - 3305, Kollmorgen, Radford, VA) was applied to the seat about ( Figure 2 . 2 A). This torque contained components that are (1) a known perturbation , and (2) a stiffness - and - damping torque ( Figure 2 . 2 B ). The perturbation signal was a pseudorandom ternary sequence ( Figure 2 . 3 A). The lower body angle (the seat angle about ) was measured by an encoder embedded in the moto r. The angle was constrained to be within ± 12 º using mechanical stops. Our preliminary data demonstrated that the encoder accurately reflected the lower body angles. The RMS difference between the pelvis angles, captured independently with an optoelectroni c motion system, and the encoder was 0.62 degrees. The upper body angle was measured by two string potentiometers (SP2 - 50, Celesco, Chatsworth, CA) attached to a chest harness and a belt at the fourth thoracic ( ) and the levels, respectively ( Figure 2 . 1 B). The time series data of both angles is shown in Figure 2 . 3 B - C for a single trial. Each subject performed balancing tests under three conditions: EO , EC, and VIB. For each condition, subjects performed seven trials (four 30 - second practice trials and three 30 - second full trials). The last three (full) trials were used for analysis in this study. According to Larivière et al. (2013) , three trials provi ded acceptable reliability. 10 (A) (B) Figure 2 . 2 : (A) The two degree - of - freedom physical model of the seated balance test utilizing a pHRI. (B) Block diagram of the pHRI system. The input is a torque signal perturbation, whereas the output is the displacement vector (lower and upper body angles). - seat/subject plant dynamics, - seat stiffness, - seat damping, - total seat torque. P + - d(t) Robot Seated Subject Body Angles, y(t) Perturbation T orque 11 (A) (B) (C) Figure 2 . 3 : The time seri es input - output data of a single trial. (A) The perturbation torque is a pseudorandom ternary sequence of 5 Nm amplitude. (B, C) The lower and upper body angles , respectively. To normalize test difficulty between subjects, both sea t stiffness ( ) and perturbation torque set to so that the seat was stable even when there was no subject. First, the sea t stiffness was normalized. The subject was instructed to balance on the seat with the stiffness provided from the electric motor, but without torque perturbations. Beginning from a high 12 stiffness value ( 250 Nm/rad), the subject was asked to maintai n balance. If the subject successfully maintained balance (i.e., seat angle stayed within ± 12 º [47] ), the operator gradually reduced the s tiffness (similar to [47] ) by 10 Nm/rad decrements until the subject could no longer maintain balance. This level of stiffness is called c ritical stiffness [47] . The procedure was repeated once again and the lowest critical stiffness was chosen as a coarse estimate of the sub stiffness. For fine estimate of critical stiffness, the procedure was repeated three times starting from the coarse estimate of critical stiffness +50 Nm/rad. In these repetitions, the stiffness was reduced by 5 Nm/rad decrements until the subject could no longer maintain balance. The average of the three results was selected as a fine estimate of the critical stiffness. Thereafter, the experimental seat stiffness was set to twice the fine critical stiffness. Second, the perturbation torque amplitude was normalized. With the experimental stiffness provided from the electric motor, a square - wave pulse torque with 1 Nm amplitude was applied to the seat for 1 s and then released. The pulse torque direction was random so that the subject could no t predict it. If the seat angle did not reach ± 8 º after the pulse, the operator increased the amplitude by 1 Nm and the pulse torque was repeated until the seat angle reached approximately ± 8 º . This result was a coarse estimate of the perturbation torque a mplitude. Fine normalization of the torque amplitude was performed because the experimental trials used an input signal that was not a simple step input as the one used in the coarse normalization. Two normalization trials were performed for this purpose w ith an input signal that was very similar to the experimental trial signal with the maximal perturbation torque amplitude equal to the coarse estimate. The duration of a normalization trial was 30 s. The subject performed these trials with 13 eyes closed and arms crossed. If the subject lost balance (seat angle = ± 12 º identified with electronic switches) during a normalization trial, the perturbation torque amplitude was decreased by 1 Nm and the trial was repeated. If the subject passed two normalization tria ls, the highest passing perturbation torque amplitude was selected to be used in the experimental trials. All data was logged using a data acquisition board (PCI - DAS6036, Measurement Computing, Norton, MA) and a quadrature encoder board (PCI - QUAD04, Measur ement Computing, Norton, MA). The data was sampled at a rate of 1000 samples/s. Testing was conducted on three different visits, separated by a minimum of 23 hours. The test order was consistent between visits: (1) EO, (2) EC, and (3) VIB. Rest periods (ap proximately 30 s) were given between all trials. 2.2.3 Data Analysis A simplified block diagram of the closed - loop system representing the seated balance test that utilizes pHRI is shown in Figure 2 . 2 B. The plant is a dynamical sys tem that models the two degree - of - freedom seat/subject dynamics including the gravity effect ( Figure 2 . 2 A). Trunk motor control generates lumbar torque about and the robot generates total torque about . The perturbation input is . The overall pHRI performance was characterized by two measures. First, the performance error, in the time domain, was expressed a s a root mean square ( ) of deviations from the upright position (both angles equal zero). Since the time - domain data includes high - frequency noise, it is beneficial to consider the error in the frequency domain within a passband frequency region to whi ch motor control mechanisms contribute. Therefore, the second performance measure was the mean bandpass signal energy of the error ( ) in the frequency domain. 14 The is a common performance measure in the time domain to assess motor control [18], [19], [43], [48] . The error or performance error ( ) used in the past studies refers to the scalar output assuming that subjects try to minimize it during a motor control test. For example, Popovich et al. (2015) and Reeves et al. (2014) used models with a single output, therefore, the scalar output was used as the performance error . In our model, there are two outputs . Thus, the motor control performance may be quantified by two errors, i.e. , and . Frequency response was used to obtain additional performance measures of the test. From the frequency - domain data of the error , we selected mean bandpass signal energy ( ) as the second performance measure. It was introduced in [19] as , ( 2 . 1 ) where is the power spectral density of the perfor mance error at frequency , is the frequency resolution, and and are the indices of the lowest and highest frequency points in the passband region, respectively. The passband region was defined as the largest contig uous frequency range containing >3% of the maximum input power spectral density ( Figure 2 . 4 A) [19] . For these tests, the passband region was . The error power spectral density of a trial is depicted in Figure 2 . 4 B - C. The summation represents the area under the curve. Intra - class correlation coefficients [49] along with confidence intervals ( ) were computed as a measure of the reliability between the three visits. was computed for and using the mean of the last three trials of each visit. was classif ied as fair ( ), good ( ), and excellent ( ) [50] . The standard error of measurement ( ) was calculated as . 15 Coefficient of multiple correlations ( ) was used to quantify within - and between - visit reliability. was modified from [51] to use the frequency - domain data for the passband region. For this purpose, the empirical transfer function estimates (ETFE) [52], [53] were obtained from to ( Figure 2 . 5 A - B and Figure 2 . 6 A - B) and from to ( Figure 2 . 5 C - D and Figure 2 . 6 C - D).The three - dimensional ETFE (real, imaginary, and frequency dime nsions) was used instead of the time - domain data used in [51] . (A) (B) (C) Figure 2 . 4 : The power spectral density of (A) the perturbation , (B) the error , and (C) the error of a trial, respectively. 16 (A) (C) (B) (D) Figure 2 . 5 : The frequency domain data represented in amplitude and phase of the ETFE for the same trial/subject/condition across three visits . The curves are plotted in the passband region of . (A, B) The ETFE from to . (C, D) The ETFE from to . 17 (A) (C) (B) (D) Figure 2 . 6 : The frequency domain data represented in amplitude and phase of the ETFE for the same trial/condition/visit across three subjects . The curves are plo tted in the passband region of . (A, B) The ETFE from to . (C, D) The ETFE from to . Performance improvement across the three visits was assessed using a repeated - measures ANOVA. P - values were considered significant. If t he p - value for a certain condition (EO, EC, VIB) was significant, a post - hoc paired t - test (without corrections) was used for pairwise comparisons between visits. 18 2.3 Results The of the performance measures presented in Table 2 . 2 and Table 2 . 3 are , which reflect excellent between - visit reliability, however, associated with were better than that of . The confidence intervals have small ranges as well. The relatively smaller i n the case might be a result of the averaging in the ETFE computation. With respect to the overall characteristics of pHRI behavior in the frequency - domain, of within - visit reliability for ETFE were between and , while of between - visit reliability ranged from to ( Table 2 . 4 ). Both reflect excellent reliability, however, associated with were better than that of . Comparing subject's performance over the three visits indicates th at some learning in the selected conditions and measures may have occurred ( Tabl e 2 . 5 and Table 2 . 6 ). The time domain performance measure ( ) improved across the three visits in the EO and VI B conditions (p=0.05 and 0.02, respectively), whereas the frequency domain performance measure ( ) improved across the visits for only the VIB condition (p=0.05). Neither the time nor frequency domain performance measures for the EC condition c hanged significantly across the three visits (p>0.05). and did not have significant changes either. The post - hoc pairwise comparisons revealed that was significantly smaller in visits 2 and 3 in comparison to visit 1 for the VIB condition. was significantly smaller in visit 3 in comparison to visit 1 for the VIB condition as well. For the EO condition, of visit 2 was significantly smaller than that of visit 1. 19 Table 2 . 2 : Between - visit reliability ( ) of seated balance test performance measures in time domain ( ). is given as mean ( , ). Condition EO EC VIB Table 2 . 3 : Between - visit r eliability ( ) of seated balance test performance measures in frequency domain ( ). is given as mean ( , ). Condition EO EC VIB Table 2 . 4 : Within - and between - v isit reliability expressed as of the empirical transfer function estimates (ETFE) from perturbation torque to lower and upper body angles ( and , respectively) in frequency domain. values are given as mean standard deviation. Condi tion EO EC VIB 20 Tabl e 2 . 5 : across visits and conditions. Values are given as mean standard deviation. Condition P - value P - value Visit 1 Visit 2 Visit 3 Visit 1 Vis it 2 Visit 3 EO 0.05 0.29 EC 0.14 0.66 VIB 0.02 0.71 # indi Table 2 . 6 : across visits and conditions. Values are given as mean standard deviation. Condition P - value P - value Visit 1 Visit 2 Visit 3 Visit 1 Visit 2 Visit 3 EO 0.06 0.37 EC 0.44 0.57 VIB 0.05 0.70 # 2.4 Discussion The objective of this study was to determine if seated bal ance performance is reliable within and between testing sessions. Three conditions (EO, EC, and VIB) that could be used to examine various sensory pathways involved in postural control were assessed for reliability. The results of this study show that the seated balance test produced reliable performance measures under the three conditions, both in the time - domain ( ) as well as the frequency - domain ( ). When comparing with other studies, there is some evidence to suggest that the pHRI 21 helps improve reliability. For instance, physical springs were used to adjust the rotational stiffness of a balance seat [48] . Their reliability was poor to fair (ICC ranged 0.29 - 0. 56). The difference in reliability between that study and the present study may be explained by the presence of a perturbation signal. With the pHRI and the use of a perturbation signal, the response of trunk control, reflecting the system output, will be correlated to the input signal. This correlation improves the measurement reproducibility, which in turn improves reliability. For a seated balance test without a perturbation input, the output of the system is dominated by random noise generated within th e neuromuscular system. Other studies that performed seated balance on a hemisphere without a perturbation signal had lower reliability (within - visit ICC range 0.03 - 0.70) [43] than our pHRI with perturbation input, which supports the u se of a pHRI along with perturbations when assessing trunk control. Another advantage of the pHRI is the ability to easily regulate test difficulty, which also appears to be important for producing reliable performance measures of trunk control [48] and identifying differences between people with and without back pain [16] . Other seated balance set - ups have the ability to regulate test difficulty by changing the radius of the hemisphere [16], [41], [43], [54] or by changing the position of the springs used to provide rotational stiffness [32], [35], [42], [48] . This adjustment, however, often requires the subject to be removed from the seat, and in the case of the hemisphere, has less resolution in changing the test difficulty. With the pHRI set - up, there is greater resolution i n changing test difficulty and seat stiffness can be adjusted with the subject on the seat. The results showed that some learning took place in the EO ( ) and VIB conditions ( ), which suggest that more practice may be required prior to 22 assessing trunk control. Future work should investigate how much practice is needed to converge on relatively stable performance. Moreover, learning may be due to modulations in the contributing pathways [46] (e.g. , intrinsic and vestibular pathways in the VIB condition). Future work should investigate such pathways by using parametric models. One of the limitations of the study is that seat st iffness and torque perturbations were not normalized to the subject anthropometrics. The seat stiffness and the perturbation amplitude were determined based on the subject's performance during the normalization phase of data collection (i.e., critical seat stiffness and maximal perturbation amplitude at which the subject could maintain balance). Because the seated balance performance measures were reliable within and between visits, the current normalization protocol is suitable for repeated measures testin g, such as assessing changes pre - post intervention. However, it is unclear if this normalization protocol could be used to assess between - group differences in trunk control, such as comparing trunk control between people with and without back pain. To allo w for such comparisons, it may be necessary to normalize test difficulty to subject anthropometrics (i.e., dynamics from body size) and not overall balance performance (i.e., dynamics from body size and control capabilities). Future work should investigate various methods for normalizing test difficulty to anthropometrics, thus isolating the control element and allowing pHRI to be used more broadly to study trunk control in various populations. 23 CHAPTER 3 INFERRING INTENT OF TRUNK MOTOR CONTROL IN SEATED BALANCE USI NG INVERSE MODEL PRE DICTIVE CONTROL Most of this chapter is published in the proceedings of the 2016 American Control Conference [55] . Abstract To understand differences in control intent between individuals with and without Lower Back Pain (LBP), it is important to develop techniques to esti mate motor control strategies. In this paper, we present an approach that infers the control intent (motor control strategy) by solving an inverse Model Predictive Control (iMPC) problem. The Model Predictive Control (MPC) structure models constraints, so it allows us to model the physiological constraints in trunk motor control. We devised an iMPC algorithm to solve iMPC problems with experimentally collected output trajectories. We tested the algorithm with a simulation example and the algorithm could ret rieve the true MPC weights. Then, we applied our algorithm on experimental data of one subject. Our experiment was a seated balance test that used a physical Human - Robot Interaction (pHRI). Results show that the estimated MPC weights reflected the task ins tructions given to the subject intent was dominated by minimizing the squared sum of the upper - body and lower - body angles. 24 3.1 Introduction Assessment of trunk mot or control is crucial for understanding basic spine function and coping mechanisms for those suffering Low Back Pain (LBP) [56] . The assessment is useful for evaluating differences in trunk motor control between those w ith and without LBP. This information can then be used to guide treatment and allows for monitoring the effects of treatment and rehabilitation [57] . For example, people with LBP often exhibit higher levels of trunk muscle coactivation [58], [59] , possibly re presenting a protective coping strategy [60] . From now on, the trunk motor control strategy will be referred to as the control i ntent [61] . Intent analysis is an emerging technique in science and engineering [61] [65] . Particularly, intent analysis is employed to assist physical Human - Robot Interactions (pHRI) or physical Human - Computer Interactions (pHCI) in understanding and predicting hu man behavior in a specific application [66] . Huang et al. [65] employed intent inference in autonomous vehicles to determine if nearby vehicles are changing lanes. Warrier and Devasia [67] proposed an iterative - learning - control appr oach to infer the intent of a novice user in tracking tasks. Farivar et al. [68] proposed intelligent controls to control eye movement in pet - type robot for social children - robot interactions . We devise an approach to infer the control intent. The long - term goals are to investigate differences in control intent between people with and without LBP and to identify subgroups of LBP with similar control intent, which can then be used to individualize treatment. We propose to capture control intent by employing a control theory approach. A feasible approach yields a solution whose simulated trajectories fit the experimental ones well, and is capable of imposing physiological const raints (e.g., maximum lateral torque of the trunk). It is common in the literature to assume 25 that human performance is optimized in different tasks [69] [72] . Therefore, we assume that the measured controlled behavior represents this optimized control. From a mathematical perspective, motor control can be modeled as an optimal controller, but the exact formulation such as the cost function is unknown [70] . On the other hand, the closed - loop system response is known from experimental measurements. Then, the inverse optimal control problem solves for the cost function knowing the measured response. However, the cost function structure must be prescribed beforehand in terms of weighted basi s functions. Thus, the inverse problem solves for the weights that reflect the relative significance of the different dynamic variables during the control task. For example, if a subject intends to penalize one variable more Two optimal control approaches were used in this study: Linear - Quadratic Regulator (LQR) and Model Predictive Control (MPC). Although MPC has been utilized in many engineering studies (e.g. [73], [74] ), it has not been utilized, to our knowledge, in biomechanical applications. The inverse problems of LQR and MPC are denoted iLQR and iMPC, respectively. An iLQR algorithm was developed in [71] and was applied to infer the control inten t, however, physiological constraints were not considered. Because MPC considers constraints [75] , we used iMPC to model the physiological constraints in trunk motor control. Our contributions are as follows. First, we propose an iMPC algorithm to solve for the MPC cost function weights using experimentally collected output traj ectories. We do not use reference control trajectories as in [70], [72], [76] , since the control is not available for measurement in our application. Second, we impose physiological constraints, which cannot be imposed with LTI 26 con trollers [24], [31] and their inverse optimal problems [71] . We show that our iMPC approach outperforms the iLQR approach in a simulation example where the reference trajectories were generated from an unstable plant controlled by MPC. Third, from the collected experimental data of the seated balance test. Preliminary work was presented in [55] without the experimental validation. 3.2 Model Predictive Control (MPC) The MPC structure is based on the state space representation of linear discrete - t ime systems. We only consider the single - input case throughout the paper. A block diagram of the closed - loop system, which uses MPC in the feedback, is shown in Figure 3 . 1 A. The standard regulating MPC problem takes the form [75] ( 3 . 1 ) where , , , and is the optimizer. is the prediction horizon, and is the terminal set. and are penalty matrices, where . It is assumed that ( 3 . 2 ) The MPC must stabilize the closed - loop system since the subject achieved the task in a stable manner. To guarantee the closed - loop system stability, we design the terminal penalty matrix 27 and the terminal set according to Theorem 4.3 in [75] as follows. First, an LQR control gain is designed using the pair , ( 3 . 3 ) where is the discrete - time algebraic Riccati equation (ARE) solution, that is ( 3 . 4 ) Then, the terminal penalty matrix is . Second, is chosen to be the maximal positive invariant set of the closed - loop system , subject to the imposed bounds and . We used MATLAB Control Toolbox along with MPT3.0.20 toolbox [77] to obtain the maximal invariant set. The MPC problem was implemented using YALMIP [78] . We selected t he quadprog 3.3 Inverse MPC (iMPC) We propose to solve the iMPC problem in a two - level approach as shown in Figure 3 . 1 B. In the lower level, the forward MPC problem i s solved for some starting - point weights . The optimal control action is then fed to the plant to obtain the simulated output . In the upper level, the weights are updated to minimize the error between the simulated an d the reference (experimental) outputs. The new weights are then passed to the lower level and so on 28 until the weights converge to their optimal values. The iMPC problem takes the form ( 3 . 5 ) (A) (B) Figure 3 . 1 : (A ) Block diagram of a generic closed - loop system u tilizing MPC in the feedback. (B ) Schematic layout of the two - level iMPC approach. Plant MPC (3.1) Disturbance Upper level: Use to solve (3.6) and update Lower level: Use to solve the MPC problem in (3.1) and obtain Close the loop with to obtain 29 where is the simulated output, and is the reference (experimental) output. is the simulation horizon and is the number of the output channels. stacks the columns of a matrix into a single vector, and returns the 2 - norm. computes the condition number ( ) of a matrix. is a scaling factor to scale the relatively large optimization is not dominated by minimizing . We assigned weight to the output - fitting term ( ) and to the condition - number term . Minimizing the condition number addresses the uniqueness problem of the solution [71] . is an estimated good starting point for . This formulation of the iMPC problem has a manifold of solutions. For example, if the optimal solution triplet is multiplied by a constant , the new triplet is also an optimal solution. To eliminate these equivalent replicates, normalization was proposed, for instance by fixing one element in the solution [70], [76 ], [79] . In our case, we normalized the solution by choosing . Thus, the inverse problem estimates only the pair . The iMPC problem becomes ( 3 . 6 ) To impose the constraint , a lower Cholesky de composition was used. It is given by , where is a lower triangular matrix. Thus, the constraint was imposed , where returns the diagonal elements of a matrix. On the other hand, the off - diagonal elements were not constraine d. 30 3.3.1 Estimating the Starting Point The iMPC problem is difficult to solve without a good initial estimate because it is a costly two - level optimization. To estimate a good initial point, first, we solved the iLQR problem to obtain the LQR weight matr ices [71] . Then, we set the LQR weights as initial - estimate MPC weights, with the MPC terminal weight matrix set to the discrete - time ARE solution. This initial estimate is the exact solution if no constraint is active [80] . To solve the iLQR problem, first, the linear control gain is required [71] . Therefore, the control was assumed to be a full - state feedback gain . It was estimated from the experimental data using nonlinear least squares optimization. Solving the iLQR problem was done using a linear matrix inequality (LMI) formulation or gradient - based least - squares minimization (when the LMIs were infeasible) [71] . If the LMIs were feasible, let the solution be , thus the starting point is If the gradient descent algorithm was used, let the solution be . In t his case, we formulated an auxiliary optimization problem to maximize the fitting of the estimated LQR gain and minimize the condition number as follows ( 3 . 7 ) where is an estimated good starting point for iMPC, returns the discrete - time LQR gain. is the same scal ing factor defined under (2). The iMPC algorithm is illustrated in Table 3 . 1 . 31 Table 3 . 1 : The proposed iMPC algorithm Input: (1) The system discrete - time state space model , where is for a scalar MPC control (2) The reference /experimental output Output: (1) The MPC weights 1: Assume the controller is a full state feedback gain , and estimate its value . 2: So lve the iLQR problem with the LMI formulation to obtain [19]. 3: if the problem is infeasible then 4: Solve the iLQR gradient descent problem [19] 5: Solve the auxiliary problem in ( 3. 7 ) to obtain . 6: else 7: Use 8: end if 9: Solve the iMPC problem in ( 3 . 6 ). 3.4 iMPC versus iLQR Example A simulation example was selected to have similar characteristics as the experiment. First, the passive spine system is unstable [33] . Therefore, the simulation example was selected as an unstable discrete - time system from [80] and given by ( 3 . 8 ) Second, the trunk mot or control is unknown, but it can be assumed to follow a control law such as PID [81] , static gain [24], [31] or LQR gain [71] . We assume it is a regulating MPC controller 32 with the following weights, ( 3 . 9 ) The weights were selected as identity so that condition number is minimal. Since we knew the true weights, we computed the residual that measures the error between estimated matrices and the true ones. The residual is given by ( 3 . 10 ) Third, the motor control is subject to physiological constraints such as the maximum lat eral bending torque [82] . Thus, we introduced control - input constraints given by ( 3 . 11 ) The closed - loop system was simulated with the initial state for 20 time steps and the horizon was 3 time steps. The simulated state trajectories served as our reference output data, i.e. . The control action saturated during the simulation as shown in Figure 3 . 2 B ( ). 33 (A) (B) Figure 3 . 2 : (A) The system st ates ( ) using a full - state, feedback gain versus the reference states ( ). (B) The control ( ) using a full - state, feedback gain versus the reference control ( ) . Given the reference output, we first estimated a full state feedback controller gain using nonlinear least squares. The estimated gain was 34 ( 3 . 12 ) The state and control fittings are shown in Figure 3 . 2 . Solving the iLQR problem in the LMI formulation was infeasible. Thus, the gradient descent iLQR problem was solved. The results were ( 3 . 13 ) Next, the auxiliary problem in ( 3. 7 ) was solved with the scaling factor . We obtained ( 3 . 14 ) Solving the iMPC problem in ( 3. 6 ) resulted in ( 3 . 15 ) Our iMPC algorithm estimated the MPC weights accurately and the corresponding residual was smaller than the iLQR residual . Thus, Figure 3 . 3 A and Figure 3 . 3 B show a perfect fit in the states and the control, respectively. For a fair comparison between iLQR and iMPC, we modified the cost function of iLQR gradient - descent algorithm to be like that of iMPC, that is 35 , ( 3 . 16 ) where the gradient was modified accordingly, and denotes the modified gradient - descent algorithm. The result was ( 3 . 17 ) The state and control of the iLQR are shown in Figure 3 . 4 A and Figure 3 . 4 B, respectively. With a static full - state feedback gain, the state fitting was good, however, the control action magnitude exceeded the bound of as shown in Figure 3 . 2 . Thus, the full - state feedback gain is not appropriate for the human motor control analysis. It may produce a good fit in the state trajectories, but violates the physiological constraints. Applying the gradient descent iLQR with the mo dified cost function resulted in an LQR gain that was far from . Therefore, the state fitting ( Figure 3 . 4 A) was worse than that in Figure 3 . 2 A. This shows that even though a good fit may be obta ined with a full - state feedback gain, the gradient descent iLQR approach may not be reliable for fitting the output data. Moreover, the corresponding LQR control magnitude exceeded the bound of . This finding indicates that whether the feedback gain is an LQR one or an arbitrary one, the control may violate the physiological constraints in trunk motor control. In this simulation example, iMPC minimized the condition number, recovered the true weights, and had a perfect fit to the states. Moreover, the co ntrol action did not violate the bounds. This last finding is the main motivation for using iMPC to infer the control intent. 36 (A) (B) Figure 3 . 3 : (A) The iMPC states ( ) versus the reference states ( ). (B) The iMPC control ( ) versus the reference control ( ) . 37 (A) (B) Figure 3 . 4 : (A) The iLQR states ( ) versus the reference states ( ). ( B) The iLQR control ( ) versus the reference control ( ) . The iLQR here is the modified gradient - descent algorithm in ( 3. 16 ) . 38 3.5 The Seated Balance Experiment 3.5.1 The pHRI Our pHRI consisted a compliant robotic platform and a sitting subjec t who performed a seated balance test ( Figure 3 . 5 A). The robot consisted of a seat that was coupled with an electric motor (C062C - 13 - 3305, Kollmorgen, Radford, VA) to apply torques to the seat about the seat pivot . The subjec t had arms crossed on chest, knees at 90 degrees, and hip strapped to the robot. The ( ) about ( ) that was assumed to be about the fourth lumbar vertebra ( ). During the experiment, both (1 ) a pseudorandom ternary torque perturbation , and (2) a stiffness - and - damping torque were applied in the coronal plane about the pivo t . The subject was instructed to balance on the robot using movements in the coronal plane. was measured by an encoder embedded in the motor, while was measured using string potentiometer s (SP2 - 50, Celesco, Chatsworth, CA) attached to a ch est harness and a belt at the fourth thoracic ( ) and levels, respectively. Data was logged using a DAQ (PCI - DAS6036, Measurement Computing, Norton, MA) and a quadrature encoder board (PCI - QUAD04, Measurement Computing, Norton, MA). Our pHRI was show n to be reliable (i.e., data is reproducible) [26] . The sampling rate was set to 1000 Hz , but data was sampled at 100 Hz for the analysis in this study . The pHRI is a two - degree - of - freedom system shown in Figure 3 . 5 B. The mathematical model is then given by 39 ( 3 . 18 ) where , and subject applies a torque abo ut , and has intrinsic lumbar stiffness and damping coefficients and , respectively. The lower body, below , was lumped with the robotic seat into a single mass and moment of inertia about the center of mass (COM1) [33] . The distance between COM1 and the pivot is . The upper body, above , was lumped into a mass and moment of inertia about its own center of mass (COM2), which is at a distance from . itself is at a distance from the pivot . The dynamical model was formulated as a closed - loop feedback control system by adding bi ological control, sensing, and muscle actuation. First, the dynamic property of the muscle actuation was assumed to be a first - order filter [31] of the form , where is the muscle time constant. Also, a sensory delay was introduced in the feedback loop [24], [83] and it was approximated by a fifth order Padé approximation [31] . The block diagram of the closed - loop system is shown in Figure 3 . 6 . The pHRI block is the linearized model of (6). The plant, therefore, has ten states: pHRI states , Padé approximation states , and the first - order filter state . 40 ( A) (B) Figure 3 . 5 : (A) The pHRI for our seated - balance experiment in the upright equilibrium position . (B) The pHRI model of our seated - bal ance experiment. The variables with subscript are for the robot model and the those with subscript are for the human model. 41 Safety is essential for a successful pHRI [84] . For safety purposes in our pHRI, limit switches were used to deactivate the perturba tion torque when the robotic seat exceeded an angle of . Both the operator and the subject held emergency switches to deactivate the perturbation torque when pressed. To ensure passive safety [84] , the impedance control (spring - damper action) was never deactivated if the subject was on the robot. The values of stiffness and damping as well as the perturbation torque amplitude were normalized at the beginning of the experiment to make the task challenging, but achievable for the s ubject. We defined a challenging and achievable task as the one in which the seat reaches angles close to . 3.5.2 Subject We used the collected data of one subject who was a 56 - year - old male with 79.7 kg body mass and 166 cm height. The experimental lower - body and upper - body angles of one trial are shown in Figure 3 . 7 . The parameters and were taken from [33] , and and were obtained from the regression relations in [85] . The model was discretized with the sampling time . Figure 3 . 6 : Block diagram of the closed - loop system of trunk motor control during seated balan cing . The experimental output measurements are the time series of . 42 Figure 3 . 7 : MPC trajectories ( dashed curves ) versus experimental trajectories ( continuous curves ) . 3.5.3 Control Intent Priess et al. [31] estimated the parameters defined as ( 3 . 19 ) such that the controller was assumed to be a linear state feedback gain with the state vector of ( 3 . 20 ) The results with only four states can be extended to our case of 10 states if we forced the structure This structure is co nsistent with the literature 43 [24], [83] that considered the reflexive feedback to have only position and velocity components. In our case, and are the position gains, and and are the velocity gains. All parameters are given in Table 3 . 2 . The iLQR problem in the LMI formulation was infeasible. Therefore, the gradient descent iLQR problem wa s used. Next, the auxiliary optimization problem was solved and the closest LQR gain to the estimated gain was . ( 3 . 21 ) The MPC pro blem was formulated with bounds on the states and inputs given by ( 3 . 22 ) The bounds on the pHRI states (positions and velocities) were derived from the experimental data, such that the bounds included the maximum measured output values. The other states were not bounded since there was no information or measurement to suggest an estimate of the ir bounds. The maximum torque was set to , which is almost half of the near maximum lateral bending torque for male subjects [82] . The horizon used was . 44 Table 3 . 2 : Physiological parameters of the model in Figure 3 . 6 . They are divided in two groups separated by a double line. The upper group include fixed parameters. The lowe r one includes estimated parameters to fit the experimental data. Parameter(unit) Value Description mass of both the lower body and the seat mass of the upper body lumbar stiffness lumbar damping robot stiffness robot damping control gain of control gain of control gain of control gain of inertia of both the lower body and the seat about COM1 inertia of the upper body about COM2 dista nce, see Figure 3 . 5 B dista nce, see Figure 3 . 5 B distance, see Figure 3 . 5 B sensory delay muscle time - constant The iMPC prob lem was solved. The fitting is shown in Figure 3 . 7 . The resulting matrix was symmetric and dense, so it could not give a clear picture of what linear combinations of the states were most penalized. Therefore, we transfor med into a diagonal matrix. Let us consider a similarity transformation given by where is a transformed state vector whose elements are linear combinations of the original states . If is selected to be an orthogonal matrix whose columns are the eigenvectors of , and given the constraint then is a diagonal matrix whose diagonal elements are the 45 eigenvalues of . The result is ( 3 . 23 ) The goodness of fit was quantified using variance accounted for ( ) [25] give n by ( 3 . 24 ) where is the number of observations. As approaches , the fitting is b etter. For our inferred control intent, the goodness of fit was ( 3 . 25 ) 3.6 Discussion The largest eigenvalue of was associated with the eigenvector given in t he last column of , that is ( 3 . 26 ) where the four states are defined in ( 3. 20 ). Ther e fore, the subject penalized the weighted 46 summation of the bod y angles and velocities such that the angles were penalized more than the velocities. Thus, the control intent was dominated by minimizing . This reflects the original task instruction to maintain balance on the robot, i.e., kee p all angles near zero. According to the positive angles convention in Figure 3 . 5 B if the robot moves in a positive direction then it is intuitive that the subject needs to move the upper body in a negative direction to stabiliz e on the robot. This intuitive action agrees with our findings because both and have the same sign in the dominant cost . In other words, this cost can be minimized only if has the opposite sign of . T he delay states, , appeared in the eigenvector corresponding to the largest eigenvalue, but their weights were relatively smaller than the pHRI states, . Moreover, the delay states had dominant weights in the eigenvector c orresponding to the second largest eigenvalue. This indicates that the delay states were second in priority in the control intent. The control penalty, , was forced to in the iMPC problem. By comparing the values of and , we find that the subject penalized the states much more than control effort . This finding is consistent with the literature [71] . A similar finding was reported in [56] , but with a different cost fun ction. The conditions required for designing an LQR controller were satisfied by the iMPC solution. That is the triplet was reachable and observable [86] , where and are the state - space matrices of t he plant shown in Figure 3 . 6 . This finding is consistent with the MPC stability conditions, where an LQR problem is solved using . 47 3.7 Conclusion We presented an iMPC algorithm to solve iMPC problems with output me asurements only. The algorithm utilized solving the iLQR problem [71] . We introduced an auxiliary optimization problem that updated the LQR weights to a good starting point for the MPC weights. Then, we developed a two - level optimizati on scheme to solve the iMPC problem. The algorithm was tested by a simulation example in which the control action saturated. Our algorithm could recover the true MPC weights. The iMPC algorithm was then applied to the experimental data collected from our p HRI to infer the control intent of one subject. We advanced the state - of - the - art work of iLQR [71] in two points. First, the plant included the pHRI along with the sensory delay and the muscle dynamics, which were not included in [71] . This gives more insight into the relative significance of the states during the task. Second, the iMPC approach is superior to the iLQR one because it models the physiological constraints of the motor control. We showed that the recov ered weight matrix reasonably represented the instructions given to the subject and yielded an accepted goodness of fit. 48 CHAPTER 4 SELECTING SENSITIVE PARAMETER SUBSETS IN DYNAMICAL MODELS WIT H APPLICATION TO BIOMECHANICAL SYSTEM IDENTIFICATION T his chapter is pub lished at the ASME Journal of Biomechanical Engineering [87] . Abstract Estimating many parameters of biomechanical systems with limited data may achieve good fit, but may also increase 95% confidence intervals in parameter estimates. This results in poor identifiability in the estimation problem. Therefore, we propose a novel method to select sensitive biomechanical model parameters that should be estimat ed, while fixing the remaining parameters to values obtained from preliminary estimation. Our method relies on identifying the parameters to which the measurement output is most sensitive. The proposed method is based on the Fisher Information Matrix (FIM) . It was compared against the nonlinear Least Absolute Shrinkage and Selection Operator (LASSO) method to guide modelers on the pros and cons of our FIM method. We present an application identifying a biomechanical parametric model of a head position - track ing task for ten human subjects. Using measured data, our method (1) reduced model complexity by only requiring five out of twelve parameters to be estimated, (2) significantly reduced parameter 95% confidence intervals by up to 89% of the original confide nce interval, (3) maintained goodness of fit measured by Variance Accounted For (VAF) at 82%, (4) reduced computation time, where our FIM method was 164 times faster than the LASSO method, and (5) 49 selected similar sensitive parameters to the LASSO method, where three out of five selected sensitive parameters were shared by FIM and LASSO methods. 4.1 Introduction The estimated parameters in models of sensorimotor control provide insight into the relative contributions of various feedback pathways for control of several body segments such as head - neck [22], [23], [88], [89] , trunk [20], [24], [90] , and arm [91] [93 ] models. However, some difficulties are associated with these models, which limit their application. Difficulties occur mainly because these models can produce a large number of parameters that need to be estimated [25] . At the sa me time, the information in experimental data is limited because subjects do not tolerate a long data collection process before fatiguing, and non - invasive external sensors are not available to measure all signals of interest. Overfitting is a consequence of estimating too many parameters in the model with limited measurements. Overfitting is well recognized in polynomial fitting, e.g., fitting a 14 th order polynomial to 15 data points yields a perfect fit but no useful information about the relevance of fi tted model parameters to the behavior tested. Similarly, as parametric models become more complex with more parameters to be estimated, goodness of fit may improve. However, the estimated parameters may have more variability [24], [31] because multiple parameter solutions yield very similar model responses [88] . This difficulty leads to poor identifiability in the estimation problem [52], [94] . Moreover, with increased model complexity, the estimation is more likely to be misled by noise possibly resulting in incorrect parameter estimates [95] . Goodness of fit must be balanced against the number of model parameters estimated to avoid overfitting with too many parameters for the in formation in data. 50 A good trade - off between number of parameters estimated and goodness of fit can be achieved by fixing some parameters to values obtained from preliminary estimation while estimating others [24], [31] . We propose only those sensitive parameters significantly affecting the model response be estimated. The selection of parameters sensitive to the data allows those parameters to be accurately estimated by that data. Sensitive parameters have been discussed in the literature in the setting of optimal experiment design, e.g., experiments were designed such that the parameters would be identifiable from the experimental data [31], [96], [97] . In this study, we consider an inverse view where sensitive parameters are identified from the experimental data to provide an optimal trade - off between number of parameters estimated and goodn ess of fit. Our method uses Fisher Information Matrix (FIM) to select sensitive model parameters using available experimental data in a novel problem formulation for biomechanical modeling. Other methods of achieving this trade - off have been reported in th e literature. An ad hoc method of selecting a model out of a few possible combinations was proposed in [24] , but the combinations [57], [98] and Information COMPlexity (ICOMP) [99] are very expensive computationally because they require every combination of the model parameters being estimated to determine their goodness of fit and there are 4083 combinations of two parameters or more chosen from a 12 - parameter model in our [100] also requires multiple estimations to update a tuning parameter. The number of estimations in LASSO is, in general, smaller than the number of parameter combinations. Other machine learning methods provide bla ck - box models that have no physiological interpretations, e.g., Gaussian process regression in gait kinematics [101] . 51 The goal of this study is to introduce and evaluate a novel FIM method to select sensitive parameters in m odeling sensorimotor control, in this case the head - neck segment. Our method improves on previously published methods because (1) it is faster as it relies on a single preliminary estimation, and (2) the full parametric model has a physiological interpreta tion. Early work of our method was presented without evaluating its performance [102] . In this study, we demonstrate s ensitive parameter selection is robust to the method used by comparing our FIM method to the LASSO method [100] . Moreover, we compare the two methods in terms of goodness of fit, parameter confidence intervals , and computation time. This information will guide modelers on the pros and cons of this method. 4.2 Methods 4.2.1 Subjects Ten healthy subjects part icipated in the experimental study ( Table 4 . 1 ). They did not have any history of neck pain lasting more than three days or any neurological motor control impairment. Health Institutional Review Board approved the test protocol. The subjec ts signed an informed consent before participating. Table 4 . 1 : standard de viation. Females Males No. of subjects Height (m) Weight (kg) Age (years) 52 4.2.2 Data Collection A head position tracking task in the horizontal, yaw, plane ( Figure 4 . 1 A) provided data on subject head response to vision - based commands. The experiment was statistically shown to be reliable [19] . Each subject wore a helmet with two attac hed string potentiometers (SP2 - 50, Celesco, Chatsworth, CA) to measure the head rotation. On a monitor (SyncMaster SA650, Samsung; height 27 cm, width 47.5 cm) located 1 meter away from the subject, two markers were displayed: reference command signal and the measured head position response . Subjects were instructed to rotate their heads in the horizontal plane to track the reference signal. The input signal was a pseudorandom sequence of steps with random step duration and amplitud e, such that the amplitudes were bounded in the range Figure 4 . 2 ). Each subject performed three 30 - second trials. Data was collected using a data acquisition system (cDAQ - 9172, National Instruments, Austin, TX ) and sampled at 60 Hz. 4.2.3 Parametric Model The block diagram of head position tracking model ( Figure 4 . 1 B) adopted published physiological and feedback control models [22], [23], [89] , where the output was modified to be angular position instead of angular acceleration. This change in output from acceleration to position does not change the behavior modelled nor the definitions of the physiological parameters defined in those previous investig ations. The model parameters (physiological descriptions given in Table 4 . 2 ) are divided into two groups. One group has the candidate parameters, out of which we would select the sensitive parameters to be estimat ed. They are ( 4 . 1 ) 53 Table 4 . 2 : The neurophysiological pa rameters of the head position tracking task along with their physiological description . The lower and upper bounds were obtained from preliminary (section 4.3 ) . Parameters Max M in Description visual feedback gain [22] vestibular feedback gain [23] proprioceptive feedback gain [23] visual feedback delay [22] lead time constant of the irregular vestibular afferent neurons [23] lead time constant of the central nervous system [23] lag time constant of the irregular vestibular afferent neurons [23] lag time constant of the central nervous system [23] first lead time constant of the neck muscle spindle [23] second lead time constant of the neck muscle spindle [23] intrinsic damping [22], [23] intrinsic stiffness [22], [23] head inertia [23] torque converter time constant [23] 54 (A) (B) Figure 4 . 1 : (A) Experimental setup for head position tracking. The reference marker is r ( t ) and the head position marker is y ( t ) . (B) The block diagram of the head position tracking model adopted from [22], [89] . The transfer functions come from [22], [23] after accounting for the change in output from acceleration to position. The second group, , was fixed to values found in the literature [23] . The bounds in Table 4 . 2 were obtained from preliminary parameter estimates discussed in Section 4.3 below . 55 Figure 4 . 2 : The best - fit model output (during a trial) when estimating the FIM - and LASSO - response. 4.2.4 Parameter Estimation Parameter estimates define physiological models from measured task response. Parameter normalization is required because physiologi cal parameters often have quite different magnitudes, which cause numerical problems in compu ter optimization. Each normalized parameter for a model parameter is given by ( 4 . 2 ) for the vector of parameters , is an offset to make the true bounds symmetric around 0, an d is a gain that normalizes the symmetric parameter bounds to . Estimates were computed 56 by minimizing least squared error, i.e., ( 4 . 3 ) where is the error at data sample , is the measured head position output, is the simulated model output computed at sample and normalized parameter , and is the number of observation samples. Goodness of fit of the simulated model data to the experimental data was measured by Variance Accounted For ( ) [24] ( 4 . 4 ) As approaches , the model response approaches the measured response and goodness of fit is better. Pa rameter 95% Confidence Interval (CI) estimates are critical for evaluating the success of any method. We used model - based bootstrap [103] ( Figure 4 . 3 ). We generated 1000 bootstrap replicates of the experimental data ( ) per subject. Then, we estimated the parameters for every replicate to get ( ). The 95% CI was computed from 2.5 to 97.5 percentiles of the empirical distribution of . Paired t - tests 05 were used to statistically compare the results between the methods. 57 Figure 4 . 3 : An outline of model - based bootstrap [103] using 1000 replicates to estimate 95% confidence interval (CI) of the parameter vec tor . 4.2.5 Sensitive Parameter Selection Sensitive Parameters were selected using FIM ( Appendix A ) that is a real matrix denoted by , where is the number of model parameters. FIM is computed at the nominal values of the model parameters , which were obtained via preliminary estimation. The diagonal elements of ( Appendix C ) scans all possible parameter combinations to ensure that the estimation error variance ( ) is below a user - defined threshold , and maximize the number of sensitive 58 parameters ( ) to improve goodness of fit. The user must provide FIM ( Appendix A ), allowable normalized variance , and the desired minimum number of sensitive parameters . LASSO ( Appendix B ) is an alternative method that shrinks insensitive normalized parameters to zero while estimating the remaining sensitive parameters. This is done by adding a new term (regularization) to the cost function in ( 4. 3). In our application, shrinking the normalized parameters to zero corres ponds to shrinking the actual parameters to the values centered between their bounds, which is likely not as correct as typical physiological values. Therefore, a modification of the regularization in LASSO was made ( Appendix C ) to ensure that parameters were pushed toward their typical values [100] . 4.3 R esults Nominal parameter values are required to compute FIM and to be the LASSO typical parameters. Adopting values from prior literature yielded a poor fit, which might be because the literature either used a simpler model [22] or di d not have a visual tracking task [23] . Therefore, w e did preliminary parameter estimation to obtain the nominal parameter values over all three tr ia ls ( ) for each subject. The cost function in ( 4. 3) becomes ( 4 . 5 ) where is the trial index, and is the error of the trial. M inimizing this cost function yields the best - fit parameters of a subject for the three trials collectively ( Table 4 . 3 ). In the FIM method, w - fit parameters to estimate paramete variance for each subject. From each FIM, a subset or more of subject - specific sensitive parameters were selected. Then, only the common parameter subsets among all subjects were selected. This 59 is a more accurate process than building one FIM from the selecting sensitive parameters based on this averaged FIM. We applied the FIM method to the experimental data with ten equally - spaced values of from 0.05 to 0.95 and ( Figure 4 . 4 ). The number of selected sensitive parameters was between 3 and 5. We selected at because this solution included a single sub set . The selected sub set was . Figure 4 . 4 : Results of applying the FIM method ( Appendix C ) w ith to the experimental data . The optimal number of sensitive parameters per subset is . The number of subsets (each with parameters) is . 60 Table 4 . 3 : Parameter esti mates of the complex model with 12 estimated parameters (Nominal) and the reduced models with five sensitive parameters selected by FIM and LASSO methods ( and , respectively). Parameter 95% confidence interval width was estimated for the complex and reduced models . Goodness of fit was measured by Variance Accounted For ( ). All values are given as mean standard deviation across subjects. Parameters Parameter Estimates 95% Confidence Interval Width (Nominal) 428 ±125 279 §, ±27 401 ±163 546 ±128 62 ±19 236 # ±122 4860 ±2205 - 8055 § ±1371 7274 ±806 - 3623 # ±1418 8 7.9 ±55.3 16.2 §, ±11.9 36.4 § ±30.5 123.0 ±49.8 18.2 ±6.0 68.7 # ±36.9 0.263 ±0.027 0.242 ±0.024 0.226 § ±0.027 0.034 ±0.010 0.020 ±0.010 0.033 ±0.010 0.099 ±0.063 0.069 ±0.013 - 0.134 ± 0.036 0.030 # ± 0.015 - 0.675 ±0.280 0.918 § ±0.077 - 0.544 ± 0.145 0.138 # ± 0.069 - 1.733 ±1.118 - 2.200 ±1.256 3.086 ±0.495 - 1.736 # ±0.534 27.2 ±16.2 - - 48.7 ± 2.4 - - 0.433 ±0.381 - - 0.932 ± 0.075 - - 0.439 ±0.325 - - 0.877 ± 0.086 - - 61 Table 4 . 3 1.690 ±1.213 - - 4.659 ± 0.144 - - 2.170 ±1.736 - - 4.758 ± 0.114 - - 82.26 ±7.67 82.65 ±7.57 81.87 ±7.56 - - - § The parameter estimate in the reduced model (with five estimated parameters selected by FIM and LASSO parameters nominal values). The parameter estimate in the FIM - the LASSO - selected parameters. # The parameter 95% confidence interval in the reduced model (from FIM or LASSO methods) was ex model. that from the LASSO method. In the LASSO method, we considered the mean of nominal parameter estimates as the LASSO typical parameters. If we us ed subject - specific parameters as the typical values, all parameters would be pushed to these typical values since they are the best - fit ones for that subject. After selecting a sensitive parameter subset for each subject, we selected the most frequent sub set among all subjects. We applied the LASSO method to the experimental data where was set to 5 to obtain comparable results. The selected subset was . shared 3 parameters (out of 5) with . The 95% CI of the three common parameters were r educed substantially more by FIM than LASSO. The CI was reduced by LASSO to 43% and by FIM to 11% of the original CI ( Figure 4 . 5 A ). The CI was reduced by LASSO to 56% and by FIM to 15% of the original CI ( Figure 4 . 5 B ). The CI was reduced by LASSO to 96% and by FIM to 59% of the original 62 CI ( Figure 4 . 5 C ). The CI of the selected sensitive parameters ( and ) were significantly narrower tha n the original CI of all parameters ( Figure 4 . 5 and Table 4 . 3 ), except the versus comparison (p=0.69). Moreover, the CI of the common parameters in were significantly narrower than those CI in ( Table 4 . 3 ). G oodness of fit measure was determined by estimating the selected sensitive parameters ( and ) wh ile fixing the remaining parameters to the mean of the nominal values ( Table 4 . 3 and Figure 4 . 2 ). FIM and LASSO ( and , respectively) were almost equal, where and across subjects for 5 selected sensitive parameters. When estimating all 12 model parameters, goodness of fit was . Statistical analysis did not show significant differences between or versus (p=0.23 and 0.33, respectively). Very different computation times ( ) of FIM and LASSO methods were measured. All Xeon E5 - 2699 v3. The average computation times to build FIM and solve the FIM method was , whi le the average LASSO computation time was . Significant differences in parameter estimates occurred between the 12 estimated parameters for the complex model and each of the five estimated parameters selected by FIM and LASSO methods for reduced models. Three out of five parameters were significantly different from their nominal values: in , and in were different from their counterparts in ( Table 4 . 3 ). Moreover, the three parameters selected by both FIM and LASSO methods had significantly different estimates. 63 (A) (B) Figure 4 . 5 : Comparing the 95% con fidence interval (CI) between estimating ALL parameters, FIM - selected parameters ( ), and LASSO - selected parameters ( - estimate and the error bars are the mean CI (across subjects) for (A) , (B) , and ( C) that were selected as sensitive parameters by both FIM and LASSO methods. 64 Figure 4 . 5 (C) 4.4 Discussion In this study, we compared our novel FIM method to the LASSO method for selecting sensitive model parameters b ased on four performance measures: (1) number of sensitive parameters that are in - common between the FIM and LASSO methods, (2) parameter CI, (3) , and (4) computation time. First , the selected sensitive parameters, using FIM and LASSO methods, were qui te similar. Three out of five selected parameters were identical . Thus, sensitive parameter selection is considered robust to the method used. Second , the 95% CI of the three common parameters were reduced significantly by both FIM and LASSO methods 65 versus the original CI, e xcept the case ( Figure 4 . 5 ). This indicates the success of both methods in selecting sensitive model parameters. Moreover, the CI were reduced significantly more by FIM than LASSO ( Table 4 . 3 ), which demonstrates that our novel FIM method is superior to the LASSO method in reducing parameter 95% CI. Third , goodness of fit measure for FIM or LASSO methods did not show significant differences from that of estimating all 1 2 model parameters. This indicates that the two methods maintained goodness of fit although they reduced the number of estimated parameters. Fourth , the average computation time of FIM method was less than 1% of that in LASSO method. Our novel FIM method i s significantly faster than the LASSO method. Significant differences in parameter estimates exist between the complex model with 12 estimated parameters ( ) and the reduced models with five estimated parameters ( and ). The reason might be that estimating all 12 parameters resulted in reaching a rather randomly found local minimum out of multiple local minima, due to poor identifiability, while estimating only sensitive parameters ensured the unique solution with good identifiability. Moreover, significant differences exist between estimates of the common parameters in and . These differences are a result of selecting dif ferent sensitive parameters whose estimates might be biased. Simulating a hypothetical subject, whose physiological parameters are known, showed that there was a bias in the estimates of the selected sensitive parameters. Once a sensitive parameter selecti and after a therapy), the bias effect will be relatively constant throughout the study. Then, our approach can provide otherwise impossible comparison results especial ly for a parametric model with poor identifiability. Future research around parametric estimation of sensorimotor control should investigate these issues. 66 4.5 Conclusion Sensitive parameter selection was robust to the method used. Both FIM and LASSO methods se lected quite similar sensitive parameters (three out of five selected parameters were identical). Both FIM and LASSO methods reduced the parameter 95% confidence intervals significantly up to 11% and 43% of the original CI, respectively. The FIM method red uced CI significantly more than LASSO. Both methods maintained the original goodness of fit. Across FIM and LASSO methods, goodness of fit measure was essentially equal (82%, Figure 4 . 2 ). The FIM method is 164 times faster than t he LASSO method. This information will guide modelers on the pros and cons of each method and indicates that FIM more efficiently reaches an appropriate reduced model parameter set. 67 CHAPTER 5 FEASIBILITY OF INCOR PORATING TEST - RETEST RELIABILITY AND MODE L DIVERSITY IN IDENTIFICATION OF KEY NEUROMUSCULAR PATHWAYS DURING HEA D POSITION TRACKING Abstract To study complex neuromuscular control pathways in human movement, biomechanical parametric models and system identification methods are employed. Although test - retest reliability is widely used to validate outcomes of motor control tasks, it was not incorporated in system identification methods. This study investigates the feasibility of incorporating test - retest reliability in our previously published method of selecti ng sensitive parameters. We consider the selected parameters via this novel approach to be the key neuromuscular parameters because they meet three criteria: reduced variability, improved goodness of fit, and excellent reliability. These criteria ensure th at parameter variability is below a user - defined value, the number of these parameters is maximized to enhance goodness of fit, and their test - retest reliability is above a user - defined value. We measured variability, goodness of fit, and reliability using Fisher information matrix, variance accounted for, and intraclass correlation, respectively. We also incorporated model diversity as a fourth optional criterion to narrow down the solution space of key parameters. We applied this approach to head position tracking in axial rotation and flexion/extension. Forty healthy subjects 68 respectively, we selected three key parameters out of twelve. The key parameters were associated with at least two neuromuscular pathways out of four modeled pathways (visual, proprioceptive, vestibular, and intrinsic), which is a measure of model diversity. The average goodness of fit was >69%. Therefore, it is feasible to incorporate rel iability and diversity in system identification of key neuromuscular pathways in our application . 5.1 Introduction Identification techniques for parametric models are needed to investigate the neuromuscular control pathways involved in human movement. U nfortun ately, we do not have direct access to individual pathways , but instead rely on gross measure s such as movement kinematics and electromyography [104] . Th erefore, neuromuscular p arametric models have been used to estimate the relative contribution of various neuromuscular control pathways [20], [24], [46], [105] . These parametric models typically have multiple ind ependent neuromuscular control signal pathways. Because the neurophysiology is complex with many pathways, these models also tend to be complex with a large number of parameters to be estimated [23], [24], [106] . For a fixed amount of information (measurement) , increased parameterization results i n increas ed parameter variability [24], [31], [107], [108] , which subsequently decrease s reliability. Reliability is usually assessed for nonparametric measures such as root mean square error and mean velocity [19], [26], [109], [110] . Reliability is rarel y assessed for estimated model parameters (we are aware of a single study that reported reliability of parameter estimates [111] ). No study, to our knowledge, considered incorporating reliability in a system identification method for estimating model parameters. Excellent parameter reliability is important for system identification method s to be useful in biomechanical and clinical applications [51], [112], [113] . Therefore, 69 reliability should be incorporated in existing system identificatio n methods for estimating model parameters. We previously developed a system identification method that estimate s only sensitive model parameters while fixing the remaining ones to preliminary values [87] . It allows selection across all possible parameter subsets, which makes it exhaustive . Because no parameter is eliminated (parameters are either estimated or fixed), all neuromuscular pathways are present in the final model. Sensitive parameters were selected based on parameter variability (measured by Fisher information matrix (FIM)) and goodness of fit. The main purpose of this study was to investigate the feasibility of incorporating test - retest paramet er reliability in our existing FIM method of selecting sensitive parameters . W e used a neuromuscular model of the head - neck system for position - tracking in axial rotation (Yaw) and flexion/extension (Pitch) . We hypothesize d that it is feasible to select se nsitive parameter subsets with excellent reliability . We consider these parameters to be the key model parameters substantially affecting head - neck control . 5.2 M ethods 5.2.1 Subjects Forty healthy subjects were recruited for this study ( Table 5 . 1 ) . They did not have neck pain lasting more than three days in the last year nor neurological or balance disorder s . The Michigan State All subjects sign ed an informed consent before participation. 70 5.2.2 Data Collection The experimental setup of head position tracking ( Fig ure 5 . 1 ) used a helmet , secured on the , with two string potentiometers (SP2 - 50, Celesco, Chatsworth, CA) attached to measure head angle [19] . The reference marker along with the head position marker were displayed on a monitor (Sy ncMaster SA650 , Samsung; 27 cm x 47.5 cm) fixed 1 m eter away from the subject. Subjects were instructed to track the reference marker as fast and accurate as possible. The reference signal was a pseudorandom sequence with random step durations (0.3 - 0.9 sec onds) and amplitudes such that the maximum amplitude was centered around (the upright head position). Subjects completed each task ( Yaw and Pitch ) on two separate visits . Each task consisted of two 15 - second practice trials and three 30 - second trial s during each visit. Data was sampled at 50 Hz and collected with a data acquisition card ( cDAQ - 9172 , National Instruments, Austin, TX) . 5.2.3 Data Analysis The parametric model s of the task s ( Fig ure 5 . 2 ) used previously published phys iological subsystems [22], [23], [89] , but the output was modified from head acceleration to head position . This modification does not alter the model behavior, nor the physiological parameter interpretation given in those studies. The physiological parameters ( Table 5 . 2 ) were divided into parameters that may be estimated (candidate key parameters) , ( 5 . 1 ) and fixed parameters . They are the physiological parameters of visual, vestibular, proprioceptive, and intrinsic pathways, torque converter, and the head - neck plant ( Fig ure 5 . 2 ). 71 ( A ) ( B ) Fig ure 5 . 1 : The experimental setup of head position tracking tasks in ( A ) axial rotation (Yaw), and ( B ) flexion - extension (Pitch). Table 5 . 1 : Demographic characteristics of the subjects given as (mean ± standard deviation). Males (n=19) Females (n=21) Height (cm) Weight (kg) Age (years) 72 * the gravity effect ( ) was used for Pitch only. Fig ure 5 . 2 : The block diagram of neuromuscular control for the head position tracking tasks. The contributing pathways are visual, vestibular, proprioce ptive, and intrinsic. Table 5 . 2 : Physiological parameters used in the neck motor control model ( Fig ure 5 . 2 ) and their bounds. Bounds were based on preliminary estimates ( Table 5 . 3 ). has model parameters that may be estimated (candidate key parameters), while has fixed parameters. Parameter s Max Min Description visual gain [22] vestibular gain [23] proprioceptive gain [23] visual delay [22] lead time constant of the irregular vestibular afferent neurons [23] lead time constant of the central nervous system [23] 73 Table 5 . 2 lag time constant of the irregular vestibular afferent neurons [23] lag time constant of the central nervous system [23] first lead time constant of the neck muscle spindle [23] second lead time constant of the neck muscle spindle [23] intrinsic neck damping [23] intrinsic neck stiffness [23] head - neck inertia in Yaw / Pitch about the center of mass torque converter time constant head - neck mass distance from C6 - C7 to head - neck center of mass Estimated model parameters identify the physiological contributions of neuromuscular control pathways . Because the parameters had different magnitudes ( Table 5 . 2 ), which may cause numerical issues in optimizatio n , the parameters were normalized . For the parameter , its normalized counterpart is . Any subset is estimated in the same way as , therefore, we use for generalization. Parameters were estimated by nonlinear least squa res , ( 5 . 2 ) where is the concatenated head position measurement of three trials, is the simulated model 74 output , and is the 2 - norm of a vector. W e used lsqnonlin 5. 2) . Preliminary estimat es of were obtained per subject per visit (nominal values, Table 5 . 3 and Table 5 . 4 unds ( Table 5 . 2 ) . For a given subset , its parameters were estimated (per subject per visit) while fixing the remaining par ameters to the mean of the nominal values ( preliminary estimates ) . Goodness of fit was measured by V aria nce A ccounted F or ( ) [24], [92] ( 5 . 3 ) As approaches 100%, the simulated data better fits the experimental data. Sensitive parameter subsets were selected using the FIM method [87] . FIM is a square matrix computed at the normalized true model parameters and denoted by true parameters are unknown, we used the nominal values ( Table 5 . 3 and Table 5 . 4 ). The diagonal elements of FIM inverse ( ) are parameter est imation error variance ( represent parameter variability). The FIM method ( 5 . 4 ) searches all possible parameter subsets to maximize the number of sensitive parameters per subset ( ) while bounding parameter variability below a user - defined value [8 7] . Maximizing the number of estimated parameters indirectly improves goodness of fit . The first constraint on the reciprocal condition number to be larger than the machine epsilon guarantees the inverse matrix operation ( ) is well - conditi oned. returns the maximum diagonal element , and is the maximum allowable parameter variability (a user - defined value) . This constraint 75 eliminates any parameter subset that has a parameter variability larger than . To be conservative, this constraint was modulated to bound across all subjects and visits below . That is, , ( 5 . 5 ) where is the subject index , is the visit index , and is the normalized nominal parameters of the subject in the visit . Sensitive parameter subsets were selected at 10 equally - spaced points of from 0.05 to 0.95. B etween - visit reliability is important in our application. Between - visit comparison of parameters can detect the effect of an intervention introduced between two visits, which is a long - term goal of our work. Therefore, b etwe en - visit reliability of the normalized sensitive parameter estimates was computed using intraclass correlation coefficient ICC(A, k) with 95% confidence interval (CI) [49] . The reliability of a normalized parameter subset was considered the smallest ICC of its parameters. That is , ( 5 . 6 ) where is the largest number of sensitive parameters in this subset. The ICCs were classified as poor (<0.40), fair (0.40 - 0.59), good (0.60 - 0.74), or excellent (0.75 - 1.00) [50] . To investigate the feasibility of incorporating reliability in the FIM method, we used three criteria to select key model parameters for estimation. They are parameter variability, goodness of fit, and parameter reliability. The first two already exist in the FIM method as discussed earlier. The reliability criterion is , ( 5 . 7 ) 76 where is the minimum allowable parameter reliability (a user - defined value). If there is more than one feasible solution (two or m ore parameter subsets) after applying the three criteria, a fourth optional criterion reflecting model diversity was applied to reduce the solution space. Table 5 . 3 : The estimates of preliminary (Nominal) and selected key model parameters and the corresponding Variance Accounted for ( ) in Yaw. Values are given as mean ± standard deviation across subjects (and visits for Nominal). Between - visit reliability was measured by ICC(A, k) with 95% confidence interv al (CI) . Yaw Parameter s Nominal Selected Visit 1 Visit 2 ICC, 95% CI 77 Table 5 . 4 : The estimates of preliminary (Nominal) and selected key model parameters and the corresponding Variance Accounted for ( ) in Pitch. Values are given as mean ± standard deviation across subjects (and visits for Nominal). Between - visit reliability was measured by ICC(A, k) with 95% confidence interval (CI) . Pitch Parameter s Nominal Selected Visit 1 Visit 2 ICC, 95% CI 78 ( A ) ( B ) Fig ure 5 . 3 : The solution space (the sensitive parameter subsets) of the FIM method for ( B ) Yaw, and ( B ) Pitch. The solution space was discretized at 10 values of (the up per bound of parameter variability). is the maximum number of sensitive parameters per subset. is the number of sensitive subsets with parameters and parameter variability . 79 Table 5 . 5 : The sensitive subsets (solution space) of the FIM method ( Fig ure 5 . 3 ). The parameters for estimation are hig hlighted. interval Yaw Pitch (0,0.05] (0.05,0.15] (0.15,0.25] (0.25,0.35] (0.35,0.45] (0.45,0.55] (0.55,0.65] (0.65,0.75] (0.75,0.85] 80 Table 5 . 5 § § Only the case of estimating all 12 model parameters is shown. Diversity in our model may be postulated as the presence of independent neuromuscular pathway signals in the model. Therefore, selecting parameters in as many pathways as possible yields a "diverse" parametric model. In our model ( Fig ure 5 . 2 ), we have four pathways: visual, vestibular, proprioceptive, and intrinsic. A diversity index reflects how many pathways ar e represented by the selected parameters. 5.3 Results Sensitive parameter subsets were selected using the FIM method. The maximum number of sensitive parameters per subset ranged from 2 to 3, and from 2 to 4 for Yaw and Pitch, respectively ( Fig ure 5 . 3 ). For the interval in the inequality (4), the number of sensitive subsets is ( Fig ure 5 . 3 ). The sensitive subsets ( Table 5 . 5 ) are listed f or the intermediate intervals such as (0.05,0.15]. 81 The ICCs of the sensitive parameters were higher (more reliable) than those of all 12 model parameters ( Table 5 . 5 ). The reliability of a sensitive subset was considered the small est ICC of its parameters. Therefore, the reliability of the sensitive subsets ranged from 0.07 to 0.85 and from 0.30 to 0.91 for Yaw and Pitch, respectively. Key model parameters for estimation were selected using the four criteria defined earlier . First, for the variability criterion in ( 5. 4), we selected , which eliminated all subsets below the 4 th row in Table 5 . 5 . Second, selecting the maximum number of sensitive parameters per subset eliminated the upper three rows in Yaw and the 1 st row in Pitc h. Third, for excellent reliability, we selected . This led to one feasible solution in Yaw (diversity index=2), and three feasible solutions in Pitch (diversity indices = 3, 2, and 2, respectively). Therefore, the selected key parameters for Yaw and Pitch were , and , respectively. The key parameters were estimated while fixing th e remaining ones to their nominal values ( Table 5 . 3 and Table 5 . 4 ) . In Yaw, the 95% CI of ICCs ranged from 0.55 to 0.9 3 . In Pitch, the 95% CI of ICCs ranged from 0.69 to 0.96. 5.4 Discussion We hypothesize d that it is feasible to incorporate parameter reliability in our FIM method to select key model parameters that substantially affect head - neck control . In our application, the reliability of the sensitive parameters ranged from very poor to excellent, howe ver, they were still more reliable than all model parameters ( Table 5 . 5 ). This indicates that the FIM method of selecting sensitive parameters improves reliability, but it does not guarantee excellent reliability. Hendershot 82 et al. found a similar wide range of parameter reliability (from 0.19 to 0.72) [111] ; however, they did not select sensitive parameters as in our case. Out of the sensitive parameter subset s ( Table 5 . 5 ), there was more than one subset with excellent reliability ( ). Therefore, it is feasible to incorporate parameter reliability in our FIM method. Reducing the solution space of the FIM method is an advantage of adding the third reliability criterion. For example, in Pitch, the solution space was reduced from 6 sensitive parameter subsets to 3. This helps modelers reach a better decision by eliminating less - reliable sensitive parameter subsets. To reduce the so lution space further, we proposed selecting subsets with the most diverse parameters, i.e., parameters from different feedback pathways (visual, vestibular, proprioceptive, and intrinsic). Using the variability, goodness of fit, reliability, and diversity criteria on data of forty subjects, we selected one subset of key model parameters in Yaw and Pitch. However, the modeler might have a specific preference that can replace this diversity criterion. For example, if a modeler expects major differences to be present in the vestibular pathway before and after an intervention, then selecting subsets with more parameters from the vestibular pathway is a better choice . Some physiological parameters and feedback pathways are relatively more important than others in head position tracking tasks. In Yaw , the visual, proprioceptive, vestibular, and intrinsic parameters were in 7, 6, 5, and 0 sensitive subsets out of 9, respectively ( Table 5 . 5 ). In Pitch , the same parameters were in 14, 9, 13, and 0 sensitive subsets out of 16, respectively. This suggests that the visual pathway is the most important feedback in head position tracking tasks. This finding validates the previously published simplified models that considered only the visual feedba ck for tracking tasks [22], [81] . The results also suggest that the vestibular pathway is more important for Pitch than for Yaw. It might be due to the effect of gravity in Pitch . By investigating the visual 83 parameters further, the visual delay was in 6 out of 9 and 14 out of 16 sensitive subsets in Yaw and Pitch , respectively. On the other hand, the visual gain was in 3 and 4 subsets in Yaw and Pitch , respectively. This may suggest that the visual delay is the most important parameter in the most important pathway contr ibuting to head position tracking performance in Yaw and Pitch. Some limitations of our study are worth noting. First, the physiological parameter bounds affect the selected sensitive parameter subsets because FIM is computed at the normalized nominal para meters. For example, if the bounds of a parameter are expanded, its normalized variability shrinks, i.e., this parameter becomes more sensitive than it really is. Therefore, reasonable parameter bounds must be used. Second, the key parameter estimates might be biased (offset from the true values), but they are more precise (with less variability) than the estimates of all model parameters [87] . This bias cannot be measured because the true subject parameters are unknown. However , for a given study, such as detecting the effect of an intervention introduced between two visits , using between - visit comparison of parameters, the bias effect is relatively constant, while the estimates are more precise. 84 APPENDICES 85 Appendix A Fisher Informat ion Matrix (FIM) Fisher information is a measure of the sensitivity of an observable random variable relative to some parameter . The experimental data value is more sensitive to the sensitive parameters. We used FIM to estimate sensitivity of each subset of the parameters as measured by estimation error variance. FIM was derived in [31] for the discrete - time system (given in state space) ( A . 1 ) where is the true parameters, is the measurement output , is a Gaussian measurement noise, and is the number of observations . The FIM computed at i s ( A . 2 ) where is the model output, and is the err or. The derived FIM formula is [31] ( A . 3 ) where is the input sequence, is the initial state, is the identity matrix, and is the Kronecker tensor product operator. 86 The estimation error converges to a normal distribution asymptotically under mild conditions [52], [114], [115] , that is ( A . 4 ) The FIM inverse is a lower bound of the covariance matrix of the estimated parameters if is an unbiased estimator, that is [52], [116] ( A . 5 ) where is known as Cramér - Rao Lower Bound (CRLB), returns the covariance matrix and the matrix ineq uality means that is positive semidefinite. This lower bound ( ) presents the achievable performance limit by any unbiased estimator. of the unnormalized parameters may be misleading because of the different parameter scales. Therefore, we use d FIM of the normalized parameters given by ( A . 6 ) 87 Appendix B L east Absolute Shrinkage and Selection Operator (LASSO) The LASSO algorithm was developed as a method to simultaneously estimate parameters and reduce model complexity by shrinking insensitive parameters to zero, thereby removing them from the model [117] . The LASSO algorithm is an ordinary least squares estimator with the addition of an norm penalization on the parameter vector [118] . The use of norm results in solutions where some parameters are shrunk exactly to zero [117] . This res ults in selecting parameters that are the most sensitive in the model characterization, whi le simultaneously providing estimates of these sensitive parameters. The LASSO minimization, in its most general form, is ( B . 1 ) where is a tuning parameter to trade - off between goodness of fit and the number of selected sensitive parameters, and is the norm. This method can be applied to our problem by tuning it should be noted that LASSO is a biased estimator. This is a result of the entire vector of param eters being compressed, resulting in compressed estimates of the parameters which remain in the model. 88 Appendix C FIM and LASSO Methods for Sensitive Parameter Selection The goal here is to find the pool of candidate subsets of sensitive parameters. Let and , then ( C . 1 ) where is a selection matri x. For example, to select the first and last parameters of , the selection matrix is ( C . 2 ) FIM Method T o compute FIM sub - matrices , where , we did not repeat the calculations in (A.6). Instead, we computed once, then obtained the sub - matrix ( C . 3 ) The FIM - based selection problem must fulfill three requirements. First, it should select sensitive parameters. Therefore, we bounded the estimation error variance from above. Second, there should be no numerical issue s with the matrix inversion. Therefore, we set a lower bound on the reciprocal condition number of [119] . Third , to consider goodness of fit, we decided to maximize the number of selected sensitive parameters. 89 Consequently, the FIM - based selection problem is ( C . 4 ) where returns the length of a vector, returns the reciprocal condition number, is the machine epsilon, and returns the maximum diagonal element of the input square matrix. and are us er - specified constants which are the minimum number of sensitive parameters and maximum allowable estimation error variance of any normalized parameter, respectively. Given and , the optimization is solved by investigat ing all possible parameter subsets ( ) to check their feasibility (satisfy the constraints). Then, the subset(s) with the largest number of parameters is selected. The solution is one or more subsets of sensitive parameters . Let cardinality of (number of different subsets) be . There are two important aspects to our solution of the problem shown in (C. 4 ) . First, the problem is formulated in a generic way to fit any application. Therefore, the user should provide the normal ized FIM , and . Second, the problem is solved by investigating all possible combinations of the parameters (combinatorial optimization). This method is feasible in our application since it has 12 parameters. If this is no t the case, the reader is referred to [120] for it s relaxation. 90 The choices of and need to be made carefully . Initially, we recommend setting as small as possible and trying different values of that are usually between 0 and 1, because the normalized pa rameters . As increases, the solution (selected parameter subsets) either (1) remains the same, (2) has more subsets, i.e. larger , while the optimal number of parameters per subset remains the same, or (3) h as a larger . If the optimization problem is infeasible, i.e. does not yield a solution, the user should decrease or increase . LASSO Method The LASSO method is ( C . 5 ) where is a vector containing the typical parameter values, is a tuning parameter to trade - off between goodness of fit and the number of selected sensitive parameters, a nd is the norm. As the estimated parameters get pushed closer to their typical values, some parameters will become the typical values. The remaining parameters, whose values are not close enough to the typical values, are considered the sensi tive parameters. This l eads to a reduction in the total number of parameters requiring estimation. The result is a single subset of sensitive parameters. The selected sensitive parameters are the ones that do not get pushed close enough to their typical va lues ( ) at a given . This process is repeated, while updating , until LASSO yields the number of sensitive parameters defined by the user beforehand. The algorithm of the LASSO method is summarized in Table C . 1 . 91 LASSO is typically applied to linear models (output is linear in parameters). In our case, the output is nonlinear in parameters, requiring a different method to solve the minimization problem. The Nelder - Mead Simplex algorithm was used to solve (C. 5 ). The algorithm was implemented through fminsearchbnd parameters. The estimated parameters from (C. 5 ) become compressed so we did not rely on these p arameter estimates. Instead, we used the LASSO method solely for sensitive parameter selection. Once selected, the sensitive parameters were estimated using nonlinear least squares (uncompressed estimates). This method is known as relaxed LASSO [121] . 92 Table C . 1 : The algorithm of sensitive parameter selection using the LASSO method. Input: (1) Experimental Data and Dynamical Model (2) Vector of Normalized P arameter Typical Values ( ), where the pre - normalization parameter vector is (3) Desired Optimal Number of Sensitive Parameters Output: (1) A Subset of the Sensitive Parameters 1: NumParams = 0 2 : = 0.13 3 : whil e NumParams ! = do 4 : NumParams = 0 5 : S (selection matrix) = zeros 6 : Solve (C.5) 7 : for i = 1 : do 8 : if > 0 . 005 then 9 : NumParams = NumParams + 1 1 0 : end if 1 1 : if then 1 2 : S ( NumParams,i ) = 1 1 3 : end if 1 4 : end for 1 5 : if NumParams ! = then 1 6 : 1 7 : end if 1 8 : end while 1 9 : 93 BIBLIOGRAPHY 94 BIBLIOGRAPHY [1] T. Vos et al. injuries 1990 - Lancet , vol. 380, no. 9859, pp. 2163 2196, 2013. [2] Spine (Phila. Pa. 1976). , vol. 21, no. 24, pp. 2820 2825, 1996. [3] Survey, 2002 (machine readable data file and docu Natl. Cent. Heal. Stat. Hyattsville, Maryland. http//www. cdc. gov/nchs/nhis. htm , 2002. [4] G. B. J. Andersson, T. Lucente, A. M. Davis, R. E. Kappler, J. A. Lipton, and S. Leurgans, ard care for patients with low N. Engl. J. Med. , vol. 341, no. 19, pp. 1426 1431, 1999. [5] mobilization for low back pain and neck pain: a systematic rev iew and best evidence spine J. , vol. 4, no. 3, pp. 335 356, 2004. [6] low back pain: a systematic review and meta - analysis of randomized controlled trial BMC Musculoskelet. Disord. , vol. 6, no. 1, p. 43, 2005. [7] Rheumatol. Int. , vol. 5, no. 4, pp. 145 148, 1985. [8] J. Cholewicki, H. S. Greene, G. K. Polzhofer, M . T. Galloway, R. A. Shah, and A. Radebold, J. Orthop. Sport. Phys. Ther. , vol. 32, no. 11, pp. 568 575, 2002. [9] P. W. Hodges, G. L. Moseley, A. Gabrielsson, and Exp. Brain Res. , vol. 151, no. 2, pp. 262 271, 2003. [10] M. Revel, C. Andre - patie Arch. Phys. Med. Rehabil. , vol. 72, no. 5, pp. 288 291, 1991. [11] Scand. J. Rehabil. Med. , vol. 28, no. 3, pp. 133 138, 199 6. [12] Spine (Phila. Pa. 1976). , vol. 22, no. 8, pp. 865 868, 1997. [13] - aged adults . 95 Subjects with healthy backs compared with subjects with low - Spine (Phila. Pa. 1976). , vol. 16, no. 3, pp. 325 330, 1991. [14] Footed and Externally Disturbed T wo Footed Postural Control in Patients With Chronic Low Back Pain and Healthy Control Subjects: A Controlled Study With Follow Spine (Phila. Pa. 1976). , vol. 23, no. 19, pp. 2081 2089, 1998. [15] E. - P. Takala, I. Korhonen, and E. Viikari - tural sway and stepping response among working population: reproducibility, long - term stability, and associations with Clin. Biomech. , vol. 12, no. 7 8, pp. 429 437, 1997. [16] A. Radebold, J. Cholewicki, G. K. Polzhofer, and H. of the lumbar spine is associated with delayed muscle response times in patients with Spine (Phila. Pa. 1976). , vol. 26, no. 7, pp. 724 730, 2001. [17] J. Neurophysiol. , vol. 48, no. 1, pp. 62 74, 1982. [18] N. P. Reeves, J. M. Popovich, M. C. Priess, J. Cholewicki, J. Choi, and C. J. R adcliffe, J. Biomech. , vol. 47, no. 1, pp. 44 49, Jan. 2014. [19] J. M. Popovich, N. P. Reeves, M. C. Priess, J. Cholewicki, J. Choi, and C. J. Radcli ffe, neck control: A test J. Biomech. , vol. 48, no. 3, pp. 549 554, 2015. [20] stabilization in J. Neurophysiol. , vol. 102, no. 1, pp. 496 512, Jul. 2009. [21] J. Neurophysiol. , vol. 88, no. 3, pp. 1097 1118, 2002. [22] K. J. Chen, E. A. Keshner, B. W. Peterson, and T. C. Hai J. Vestib. Res. , vol. 12, no. 1, pp. 25 33, Jan. 2002. [23] Biol. Cybern. , vol. 7 5, pp. 309 319, 1996. [24] P. van Drunen, E. Maaswinkel, F. C. T. der Helm, J. H. van Dieën, and R. Happee, - J. Biomech. , vol. 46, no. 8, pp. 1440 1446, 2013. [25] D. C. Lin and J. Biomech. Eng. , vol. 125, no. 1, pp. 132 140, Feb. 2003. [26] A. Ramadan, J. Cholewicki, C. J. Radcliffe, J. M. Popovich, N. P. Reeves, and J. Choi, ontrol during seated balancing using a physical human - J. Biomech. , vol. 64, pp. 198 205, 2017. 96 [27] Exp. Brain Res. , vol. 171, no. 2, pp. 231 250, May 2006. [28] Neurosci. Lett. , vol. 366, no. 1, pp. 63 66, 2004. [29] S. Brumagne, R. Lysens, S. Swin Spine (Phila. Pa. 1976). , vol. 24, no. 13, pp. 1328 1331, 1999. [30] - body v ibration on Clin. Biomech. , vol. 23, no. 4, pp. 381 386, 2008. [31] M. C. Priess et al. - Domain Optimal Experimental Design in Human Seated Postural J. Dyn. Syst. Meas. Control , vol. 137, no. 5, pp. 545011 545017, 2015. [32] coactivation is tuned to changes in task dynamics to improve responsiveness in a seated J. Electromyog r. Kinesiol. , vol. 25, no. 5, pp. 765 72, Oct. 2015. [33] J. Biomech. , vol. 42, no. 2, pp. 164 170, 2009. [34] U. Van Daele, F. in balance strategies between nonspecific chronic low back pain patients and healthy control Spine (Phila. Pa. 1976). , vol. 34, no. 11, pp. 1233 123 8, 2009. [35] Clin. Biomech. , vol. 30, no. 9, pp. 933 939, 2015. [36] ter of pressure trajectories, trunk Gait Posture , vol. 38, no. 4, pp. 625 630, 2013. [37] deficits in acute Spine J. , vol. 15, no. 8, pp. 1772 1782, 2015. [38] sway in unstable s Spine (Phila. Pa. 1976). , vol. 35, no. 7, pp. 812 817, 2010. [39] S. Brumagne, L. Janssens, S. Knapen, K. Claeys, and E. Suuden - European Spine Journal , vol . 17, no. 9. Springer - Verlag, Berlin/Heidelberg, pp. 1177 1184, 2008. [40] Spine (Phila. Pa. 1976). , vol. 29, no. 6, 97 pp . E107 E112, 2004. [41] J. Biomech. , vol. 33, no. 12, pp. 1733 1737, Dec. 2000. [42] C. Larivière, H. Mecheri, A. Shahvarpour, D. Gagnon, and A. Shirazi - A validity and between - day reliability of an inertial - sensor - based trunk postural stability test J. Electromyogr. Kinesiol. , vol. 23, no. 4, pp. 899 -- 907, 2013. [43] J. H. van Dieën, L. L. J. Koppes, and J. W. R. Twis Gait Posture , vol. 31, no. 1, pp. 42 46, Jan. 2010. [44] D. R. Cruise et al. - time control of stiffn ess and time - J. Biomech. , vol. 60, pp. 48 56, 2017. [45] during sagittal pelvic tilt: from trunk - on - pelvis to trun k - in - space due to vestibular and visual J. Neurophysiol. , vol. 115, no. 3, pp. 1381 1388, 2016. [46] P. van Drunen, Y. Koumans, F. C. T. van der Helm, J. H. van Dieën, and R. Happee, - ba ck stabilization due to vision, Exp. Brain Res. , vol. 233, no. 3, pp. 735 749, Mar. 2015. [47] J. B iomech. , vol. 42, no. 8, pp. 1017 1022, May 2009. [48] Clin. Biomech. , vol. 23, no. 6, pp. 735 742, Jul. 2008. [49] es about some intraclass correlation Psychol. Methods , vol. 1, no. 1, pp. 30 46, 1996. [50] Psychol. Asse ss. , vol. 6, no. 4, pp. 284 290, 1994. [51] J. Orthop. Res. , vol. 7, no. 6, pp. 849 860, 1989. [52] L. Ljung, System identification: theory for the user . PTR Prentice Hall, Upper Saddle River, NJ, 1999. [53] D. R. Billinger, Time Series: Data Analysis and Theory . SIAM, 2001. [54] - activation is associated with impaired neuromuscular Exp. Brain Res. , vol. 188, no. 3, pp. 457 463, 2008. 98 [55] 2016 American Control Conference (ACC) , 2016, pp. 5791 5796. [56] Y. Xu, J. Biomech. Eng. , vol. 132, no. 5, pp. 051004 - 1 10, 2010. [57] Ann. Biomed. Eng. , vol. 39, no. 8, pp. 2263 2273, 2011. [58] - MVC EMG normalization J . Electromyogr. Kinesiol. , vol. 11, no. 1, pp. 11 18, 2001. [59] Spine (Phila. Pa. 1976). , vol. 2 8, no. 8, pp. 834 841, 2003. [60] Clin. Biomech. , vol. 22, no. 3, pp. 266 274, 2007. [61] ntrol Intent Using Inverse ASME 2013 Dynamic Systems and Control Conference , 2013. [62] via predictive inverse linear - AA AI , pp. 3672 3678, 2015. [63] - time intent recognition of a powered Biomed. Eng. IEEE Trans. , vol. 57, no. 3, pp. 542 551, 2010. [64] G. Wasson, P. Sheth, M. Alwan, K. Granata, A Intelligent Robots and Systems, 2003.(IROS 2003). Proceedings. 2003 IEEE/RSJ International Conference on , 2003, vol. 3, pp. 2962 2967. [65] R. Huang, H. Li 2015 IEEE International Conference on Robotics and Biomimetics (ROBIO) , 2015, pp. 1465 1470. [66] H. El - Hu Eng. Appl. Artif. Intell. , vol. 50, pp. 115 124, 2016. [67] Demonstrations for IEEE Trans. Human - Machine Syst. , vol. 46, no. 4, pp. 510 521, Aug. 2016. [68] intelligent control of human prosthetic eye movements s ystem for the emotional support by a huggable pet - J. Franklin Inst. , vol. 349, 99 no. 7, pp. 2243 2267, 2012. [69] Nat. Neurosci. , vol. 7, no. 9, pp. 9 07 915, 2004. [70] K. Mombaur, A. Truong, and J. - --- Auton. Robots , vol. 28, no. 3, pp. 369 383, Apr. 2010. [71] M. C. Priess, R. Conway, J. Choi, J. M. Popovich, and C. IEEE Trans. Control Syst. Technol. , vol. 23, no. 2, pp. 770 777, Mar. 2015. [72] A. - S. Puydupin - e optimal 2012 IEEE International Conference on Robotics and Automation (ICRA) , 2012, pp. 531 536. [73] feedba J. Franklin Inst. , vol. 351, no. 3, pp. 1335 1355, 2014. [74] - agent systems J. Franklin Inst. , vol. 354, no. 4, pp. 2068 2085, Mar. 2017. [75] F. Borrelli, Constrained optimal control of linear and hybrid systems York: Springer, 2003. [76] M. Johnson, N. Aghasadeghi, and T. Bretl, Inverse optimal control for deterministic continuous - time nonlinear systems . IEEE, 2013, pp. 2906 2913. [77] Multi - Proc.~of the European Control Conference , 2013, pp. 502 510. [78] - 2012. [79] 2011 IEEE International Symposium on Intelligent Control (ISIC) , 2011, pp. 613 619. [80] IEEE Trans. Automat. C ontr. , vol. 55, no. 1, pp. 185 190, Jan. 2010. [81] M. C. Priess et al. - ASME 5th Annual Dynamic Systems and Control Conference joint with the JSME 11th Motion and Vibration Conference , 2012, pp. 503 510. [82] - dimensional model to predict loads J. Biomech. , vol. 25, no. 4, pp. 395 414, 1992. [83] H. van der Kooij, E. 100 J. Neurosci. Methods , vol. 145, no. 1, pp. 175 203, 2005. [84] al human - robot Mech. Mach. Theory , vol. 43, no. 3, pp. 253 270, 2008. [85] V. M. Zatsiorsky, Kinetics of human motion . Human Kinetics, 2002. [86] P. J. Antsaklis and A. N. Michel, A linear systems primer . Birkhäuser, 2007. [87] A. Ramadan et al. J. Biomech. Eng. , vol. 140, no. 7, pp. 074503 - 1 8, 2018. [88] P. A. Forbes, E. De Bruijn, A. C. Schouten, F. C. T. Van Der Helm, and R . Happee, - Exp. brain Res. , vol. 226, no. 1, pp. 1 14, 2013. [89] namic and Ann. N. Y. Acad. Sci. , vol. 942, no. 1, pp. 381 393, Jan. 2001. [90] muscle stiffness alone is insufficien J. Biomech. , vol. 40, no. 5, pp. 1058 1065, Jan. 2007. [91] J. Neu rosci. Methods , vol. 119, no. 1, pp. 1 14, Sep. 2002. [92] IEEE Trans. Biomed. Eng. , vol. 55, no. 1 , pp. 311 321, 2008. [93] a musculo - Comput. Methods Programs Biomed. , vol. 114, no. 3, pp. e46 e59, May 201 4. [94] Comput. Methods Programs Biomed. , vol. 114, no. 3, pp. e60 e69, May 2014. [95] J. Control. Autom. Electr. Syst. , vol. 24, no. 1 2, pp. 3 10, 2013. [96] for system Automatica , vol. 43, no. 6, pp. 993 1008, 2007. [97] Am. J. Physiol. Circ. 101 Physiol. , vol. 261, no. 3, pp. H929 H949, 1991. [98] S. Zeinali - J. Biomech. , vol. 42, no. 4, pp. 524 530, 2009. [99] nd recent developments in information J. Math. Psychol. , vol. 44, no. 1, pp. 62 91, 2000. [100] Electr. Power Syst . Res. , vol. 88, pp. 1 8, Jul. 2012. [101] Y. Yun, H. - J. Biomech. , vol. 47, no. 1, pp. 186 192, 2014. [102] A. Ramadan, J. Choi, C. J. Radcliffe, J. Cholewicki, N. P. Reeves, and J. M. Popovich, 2017 14th International Conference on Ubiquitous Robots and Ambient Intelligence (URAI) , 2017, pp. 174 1 78. [103] Automatica , vol. 46, no. 5, pp. 785 795, 2010. [104] Biol. Cybern. , vol. 57, no. 4 5, pp. 217 231, Nov. 1987. [105] IEEE Eng. Med. Biol. Mag. , vol. 22, no. 2, pp. 63 68, Mar. 2003. [106] - Neck Complex in Response J. Biomech. Eng. , vol. 125, no. 4, pp. 533 539, Aug. 2003. [107] J. R. Soc. Interface , vol. 3, no. 10, pp. 603 616, Oct. 2006. [108] Bioinformat ics , vol. 30, no. 10, pp. 1440 1448, 2014. [109] - Retest IEEE Trans. Neural Syst. Rehabil . Eng. , vol. 22, no. 5, pp. 1020 1029, 2014. [110] Experimental Paradigm to Assess Postural Stabilization: No More Movement and Not Yet IEEE Trans. Neural Syst. Rehabil. Eng. , vol. 19, no. 4, pp. 420 426, 2011. [111] - and 102 between - day reliability of trunk mechanical behaviors estimated using position - controlled J. Biomech. , v ol. 45, no. 11, pp. 2019 2022, Jul. 2012. [112] Pediatr. Phys. Ther. , vol. 22, no. 3, pp. 246 257, 2010. [113] G Individuals with and without Low Back Pain: Intratester Reliability, Clinical Applicability, J. Orthop. Sport. Phys. Ther. , vol. 32, no. 7, pp. 327 335, Jul. 2002. [114] - Automatica , vol. 41, no. 3, pp. 393 438, Mar. 2005. [115] 45t h IEEE Conference on Decision and Control , 2006, pp. 163 168. [116] Meas. Sci. Technol. , vol. 9, no. 6, pp. 864 876, Jun. 1998. [117] J. R. Stat. Soc. Ser. B , vol. 58, no. 1, pp. 267 288, 1996. [118] J. R. Stat. Soc. Ser. B (Statistical Methodol. , vol. 73, no. 3, pp. 273 282, 2011. [119] G. Golub and C. Van Loan, Matrix computations , 4th ed. Baltimore: JHU Press, 2013. [120] IEEE Trans. Signal Process. , vol. 57, no. 2, pp. 451 462, 2009. [121] Comput. Stat. Data Anal. , vol. 52, no. 1, pp. 374 393, 2007.