Validation of Dynamic Bayesian Optimization for a Non-Stationary Human-in-the-Loop Optimization Problem.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice IEEE Robot Autom Lett . Author manuscript; available in PMC: 2026 Apr 11. Published in final edited form as: IEEE Robot Autom Lett. 2026 Feb 16;11(5):5733–5740. doi: 10.1109/lra.2026.3665072 Search in PMC Search in PubMed View in NLM Catalog Add to search Validation of Dynamic Bayesian Optimization for a Non-Stationary Human-in-the-Loop Optimization Problem GilHwan Kim GilHwan Kim 1 Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA Find articles by GilHwan Kim 1 , Fabrizio Sergi Fabrizio Sergi 2 Department of Mechanical and Biomedical Engineering, University of Delaware, Newark DE, 19713, USA Find articles by Fabrizio Sergi 2 Author information Article notes Copyright and License information 1 Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA 2 Department of Mechanical and Biomedical Engineering, University of Delaware, Newark DE, 19713, USA ✉ Email: [email protected] Issue date 2026 May. All rights reserved, including rights for text and data mining, and training of artificial intelligence and similar technologies. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. PMC Copyright notice PMCID: PMC13056379  NIHMSID: NIHMS2161189  PMID: 41953283 The publisher's version of this article is available at IEEE Robot Autom Lett Abstract Human-in-the-Loop Optimization (HILO) has demonstrated efficacy in achieving a plethora of assistive or augmentative effects. However, conventional optimizers such as Bayesian optimization (BO) do not account for the non-stationary aspects of the human-robot system, and may thus be limited in the domains of robot-assisted training or rehabilitation. In this study, we implemented HILO using dynamic Bayesian optimization (DBO) to define the optimal value of a single control parameter to target a desired effect in propulsion mechanics, specifically in the maximum hip extension (HE) angle during stance. Sixteen healthy participants received unilateral hip torque pulses via a hip exoskeleton for about 15 minutes while walking on a treadmill. Exoskeleton torques were determined via HILO using DBO or BO. Unbeknownst to the optimizers, treadmill speed was gradually increased to amplify the non-stationary behavior of the system. Experimental results revealed that DBO outperformed BO in the later stages of training, leading to improved cost and to final values of torque closer to the true value. The primary difference between optimizers arose from their ability to model the history-dependent relationship between applied torque and output. BO associated later deviations between expected and measured output to process noise, while DBO developed a more accurate predictive model of the input-output relationship that gave proper weighing to past datapoints for making predictions. Index Terms—: Human Factors and Human-in-the-Loop, Machine Learning for Robot Control, Wearable Robotics I. INTRODUCTION Robot-assisted gait training offers the advantage of delivering consistent and repeatable mechanical actions during walking, compared to standard training methods [ 1 ], [ 2 ]. In this context, training effects can be achieved by applying an appropriate level of robotic assistance to the user. Typically, the level of robotic assistance is determined by identifying the user’s input-output relationship or by relying on group-averaged responses. However, estimating individual responses using group-averaged model eliminates inter-individual variability during construction of a response model, which can result in sub-optimal level of assistance for specific participants whose optimal assistance is substantially different from the group average. Alternatively, identifying an accurate input-output relationship model for an individual requires a large number of observations, a number that increases further when the robotic device has large number of control parameters, with a consequent high risk of inducing fatigue. Therefore, efficiently gathering sufficient observations from individuals to construct such input-output relationship model remains a primary challenge in robot-assisted gait training. To simultaneously account for variability of individual responses and reduce the number of required observations to generate a specific individual response model, human-in-the-loop optimization (HILO) has been introduced [ 3 ], [ 4 ]. By iteratively tuning control parameters during training, HILO enables efficient identification of a individual-specific response model and corresponding optimal input, using different optimization methods that include gradient descent, Bayesian optimization (BO), and covariance matrix adaptation evolution strategy [ 5 ], [ 6 ], [ 7 ]. Previous implementations of HILO in human-robot interaction scenarios targeting walking have primarily targeted the minimization of metabolic cost. While effective for human augmentation purposes, this approach requires relatively longer bouts of walking to establish a reliable estimation of the cost function at each iteration compared to other metrics such as propulsion mechanics or walking speed. Moreover, metabolic cost is not a primary outcome in motor rehabilitation, which may target more directly aspects of assisted or unassisted biomechanical performance, including walking speed or propulsion mechanics. Although HILO has the potential to deal with non-stationary system, such as changes in individual response due to neuromotor adaptation or learning, optimizer forms used in prior implementations have relied on the assumption that the human-robot system remains stable over time or across interventions. This assumption disregards the influence of training history on the response of an individual, leading the optimizer to expect similar outcomes when applying identical torque/force assistance regardless of the timing or prior training. Such assumption limits the ability of previous implementations of HILO to account for the changes in participant input expected during training. Dynamic Bayesian optimization (DBO) was introduced in previous studies to support identification of optimal inputs in problems where the cost function could change over time. In DBO, this functionality is achieved by explicitly accounting for the time of inputs and output in the model covariance function [ 8 ], [ 9 ]. While DBO has been validated in prior studies that aimed to identify the optimal input for time-varying systems in simulations [ 8 ], [ 9 ], [ 10 ], DBO has not been tested in the context of HILO for robot-assisted gait training. Previous studies showed the implementation of DBO in HILO simulations involving neuromotor learning [ 11 ] and in a pilot experimental study involving robot-assisted gait training [ 12 ], showing performance that was similar to or exceeded that of standard BO. For accurate comparison of controller performance across different implementations and thus HILO runs, it is critical to compare outcomes not only during optimization, but also in separate cross validation sessions. In absence of such validation trials, hypothetical differences in performance measured across optimizers may arise due to the inherently different exploration-exploitation trade-offs pursued by different optimizers, a trade-off that is impossible to control over different optimizer implementations. In this study, we implemented HILO on a hip exoskeleton using DBO to define the optimal value of a single control parameter targeting a desired effect in propulsion mechanics over a single session lasting approximately 15 minutes, and compared performance between DBO and a standard BO algorithm. We hypothesized that DBO would be more suitable for identifying the optimal solution in a training protocol featuring non-stationary aspects. In a non-stationary problem, the target outcome measure would need to exhibit substantial variability due to neuromotor adaptation and learning, or the participant would have to experience training for a sufficiently long time to result in a significant change in individual response. To establish the feasibility of DBO for such a non-stationary problem in a convenient experimental model, while avoiding fatigue, we decided to amplify the non-stationary aspect of the human-robot system by gradually increasing treadmill speed during training. To compare the performance across optimizers, we incorporated validation iterations in which the optimizer tested its best estimated input value, expected to minimize cost function value, an experimental design choice that was not used in previous studies using DBO [ 12 ]. Such an experimental design enabled a properly controlled comparison of the optimizers’ performance. Moreover, because the true optimal solution was expected to change as treadmill speed increased, this allowed us to evaluate each optimizer’s ability to accurately model the participant’s input–output relationship. Experimental results showed that the DBO outperformed BO during the later stage of training. II. Methods A. Dynamic Bayesian Optimization DBO is a modified version of BO. A detailed description of the method is provided in previous studies [ 10 ], [ 11 ], but a general overview of the DBO framework is provided in this section for the purposes of following the experimental work presented in this paper. In general, BO uses a Gaussian process (GP) model to estimate outputs corresponding to unobserved input values. Based on this probabilistic estimation model, BO iteratively searches for the input value that minimizes a cost function that is evaluated under measurement noise [ 13 ]. The GP model is regenerated iteratively when new observation is provided. The input value (or set of values) to be tested at the next iteration is determined by maximizing the acquisition function. This acquisition function scores candidate input values based on the predicted performance and the related uncertainty, so to balance exploration and exploitation. The conventional BO framework assumes that the system is stationary; i.e., BO expects similar outputs when inputs are repetitively applied to the system. However, this assumption becomes inaccurate when the system includes a dynamic component, as it is the case when human factors are involved due to neuromotor learning and adaptation. In standard BO, the GP model relies on a stationary covariance function k ( u , u ′) between inputs u and u ′, which leads to similar predictions for repetitive inputs. To account for the dynamics of the system, in DBO, the GP model was modified by introducing the time variable t into the covariance function, resulting in k (( u , t ), ( u ′, t ′)), where t and t ′ represent the time instances associated with the inputs u and u ′ respectively. This covariance function has been simplified in previous implementations of DBO [ 9 ], [ 14 ] by assuming that the components of the covariance function are separable, allowing the overall covariance function k (( u , t ), ( u ′, t ′)) to be expressed as a multiplication of a static k ( u , u ′) and dynamic component k ( t , t ′): k ( u , t ) , u ′ , t ′ = k u u , u ′ ⋅ k t t , t ′ = k u u , u ′ ⋅ α t − t ′ . (1) In this work, k u ( u , u ′) was implemented as the automatic relevance determination squared exponential function [ 15 ]. Hyperparameter α ∈ (0, 1] defines the time scale in the sense that the covariance between measurements collected at times t and t ′ is scaled by a factor α | t − t ′ | ( Fig. 1b ). Fig. 1. Open in a new tab (a) Experimental setup including the active pelvis orthosis (APO), instrumented treadmill, motion capture system, and safety harness. (b) Dynamic component of the covariance function. Each line represents the covariance function with different hyperparameter α . (c) Hip torque profile in the APO Gait Cycle domain where 0% is peak of hip flexion. Heel-strike and toeoff timings are shown as distributions, with the mean and standard deviation indicated in blue and red, respectively. The control parameters defining the torque profile (duration: D , timing: P , and amplitude: T ) are labeled. B. Experimental Methods 1). Equipment: A lower extremity robotic exoskeleton, the Active Pelvis Orthosis (APO), was used to provide torque pulses to participants while walking on the treadmill ( Fig. 1a ). The APO (IUVO, Pisa, Italy) is a bilateral hip exoskeleton with a weight of 6 kg including battery [ 16 ]. An onboard encoder in the APO was used to measure the hip angle with a resolution of 0.015 deg and sampling rate of 100 Hz. Based on the hip joint angle measurements, the APO estimated the participant’s gait phase via adaptive oscillators [ 17 ]. These estimates were computed individually for the right and left legs. Using the predicted gait phase and hip kinematics, hip torques were applied in the form of a Gaussian function shaped pulse profile, defined by three control parameters: torque pulse amplitude ( T ), duration ( D ), and peak timing ( P ) ( Fig. 1c ). In this experimental study, a single hip torque pulse was applied to the participant’s right hip during each stride. The hip torque pulse was applied near the end of the stance phase, with the duration ( D ) of 18% of the gait cycle (full width at half-maximum: 7.07%), and the timing of peak torque amplitude ( P ) of 65% in the APO gait cycle (0% at peak of hip flexion), resulting in the peak torque being applied at the end of the stance phase. The torque pulse amplitude ( T ) was iteratively updated in real-time via Simulink (MathWorks Inc., Natick, MA, USA). The Simulink controller communicated with the APO via a Serial Peripheral Interface connection. To ensure proper execution of the updated torque profile, the new amplitude command was sent to APO during the swing phase, so that the new torque pulses would be applied in the following stride during the stance phase. Analog/force plate data from both the right and left side were collected via an instrumented split-belt treadmill (Bertec Corp., Columbus OH, USA) that houses two independent force plates located beneath each treadmill belts. A Simulink-based controller was used to collect signals from the force plates and to send speed commands to the treadmill at 400 Hz. Leg kinematics were collected in real-time via a 10-camera motion capture system; T40-S (Vicon Motion Systems Ltd, Oxford, UK) with primary goal of tracking the trailing limb angle (TLA) of both legs. Cameras were used to track four retroreflective markers placed on the participant: right and left side of malleoli, and on both side of the rotational axis of APO which was aligned with the participant’s hip joint. Marker locations were recorded at a sampling rate of 100 Hz. For safety, participants were secured using a harness system (Solo-Step Inc., North Sioux City, SD, USA) connected to an overhead railing track during experiment. 2). Study Participants: 16 healthy participants (8 females; age: 25.6 ± 2.1 yrs, height (mean ± std): 170.3 ± 7.0 cm, and mass: 73.4 ± 10.3 kg) were recruited for this experiment (protocol no. 1755609–7, approved by the University of Delaware Institutional Review Board). Each participant conducted two training sessions on the same day, with a rest period of at least fifteen minutes between sessions to prevent fatigue and minimize potential remaining training effects of the first training session. The two sessions differed in the optimization algorithm used (DBO vs. BO), and the order of optimizer used for a participant was pseudo-randomized and evenly balanced across participants. 3). Pre-training Procedures: The experimental protocol consisted of a pre-training phase followed by two training sessions. During the pre-training phase, retroreflective markers were placed on the participant and APO. The participant donned the APO, aligning the rotational axis of the APO with the participant’s hip axis in sagittal plane. Once the participant confirmed comfort, the device and safety harness was secured. The participant then walked on the treadmill at low speed (0.55 m/s) for approximately 2 minutes with the APO in transparent mode, where the APO sought to minimize the interaction force/torque between the device and user. Treadmill speed was then gradually increased by 0.05 m/s every 5 strides, maintaining a constant speed over such intervals, until reaching 1.2 m/s. The maximum hip extension (HE) angle of the right leg during stance phase from the last 3 strides at each treadmill speed was collected and used to estimate the linear relationship between HE and treadmill speed for each participant. Based on the estimated linear model, treadmill speed required to induce the target change (5 degrees) in participant’s HE was calculated. Immediately after collecting HE data at different treadmill speed, the treadmill speed was set to 1 m/s, and the familiarization session was conducted. During this session, participant experienced a series of torque pulses (extension: 1, 3, 5, 7, 9 Nm, flexion: 1, 3, 5 Nm) applied in ascending order for ten strides each (duration D = 18% and timing P = 65%, as tested during optimization). Following the pre-training phase, participants were given a rest period approximately 5 minutes. 4). Treadmill Speed Increment During Training: Given their formulation, DBO and BO may perform differently if the system response is non-stationary [ 9 ], [ 10 ], [ 11 ]. In practical training applications, such difference may emerge when the target outcome is highly non-stationary or when participants are exposed to training over a sufficient number of iterations, possibly over multiple repeated sessions. However, testing a new algorithm on a multi-session study conducted in natural walking conditions would be an inefficient experimental design that would be overly burdening for participants. As an initial step towards validating the feasibility of a optimizer algorithm in HILO that can account for the non-stationary properties of the human response to training, we opted to intentionally induce a controlled amount of non-stationary behavior in the measured response. Because most metrics of walking mechanics are speed-dependent, we opted to gradually increase treadmill speed during training, a condition that was completely unknown to the optimizer. As treadmill speed increased, participants altered their propulsion mechanics [ 18 ], naturally choosing to walk with longer strides, greater cadence, and greater peak hip extension. In this condition, the optimizer would be required to realize that the optimal input in the initial stages of the experiment would be substantially different from the one that would be required in the late stages of the experiment. Moreover, because the treadmill speed conditions were consistent across optimizer/session for a given participant (BO vs. DBO), this experimental design allowed us to directly compare outcomes between optimizers during training. The overall speed range for the training session was determined based on the outcomes of the pre-training session. The maximum (final) speed was fixed at 1.2 m/s, and each participant’s initial speed (low speed) was calculated by subtracting the individual speed range from 1.2 m/s. Based on the observation that maximum HE increases with walking speed, the speed range was derived from a linear model based on pre-training session data to produce a 5-deg change in the target outcome. If the calculated low speed fell below 0.55 m/s, it was set to 0.55 m/s. The speed range was defined such that the ideal torque input was expected to decrease progressively over the course of training, ultimately reaching zero torque amplitude by the final iteration (80 th iteration). Treadmill speed was initially set to the low speed for baseline and the first iteration. Speed was then increased by a fixed increment at each iteration to reach 1.2 m/s by the 80 th iteration. Treadmill speed remained constant within an iteration and was updated right after cost function evaluation. 5). Training Procedures: The training session consisted of a baseline phase followed by a torque intervention phase. During the baseline phase, participants walked at low speed in transparent mode of APO, without any torque pulses, for 110 strides. The torque intervention phase started right after. During this phase, torque pulses were applied to the right hip joint at the end of the stance phase. The torque intervention session consisted of 80 continuous optimizer iterations, where each iteration consisted of either 6 or 7 strides, depending on the time it took for the optimizer to identify the input for the next iteration. Within each iteration, the APO applied the same torque pulse continuously from the first stride until the optimizer updated the control parameters for the next iteration. Based on the participant response, the optimizer updated the control parameter T . HE was measured over the first 5 strides of each iteration, and the average value of these 5 measurements was used to calculate the cost function value c ( i ), which was sent to the optimizer. c ( i ) was defined as the absolute difference between the desired HE and the average HE measured at that iteration. After completing the first training session, the participant removed all equipment, including exoskeleton and harness, and was asked to rest for at least 15 minutes before starting the second training session with the different optimizer. 6). Optimizer Setup: The objective of both optimization methods (DBO and BO) used in each training session was to increase maximum HE of the right leg during the stance phase by 5 degrees relative to baseline. The baseline HE ( HE BL ) was calculated by averaging maximum HE during stance phase over the last 20 strides of the baseline session. At each iteration i , HE ( HE ( i )) was calculated as average of the first 5 consecutive strides within that iteration. Accordingly, the cost function c ( i ) was defined as absolute difference between the HE target ( HE BL + 5°) and measured HE values, i.e., c ( i ) = | HE BL + 5° − HE ( i )|. During each iteration, the participant experienced the same torque input for each stride until the optimizer determined the new control input to test. This update took an additional 1–2 strides in this experiment. This experiment featured a single-parameter optimization method, targeting the parameter of torque pulse amplitude T . The range of T was from −5 Nm (maximum flexion) to 9 Nm (maximum extension). The input values for the first three iterations were predefined as 5 Nm, 7 Nm, and 3 Nm, respectively. The acquisition function was expected improvement, selected based on its highest performance in previous simulation results [ 19 ]. The exploration-exploitation ratio (e-ratio) was set to 0.1, favoring exploitation. This optimizer value was set based on previous simulation results [ 11 ], which showed no significant difference between different e-ratios ranging from 0.1 to 0.8 when using the same optimizer (DBO or BO). Pilot testing indicated that optimization runs with an e-ratio of 0.1 were less overly variable and more intuitive for the user. A total of 80 iterations were performed in each training session, with validation iterations conducted every 10 iterations. This procedure allowed to validate the accuracy with which the optimizer was able to form a model of the response between input and output at the estimated optimal input, and compare between optimizers without confounds related to different instances of exploration pursued by different optimizers or optimization runs. The GP model was constructed at each iteration by selecting hyperparameter values, including α , based on past observations, including validation iteration, to maximize the log-likelihood. During validation iterations, the input with the minimum upper confidence bound of the cost function value was selected. 7). Statistical Analysis: Individual participant responses at each validation iteration were collected and analyzed, and paired comparisons of outcome measures were used to assess the performance of the two optimizers. Outcomes subjected to statistical analysis included the following primary outcomes, which were expected to be directly influenced/modulated by the optimizer: change in HE relative to baseline, cost function value, and applied torque. Moreover, secondary outcomes such as cadence, normalized peak anterior ground reaction force (nAGRF), normalized propulsive impulse (nPI) integrated over the positive portion of the AGRF curve, the maximum trailing limb angle of the right leg during stance (max. TLA), were also analyzed in a similar fashion to assess generalization to other components of walking not directly targeted by the optimizer. Finally, the prediction error of the best-estimated response at the cross-validation iterations was also calculated to assess the accuracy of the model developed by BO and DBO longitudinally across an experiment. To assess the need of parametric testing, data distributions for each outcome across the sixteen participants at each validation iteration were evaluated for normality using the Shapiro–Wilk test. Pairwise t-tests were conducted to compare outcome differences between optimizers at each validation iteration. In addition, two two-way repeated-measures ANOVAs were performed to examine the effects of optimizer type (DBO vs. BO) and iteration number (20 th vs. 80 th ) on cost function values and applied torque, with the 20 th and 80 th iterations representing early and late training phases, respectively. A repeated-measures design was chosen to account for within-subject variability, as the same participants experienced both optimizers across multiple iterations. When significant main effects or interactions were identified, post-hoc pairwise t-tests were performed to compare each combination of two factors: optimizer and iteration number. We reported effect sizes by computing the standardized mean differences (Cohen’s d z coefficient for the paired t-tests, or Cohen’s d coefficient for comparisons to a constant) and partial eta squared η p 2 for the repeated-measures ANOVA. III. Results A. Relationship between Treadmill Speed and Maximum Hip Extension Angle during Stance Phase The relationship between maximum HE during stance and treadmill speed was assessed during the pre-training phase ( Fig. 2 ). As shown in the figure, all participants exhibited a positive linear relationship between HE and treadmill speed. Based on these observations, the baseline treadmill speed for the main training sessions was selected to produce a 5 degrees change in HE. Fig. 2. Open in a new tab Individual participants’ HE values are shown with respect to treadmill speed during the pre-training phase. Each color represents a different participant, circular markers indicate the mean measured HE at each speed, and solid lines show the corresponding linear fits based on observed data. B. Performance difference between DBO and BO Group averaged cost function values at each iteration are shown in Fig. 3 . Results from each optimizer (DBO and BO) were subject to paired comparisons at the validation iterations, where each optimizer tested the best estimated input value. As shown in Fig. 3 , the performance of the two optimizers was indistinguishable until the final iterations (specifically, iterations 70 and 80). At the 70 th iteration, DBO resulted in a cost function value of 1.5 ± 0.25° (mean ± s.e.m.), while the cost for BO was 3.0 ± 0.3°, with a significant difference between DBO and BO (1.5 ± 0.4°, p = 0.002, d z = 0.999). The difference between DBO and BO remained significant also at the final validation iteration (1.7 ± 0.5°, p = 0.002, d z = 0.966), with DBO resulting in a cost function value of 1.8 ± 0.35°, and BO resulting in a cost function value of 3.5 ± 0.33°. Fig. 3. Open in a new tab Group level cost function values during the experiment. The solid line indicates the mean cost function value for 16 participants at each iteration, and the shaded area extends by one standard deviation. The cyan dashed line indicates validation iterations where the estimated best input was tested, corresponding cost function values are shown with a circle marker. An asterisk represents significant difference between DBO and BO (p < 0.05). Group averaged torque amplitude at each iteration is shown in Fig. 4 . In BO, the mean torque amplitude remained relatively constant from the 10 th iteration to the end of experiment. In contrast, DBO showed a progressive decrease in average torque amplitude over iteration, a trend closer to the one expected due to the fact that speed progressively increased during training. Significant differences in torque amplitude between DBO and BO were observed at the 70 th (1.0 ± 0.4°, p = 0.019, d z = 0.680) and 80 th (1.4 ± 0.4°, p = 0.001, d z = 1.023) validation iteration. At the 70 th iteration, DBO applied a significantly lower torque amplitude (1.6 ± 0.4 Nm) compared to BO (2.6 ± 0.25 Nm). The difference remained significant at the 80 th iteration (DBO: 1.2 ± 0.33 Nm, BO: 2.6 ± 0.25 Nm). Fig. 4. Open in a new tab Group level torque amplitude at each iteration. The black dashed line indicates the zero value, which is the expected value that the optimizer should reach at the end of the experiment. The solid line indicates the mean and the shaded area represents one standard deviation at each iteration. Torque amplitude is positive in the direction of hip extension. The two-way ANOVAs were conducted to analyze longitudinal changes in cost function value and torque. The ANOVA for cost highlighted that both optimizer type ( p < 0.0001, η p 2 = 0.404 ) and iteration number ( p = 0.0004, η p 2 = 0.191 ) significantly affected cost function value, as did the interaction between optimizer type and iteration number ( p = 0.004, η p 2 = 0.132 ). As shown in Fig. 5 , post-hoc tests indicated that cost increased slightly between iterations for both DBO and BO (DBO: 0.9 ± 0.4°, BO: 2.4 ± 0.3°, but the increase in cost between iterations was greater for BO than it was for DBO (1.5 ± 0.5°, p = 0.005, d z = 0.851). The ANOVA for torque highlighted a significant effect of both factors (optimizer type: p = 0.011, η p 2 = 0.104 ; iteration number: p = 0.001, η p 2 = 0.166 ), but no significant interaction effect. Similarly as for cost, the amount of applied torque decreased slightly between iterations for both DBO and BO (DBO: −1.2 ± 0.4 Nm, p = 0.008, d z = 0.790; BO: −0.2 ± 0.06 Nm, p = 0.022, d z = 0.660), and remained significantly different from zero at iteration 80 (DBO: 1.2 ± 0.3 Nm, p = 0.003, d = 0.916; BO: 2.6 ± 0.3 Nm, p < 0.0001, d = 2.724), however the decrease in torque was significantly greater in magnitude for DBO compared to BO (1.1 ± 0.4 Nm, p = 0.02, d z = 0.650). All individual responses for cost function value and applied torque from the experiment are shown in Fig. 7 . Fig. 5. Open in a new tab Paired comparisons broken down by iteration number and optimizer type. Colored dots indicate each participant’s data. An asterisk represents a significant difference based on paired comparisons. Fig. 7. Open in a new tab Individual participant responses. The solid line shows the cost function value at each iteration, and dashed line indicates the corresponding torque amplitude applied to the participant. Changes in group-averaged maximum HE at each iteration relative to baseline are shown in Fig. 6 . Consistent with the cost function results ( Fig. 3 ), DBO produced significantly different outcomes compared to BO at iteration 70 (1.6 ± 0.5°, p = 0.004, d z = 0.883) and 80 (1.9 ± 0.5°, p = 0.001, d z = 1.033). DBO maintained HE values close to the desired increase relative to baseline (70 th iteration: 6.3 ± 0.3°, 80 th iteration: 6.7 ± 0.4°). In contrast, BO resulted in a larger deviation from the target 5° increase, with observed values of 7.9 ± 0.4° at the 70 th iteration and 8.5 ± 0.3° at the 80 th iteration. Fig. 6. Open in a new tab Group level maximum HE change relative to baseline at each iteration. The solid line indicates the mean and the shaded area represents one standard deviation at each iteration. The optimizers have a target of a 5-deg increase in maximum HE relative to baseline (black dashed line). C. Training-induced changes in secondary outcomes of propulsion mechanics Group averaged results for secondary outcomes quantifying propulsion mechanics of the leg exposed to exoskeleton torque are reported in Fig. 8 . As expected, all of these outcomes increase as the imposed treadmill speed increases. Unlike the primary metrics reported above, there was no significant difference between the results from DBO and BO. Fig. 8. Open in a new tab Group level change in secondary outcomes at each iteration relative to baseline. The solid line indicates the mean of all participant’s average outcome change at each iteration ( n = 16 for cadence, nPI, and nAGRF, and n =13 for max. TLA), and the shaded area represents one standard deviation. IV. Discussion The main goal of this paper was to validate DBO in a non-stationary HILO problem in the context of robot-assisted training of propulsion mechanics. An experiment involving a robotic hip exoskeleton and sixteen healthy participants was conducted to compare the outcomes of DBO with those of a standard BO algorithm. Experimental results indicate that DBO exhibited improved performance relative to BO, primarily because it was more capable than BO of accounting for non-stationary aspects in the human-robot system. This is clearly demonstrated in Fig. 3 and 6 , where the results from DBO show smaller values of cost in the later stages of the walking experiment, despite the imposed increase in treadmill speed. In contrast, the results from BO exhibit a form of drift of the cost function value as iterations progress. Statistical analysis revealed a significance difference between DBO and BO results starting from the later stages of training (iteration 70), which persisted through the final validation iteration. The non-stationary aspect of human-robot system in this experiment was amplified by gradually increasing treadmill speed, with adjustments individualized based on each participant’s pre-training responses. Since participants were expected to require progressively less torque as training proceeded, ideally reaching zero torque at the end of training, the desired torque amplitude at each iteration was expected to decrease during training. As shown in Fig. 4 , group average torque amplitudes from DBO and BO diverged significantly at the 70 th and 80 th validation iterations. Both optimizer achieved the target HE change early in the training (around iteration 10 and 20, Fig. 6 ). This early achievement resulted in BO applying similar torque amplitude longitudinally across training (20 th iteration: 2.8 ± 0.2 Nm, 80 th iteration: 2.6 ± 0.3 Nm). In contrast, DBO escaped from keeping such a constant torque amplitude pattern around the 40 th iteration and resulted in a larger overall decrease in torque amplitude during training compared to BO (20 th iteration: 2.4 ± 0.2 Nm, 80 th iteration: 1.2 ± 0.3 Nm). As shown in Fig. 5 , both DBO and BO showed a significant decrease in torque amplitude from the early (20 th iteration) to late stage of training (80 th iteration), but the reductions were insufficient to reach zero torque at the last validation iteration. As visible in Fig. 6 , the applied torques resulted in the outcome measure exceeding the target, implying that both optimizers consistently estimated the optimal input to be greater than the true value. This phenomenon is likely due to the fact that both optimizers have some unavoidable delay in estimating the GP model accurate at the current iteration. Additionally, asymmetry in the participant’s response may occur due to the unilateral torque intervention, which could cause the true optimal input to be greater than zero. Given the significant differences between DBO and BO, such delay was considerably smaller in DBO compared to BO. However, the observation that both optimizers applied a torque that was significantly greater than zero (expected value based on the pre-training session) at iteration 80 th indicates that such a delay was non-zero for both optimizers. The performance difference between DBO and BO reflects their differing ability to account for the relevance of past observations to current conditions, which is closely tied to the model used by each optimizer to define an optimal input at each iteration during training. An example describing the training history for a selected participant is shown in Fig. 9 . Initially (e.g., iterations 10–20), both optimizers applied small torque amplitude required to achieve the target increase in HE. As training continued, because the treadmill speed increased, the amount of torque necessary to achieve the same outcome should have decreased over iteration, leading to DBO and BO observing outcomes that differed from expected values. While BO interpreted those differences to arise from the unstructured variability of the human-robot system (process noise + measurement noise), leading to an increase in the confidence interval of the estimated GP model at later iterations, DBO was able to correctly assign the source for such a mismatch to the fact that there is a time-dependent co-variance between measurements collected at different time points. As such, DBO was able to correctly down-weigh the observations collected in the early phases of the experiment, in favor of more recent observations that led to better model predictions, resulting in the application of torque inputs that were closer to the expected value of zero, leading to a smaller cost compared to BO. While this is intuitively visible in the single-subject plot ( Fig. 9 ), the prediction errors obtained via DBO and BO are significantly different also at a group level ( Fig. 10 ). The superior performance of DBO in handling non-stationary systems with continuous changes demonstrates its potential for applications that require ongoing adjustments. This makes it particularly promising for tasks involving learning and adaptation, fatigue management, or assist-as-needed control paradigms, where continuous intervention and response changes are essential. Fig. 9. Open in a new tab Gaussian process models from participant 1 during training. Models were constructed using the response history before the validation sessions at iterations 10, 20, 70, and 80. The solid line shows the estimated mean cost function across torque values, with dashed lines indicating the 5th to 95th percentile confidence interval. Circles and asterisks denote the estimated mean and actual response at each validation session, respectively. Dots represent past responses, with greater transparency indicating earlier iterations. Fig. 10. Open in a new tab Group-averaged prediction error at each validation iteration. Prediction error was calculated based on the absolute difference between estimated cost value from optimizer and measured cost value. Markers indicate the mean, and the shaded area indicates one standard deviation. V. Conclusion This study validates the use of DBO in a non-stationary HILO problem in the context of robot-assisted training of propulsion mechanics. Performance differences between DBO and BO were evaluated using validation iterations, where the estimated best solutions were tested. The results showed that DBO outperformed BO in the later stage of training. There are some limitations to this study. In real-world application for training propulsion mechanics, the dynamic changes observable within a similar time-span (10–15 minutes of robot-assisted gait training) may be less pronounced than those observed in this study, where a non-stationary component was intentionally introduced by increasing treadmill speed. Consequently, the performance difference between optimizers may be smaller in other settings. It is worth noting that some domains of application of human-robot interaction involve sessions lasting hours and include the repetition of multiple sessions [ 20 ], making this method potentially a good candidate for extension to scenarios of long-term human-robot interaction. Additionally, the cost function used in this study, defined as the absolute difference between the desired and current change in HE, is not expected to be a smooth function since it exhibits a discontinuity of the derivative at the optimal input. This discontinuity is also expected to widen the confidence intervals of the GP model. This characteristic limits the optimizer performance to precisely model the cost function relative to torque using GP modeling. In future work, this study will be expanded to a larger number of control parameters targeting biomechanically relevant outcomes such as walking speed. A key challenge in this expansion is the increased computational complexity, potentially resulting in longer iteration times and reduced GP modeling accuracy. Acknowledgments This work is supported in part by NIH-R01HD111071 and in part by American Heart Association under Grant TPA – 947225. Contributor Information GilHwan Kim, Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA. Fabrizio Sergi, Department of Mechanical and Biomedical Engineering, University of Delaware, Newark DE, 19713, USA. References [1]. Esquenazi A and Packel A, “Robotic-assisted gait training and restoration,” American journal of physical medicine & rehabilitation, vol. 91, no. 11, pp. S217–S231, 2012. [ DOI ] [ PubMed ] [ Google Scholar ] [2]. Mehrholz J, Thomas S, Kugler J, Pohl M, and Elsner B, “Electromechanical-assisted training for walking after stroke,” Cochrane database of systematic reviews, no. 10, 2020. [ Google Scholar ] [3]. Felt W, Selinger JC, Donelan JM, and Remy CD, ““body-in-the-loop”: Optimizing device parameters using measures of instantaneous energetic cost,” PloS one, vol. 10, no. 8, p. e0135342, 2015. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [4]. Slade P, Atkeson C, Donelan JM, Houdijk H, Ingraham KA, Kim M, Kong K, Poggensee KL, Riener R, Steinert M et al. , “On human-in-the-loop optimization of human–robot interaction,” Nature, vol. 633, no. 8031, pp. 779–788, 2024. [ DOI ] [ PubMed ] [ Google Scholar ] [5]. Koller JR, Gates DH, Ferris DP, and Remy CD, “‘body-in-the-loop’optimization of assistive robotic devices: A validation study.” in Robotics: Science and Systems, vol. 2016, 2016, pp. 1–10. [ Google Scholar ] [6]. Zhang J, Fiers P, Witte KA, Jackson RW, Poggensee KL, Atkeson CG, and Collins SH, “Human-in-the-loop optimization of exoskeleton assistance during walking,” Science, vol. 356, no. 6344, pp. 1280–1284, 2017. [ DOI ] [ PubMed ] [ Google Scholar ] [7]. Ding Y, Kim M, Kuindersma S, and Walsh CJ, “Human-in-the-loop optimization of hip assistance with a soft exosuit during walking,” Science robotics, vol. 3, no. 15, p. eaar5438, 2018. [ Google Scholar ] [8]. Aglietti V, Dhir N, González J, and Damoulas T, “Dynamic causal bayesian optimization,” Advances in Neural Information Processing Systems, vol. 34, pp. 10 549–10 560, 2021. [ Google Scholar ] [9]. Bogunovic I, Scarlett J, and Cevher V, “Time-varying gaussian process bandit optimization,” in Artificial Intelligence and Statistics. PMLR, 2016, pp. 314–323. [ Google Scholar ] [10]. Nyikosa FM, Osborne MA, and Roberts SJ, “Bayesian optimization for dynamic problems,” arXiv preprint arXiv:1803.03432, 2018. [ Google Scholar ] [11]. Kim G, Chishty HA, and Sergi F, “Using dynamic bayesian optimization to induce desired effects in the presence of motor learning: a simulation study,” bioRxiv, 2024. [ Google Scholar ] [12]. Kim G, Chishty H, and Sergi F, “Pilot testing of dynamic bayesian optimization for exoskeleton-assisted training of propulsion mechanics,” in 2025 International Conference On Rehabilitation Robotics (ICORR). IEEE, 2025, pp. 906–911. [ Google Scholar ] [13]. Brochu E, Cora VM, and De Freitas N, “A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning,” arXiv preprint arXiv:1012.2599, 2010. [ Google Scholar ] [14]. Gneiting T, Genton MG, and Guttorp P, “Geostatistical space-time models, stationarity, separability, and full symmetry,” Monographs On Statistics and Applied Probability, vol. 107, p. 151, 2006. [ Google Scholar ] [15]. Van Gestel T, Suykens J, De Moor B, and Vandewalle J, “Automatic relevance determination for least squares support vector machine regression,” in IJCNN’01. International Joint Conference on Neural Networks. Proceedings (Cat. No. 01CH37222), vol. 4. IEEE, 2001, pp. 2416–2421. [ Google Scholar ] [16]. d’Elia N, Vanetti F, Cempini M, Pasquini G, Parri A, Rabuffetti M, Ferrarin M, Molino Lova R, and Vitiello N, “Physical human-robot interaction of an active pelvis orthosis: toward ergonomic assessment of wearable robots,” Journal of neuroengineering and rehabilitation, vol. 14, no. 1, p. 29, 2017. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [17]. Ronsse R, Vitiello N, Lenzi T, Van Den Kieboom J, Carrozza MC, and Ijspeert AJ, “Human–robot synchrony: flexible assistance using adaptive oscillators,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 4, pp. 1001–1012, 2010. [ DOI ] [ PubMed ] [ Google Scholar ] [18]. Fukuchi CA, Fukuchi RK, and Duarte M, “Effects of walking speed on gait biomechanics in healthy participants: a systematic review and meta-analysis,” Systematic reviews, vol. 8, no. 1, p. 153, 2019. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [19]. Kim G and Sergi F, “Using bayesian optimization to identify optimal exoskeleton parameters targeting propulsion mechanics: A simulation study,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2021, pp. 6225–6231. [ Google Scholar ] [20]. Gil-Agudo Á, Megía-García Á, Pons JL, Sinovas-Alonso I, Comino-Suárez N, Lozano-Berrio V, and Del-Ama AJ, “Exoskeleton-based training improves walking independence in incomplete spinal cord injury patients: results from a randomized controlled trial,” Journal of neuroengineering and rehabilitation, vol. 20, no. 1, p. 36, 2023. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] ACTIONS View on publisher site PDF (2.8 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top