Automatic Electromechanical Perturbator for Postural Control Analysis Based on Model Predictive Control

Abstract

Objective clinical analyses are required to evaluate balance control performance. To this outcome, it is relevant to study experimental protocols and to develop devices that can provide reliable information about the ability of a subject to maintain balance. Whereas most of the applications available in the literature and on the market involve shifting and tilting of the base of support, the system presented in this paper is based on the direct application of an impulsive (short-lasting) force by means of an electromechanical device (named automatic perturbator). The control of such stimulation is rather complex since it requires high dynamics and accuracy. Moreover, the occurrence of several non-linearities, mainly related to the human–machine interaction, signals the necessity for robust control in order to achieve the essential repeatability and reliability. A linear electric motor, in combination with Model Predictive Control, was used to develop an automatic perturbator prototype. A test bench, supported by model simulations, was developed to test the architecture of the perturbation device. The performance of the control logic has been optimized by iterative tuning of the controller parameters, and the resulting behavior of the automatic perturbator is presented.

1. Introduction

Postural control can be enlisted as one of the most critical functionalities of the human body. It is fundamental to perform everyday life activities, but also for managing random perturbations coming from the interaction with the environment. The central nervous system implements different postural strategies thanks to the information provided by the somatosensory, vestibular, and visual systems. The apparatus continuously commands the contraction and relaxation of the body’s muscles to correct the posture and reduce the risk of falling.

In order to get objective results about postural control performance, it is possible to perform balance analyses through dedicated instrumentation and protocols. Posturographic trials are performed in static (unperturbed) and dynamic (perturbed) conditions, depending on the scope of the analyses. In general, unpredictable mechanical perturbations are applied to the body of the patient to generate a loss of balance that is objectively evaluated and correlated to the features of the stimuli. Some devices directly impart a force by pushing/pulling the subject to/from a chosen point of application [1,2,3,4,5,6]. This approach can ease the configuration of the stimuli, however, it presents several complications regarding the control of the perturbation’s magnitude and duration. Other devices exert the perturbation through sliding/tilting of the base of support (BoS) [7,8,9] or by means of obstacles placed on the path that disturb the patient during walking [10,11], which translate into a disturbance to the body’s center of mass (CoM). This kind of stimulation, however, does not include the possibility to replicate a rather simple and typical balance perturbation as the one defined by a collision between the patient and an external body. Moreover, it does not allow to freely configure the characteristics of the stimuli, e.g., in terms of point of application and direction. The postural response to the stimuli is measured by center of pressure (CoP) related measurements [12,13,14,15], motion analysis [16,17,18,19], and/or electromyography [20,21,22].

A posturographic trial involves the interaction between a perturbation system (i.e., the perturbator) and the human being, i.e., the patient. Moreover, the perturbation system could be directly maneuvered by another human being, i.e., the operator. Therefore, it is particularly relevant to study the complex human–machine interfaces resulting from this clinical task. Mechatronic systems dealing with the human body are typically the object of research activities in a wide set of applications, as collaborative robots (cobots) environment [23,24] or robotics for rehabilitation [25,26] and surgery [27]. In many applications, it is critical to keep the interaction forces below specific thresholds or to control the actual evolution of forces rising due to human–machine collisions over time. For these reasons, active control solutions based on the real-time monitoring of physical quantities from the system (i.e., the plant) are often implemented. Given a certain application, control schemes with different complexity are designed in order to achieve reasonable performance in these tasks, which shall require the ability to track desired force profiles during sudden contacts, thus assuring the safety required. In previous works from the authors [28,29,30], a pneumo-tronic perturbator device for dynamic posturography has been developed, aimed at providing controlled stimuli with defined amplitude and duration to the body of a patient. Based on a linear pneumatic actuator, this device has only partially allowed fulfilling the application’s specifications, especially regarding the duration of the stimuli. Since it is advisable to end each stimulation before the actual balance control intervention, to reduce any overlap between active neuromuscular response and the force imparted by the device, it is particularly relevant to limit the duration approximately below 75 ms [31]. This specification can be hardly achieved with pneumatic actuated devices and different actuators (e.g., electric or electromechanical) have to be explored to optimize the performance of the system. Theoretically, using a linear electric motor provides several advantages over pneumatics: the nonlinearities due to airflow and friction of pneumatic actuators are no longer present, resulting in more robust control strategies. Moreover, due to the lack of compressed air supply, the system’s compactness and portability can be improved.

The implementation of control architectures in human–machine interface applications can rely either on linear controllers or on complex, non-linear controllers, depending on the complexity of the systems considered and on the level of performance to achieve. Examples of control architectures developed for similar applications are hybrid force/position [32,33] and impedance [34] or admittance [35] controllers. In addition, more complex control schemes have been developed to enhance the performance of traditional control schemes [36,37]. Although a more refined, articulated control logic can be particularly more effective with respect to simpler architectures, the ease of implementation and reliability of such structures have to be considered during the design process. With respect to pneumatics, electric actuation includes different aspects to be investigated, that may require the definition of more robust control logics, as Model Predictive Control (MPC) [38,39], Generalized Predictive Control [40], Sliding Mode Control [41,42], Adaptive Control [43]. In these control architectures, the quality of the results can be directly influenced by the accuracy of the plant model. Moreover, their implementation can also require a deep knowledge about the state of the system, which brings the necessity for additional transducers (if compared to simpler, linear controllers) to monitor the system’s variables. In some cases, when such physical quantities cannot be directly measured in a convenient way, software observers aimed at calculating the unknown variables on the basis of the measured ones can be developed. Among the different approaches available, MPC allows designing a controller that is able to operate the system with proficiency within a certain set of constraints (e.g., the maximum peak force achievable by the actuator under safety conditions) by implementing the ability to predict the system’s behavior. With respect to other approaches, this solution can help to avoid saturation problems of the control input, and simultaneously taking into account the variability of unknown parameters of the plant through proper state observers.

The aim of this paper is to investigate if electric actuation, in combination with MPC, can be used to develop an automatic perturbation device, able to impart accurate short-lasting forces (that we refer to as impulsive forces, with a limited duration of about 75 ms), to be used in dynamic posturographic trials to assess balance recovery skills. With respect to similar systems presented in the literature, which typically are bulky and provide a loss of balance through the BoS, this novel handheld device should enable customizable and repeatable perturbations by regulating the contact force raising from the impact between the human body and the perturbator itself. By accurately setting both the point of application and the magnitude of the perturbation, clinicians can adapt the characteristics of the postural disturb to the physical characteristics of the patient, define wide conditions of investigation, and explore the patterns of postural reactions in a comprehensive and clinically systematic way. The linearized model required to design the control logic is briefly described, as well as the test-bench used to test the control architecture. The model simulations and the experimental trials performed are discussed, with a particular focus on the accuracy and repeatability of the perturbations imparted.

2. Materials and Methods

2.1. Test-Bench and Automatic Perturbator Architecture

This study presents the development of an electromechanical perturbation device used in dynamic posturographic analysis. As shown in previous studies [28,29,30], the objective of the application is to analyze the importance of the force impulse (i.e., the time integral of the force signal over time, FI) in the subsequential balance recovery action: the postural response, expressed by the maximal displacement of the CoP, is proportional to the impulse rather than to the peak force of the perturbation [13]. The device focuses on a lightweight and compact design solution that can be easily handled by an operator. During the test, the operator maintains the perturbator’s interface at a distance of about 2 cm from the subject’s body and decides in an arbitrary way which body part to strike, as shown in Figure 1c. Throughout the experiment, the subject has to maintain the orthostatic position without losing contact with a BoS, which outputs, through proper load cell sensors measurements, the CoP displacement to be analyzed for the clinical assessment of postural control. The actuation is based on a linear electric motor, whose rod can be triggered by a push-button, which is responsible to impress, once the target body is reached, a rectangular contact force profile. Among the different types of linear motors, a tubular linear synchronous motor was chosen, since it allows accurate control of the rod’s motion while developing the high accelerations necessary to generate force up to 100 N in briefly time intervals, lasting about 75 ms or, in any case, less than 300 ms. By taking into account the expected initial distance between the perturbator and the patient, as well as the magnitude of his oscillation occurring during the perturbation, the stroke of the actuator was selected to be not lower than 100 mm. The development and test of the electromechanical perturbation device have been made by means of an appropriate Hardware-In-the-Loop (HIL) equipment which is shown in Figure 1a,b. A similar solution, based on pneumatic actuation, has been already presented in the literature [30]. The linear electric actuator (1) and the target body weight (8) are placed on two independent plates (7) able to slide on linear guides (10) and are connected to the fixed frame by means of two viscoelastic dampers (9), which model the passive reactions of the operator and of the subject, respectively. The displacement transducers (11) are used to measure the absolute positions of the electric actuator and mass. The laser sensor (3) is fixed on the top of the actuator’s barrel and is used to measure the stroke length. Finally, a load cell (4) is fixed to the rod end to measure the contact force between the actuator and the target body weight. The list of HIL key components is reported in

Table 1. List of test bench key components for HIL testing. Numbers refer to labels in Figure 1.

Number	Component
1	GD160Q motor (NiLAB GmbH, Klagenfurt am Wörthersee, Austria)
2	SLVD1N driver (Parker Hannifin Corp., Cleveland, USA)
3	Q4XTULAF400-Q8 optical sensor (Banner Engineering Corp., Plymouth, MN, USA)
4	UMM 50 kgf-ranged load cell (Dacell Co. Ltd., Cheongju, Korea)
5	Expanded polyethylene interface
6	Baseline Real-Time target machine (Speedgoat Inc., Natick, MA, USA)
7	Aluminum sliding plates
8	Target body weight
9	Custom made viscoelastic dampers
10	C-SHR28-1000-B4 linear guides (MISUMI, Europa GmbH, Frankfurt, Germany)
11	PZ-34-A-100 displacement transducers (GEFRAN, Provaglio d’Iseo, Italy)
12	DEWE-RACK-4 (Dewetron GmbH, Grambach, Austria)
13	Connection box for Real-Time target system

.

2.2. Model and Control System Design

The system’s control logic was developed by the finite state machine (FSM) criterium as in the previous prototypes [28,29,30]. The FSM involved the following phases:

• Idle: the actuator’s rod remains still in the retracted position;
• Approach: the piston nears the target body following a constant speed reference signal. It is triggered by the occurrence of the start command delivered by the operator;
• Strike: the actuator meets the target body and force reference control is issued. It is triggered by a 2.5 N threshold on the load cell signal;
• Retraction: the actuator’s rod returns to the retracted position with constant speed after having completed or failed to enact the strike phase.

The rod’s idle, approach, and retraction phases, concerning only the actuator part, are managed by the driver’s embedded PI speed controller. The strike phase, of higher importance due to the phenomena involving the human–machine interaction, is controlled by the real-time machine through an MPC-based force control design. The driver’s embedded PLC is responsible for the toggling between the two types of control.

System modeling has thus exclusively involved the strike phase. The human body recovers low balance disturbances in the orthostatic position configuration mainly by applying a correcting torque at the ankle (Tc). Dynamic equilibrium is reached once the body’s CoM is vertically aligned with the CoP, i.e., the resultant of the gravity forces ($-m_{tot}g$) lies on the same (vertical) axis of the ground reaction force (${\mathrm{R}}_{\mathrm{y}}$). This approximation allows studying human balance control in the sagittal plane through a simple inverse pendulum model, as shown in Figure 2a. The same approach can be used to describe the oscillations of the operator about the ankle joint while maneuvering the automatic perturbator.

As for the nonlinearities involved in the strike phenomenon, the literature offers many solutions in order to efficiently model the contact force, mainly depending on the properties of the materials involved and on the definition of boundary conditions. Nonlinearities generally refer to polynomial forms in the bodies’ relative penetration $\delta$, its respective speed $\dot{\delta }$, and a hysteresis damping factor $\chi$ [44]. Due to the difficult measurement of these parameters, said aforementioned nonlinear relations were approximated by adopting a lumped parameter system that models the actuator’s and target’s bodies as directly connected, therefore modeling the contact force as a stress. This approximation, in addition to linearizing the angular configuration of the phenomenon, as in Figure 2a, through Taylor’s series, leads to the linearized lumped parameter model shown in Figure 2b. The system includes the appropriate transducers required to monitor relevant state variables, as shown in the test-bench design (see Figure 1). The modeled system neglects the CoP displacement signal, which is the target of biomedical analysis only and therefore is not used in the control loop computations. All the system parameters are listed in

Table 2. Modeled system key parameters list. Symbols refer to Figure 2b. All displacements are referred to the respective equilibrium positions.

Symbols	Description
${\mathrm{m}}_{\mathrm{a}},{\mathrm{x}}_{\mathrm{a}}\left(\mathrm{t}\right),{\mathrm{F}}_{\mathrm{a}}\left(\mathrm{t}\right)$	Actuator’s mass, displacement, and issued force
${\mathrm{c}}_{\mathrm{a}},{\mathrm{k}}_{\mathrm{a}}$	Operator response’s damping and stiffness
${\mathrm{m}}_{\mathrm{b}},{\mathrm{c}}_{\mathrm{b}},{\mathrm{k}}_{\mathrm{b}},{\mathrm{x}}_{\mathrm{b}}\left(\mathrm{t}\right)$	Target body’s mass, linear damping, stiffness, and displacement
${\mathrm{m}}_{\mathrm{r}},{\mathrm{x}}_0$	Rod’s mass, elongation at impact
${\mathrm{c}}_{\mathrm{i}},{\mathrm{k}}_{\mathrm{i}},{\mathrm{F}}_{\mathrm{i}}\left(\mathrm{t}\right),\mathrm{x}\left(\mathrm{t}\right)$	Interface’s damping, stiffness, contact force, and displacement

.

The Lagrange approach, adopting the displacements displayed in Table 2 as degrees of freedom q(t), leads to the following dynamic equations:

$$\left[\mathrm{M}\right]\stackrel{\cdot \cdot }{\mathrm{q}}\left(\mathrm{t}\right)+\left[\mathrm{C}\right]\stackrel{\cdot }{\mathrm{q}}\left(\mathrm{t}\right)+\left[\mathrm{K}\right]\mathrm{q}\left(\mathrm{t}\right)=\left\{\mathrm{F}\left(\mathrm{t}\right)\right\}, \mathrm{q}\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\mathrm{x}}_{\mathrm{a}}\left(\mathrm{t}\right)\\ \mathrm{x}\left(\mathrm{t}\right)\\ {\mathrm{x}}_{\mathrm{b}}\left(\mathrm{t}\right)\end{array}\right\}$$

(1)

$$\left[\mathrm{M}\right]=\left[\begin{array}{ccc}{\mathrm{m}}_{\mathrm{a}}& 0& 0\\ 0& {\mathrm{m}}_{\mathrm{r}}& 0\\ 0& 0& {\mathrm{m}}_{\mathrm{b}}\end{array}\right], \left[\mathrm{C}\right]=\left[\begin{array}{ccc}{\mathrm{c}}_{\mathrm{a}}& 0& 0\\ 0& {\mathrm{c}}_{\mathrm{i}}& -{\mathrm{c}}_{\mathrm{i}}\\ 0& -{\mathrm{c}}_{\mathrm{i}}& {\mathrm{c}}_{\mathrm{i}}+{\mathrm{c}}_{\mathrm{b}}\end{array}\right], \left[\mathrm{K}\right]=\left[\begin{array}{ccc}{\mathrm{k}}_{\mathrm{a}}& 0& 0\\ 0& {\mathrm{k}}_{\mathrm{i}}& -{\mathrm{k}}_{\mathrm{i}}\\ 0& -{\mathrm{k}}_{\mathrm{i}}& {\mathrm{k}}_{\mathrm{i}}+{\mathrm{k}}_{\mathrm{b}}\end{array}\right]$$

(2)

$$\left\{\mathrm{F}\left(\mathrm{t}\right)\right\}=\left\{\begin{array}{c}-{\mathrm{F}}_{\mathrm{a}}\left(\mathrm{t}\right)\\ {\mathrm{F}}_{\mathrm{a}}\left(\mathrm{t}\right)\\ 0\end{array}\right\}$$

(3)

Standard dynamic to state-space (SS) representation transformation is performed through the reduction of order technique by the following equations:

$$\{\begin{array}{c}\left[\mathrm{M}\right]\stackrel{\cdot }{\mathrm{v}}+\left[\mathrm{C}\right]\mathrm{v}+\left[\mathrm{K}\right]\mathrm{q}=\left\{\mathrm{F}\left(\mathrm{t}\right)\right\}\stackrel{\mathrm{def}}{=}\left\{\mathrm{T}\right\}\mathrm{u}\left(\mathrm{t}\right)\\ \stackrel{\cdot }{\mathrm{q}}\left(\mathrm{t}\right)=\mathrm{v}\left(\mathrm{t}\right)\end{array},\mathrm{z}\left(\mathrm{t}\right)=\left\{\begin{array}{c}\mathrm{v}\left(\mathrm{t}\right)\\ \mathrm{q}\left(\mathrm{t}\right)\end{array}\right\},\mathrm{T}=\left\{\begin{array}{c}-1\\ 1\\ 0\end{array}\right\}$$

(4)

$$\stackrel{\cdot }{\mathrm{z}}\left(\mathrm{t}\right)=\left[\begin{array}{cc}-{\left[\mathrm{M}\right]}^{-1}\left[\mathrm{C}\right]& -{\left[\mathrm{M}\right]}^{-1}\left[\mathrm{K}\right]\\ {\mathrm{I}}_{3\times 3}& {\mathrm{O}}_{3\times 3}\end{array}\right]\mathrm{z}\left(\mathrm{t}\right)+\left[\begin{array}{c}{\left[\mathrm{M}\right]}^{-1}\left\{\mathrm{T}\right\}\\ 0_{3\times 1}\end{array}\right]\mathrm{u}\left(\mathrm{t}\right)\stackrel{\mathrm{def}}{=}\left[\mathrm{A}\right]\mathrm{z}\left(\mathrm{t}\right)+\left[\mathrm{B}\right]\mathrm{u}\left(\mathrm{t}\right)$$

(5)

$$\mathrm{y}\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\mathrm{F}}_{\mathrm{i}}\left(\mathrm{t}\right)\\ \mathrm{x}\left(\mathrm{t}\right)-{\mathrm{x}}_{\mathrm{a}}\left(\mathrm{t}\right)\end{array}\right\}=\left[\begin{array}{ccc}0& -{\mathrm{c}}_{\mathrm{i}}& {\mathrm{c}}_{\mathrm{i}}\\ 0& 0& 0\end{array}\begin{array}{ccc}0& -{\mathrm{k}}_{\mathrm{i}}& {\mathrm{k}}_{\mathrm{i}}\\ -1& 1& 0\end{array}\right]\mathrm{z}\left(\mathrm{t}\right)+\left[0_{2\times 1}\right]\mathrm{u}\left(\mathrm{t}\right)\stackrel{\mathrm{def}}{=}\left[\mathrm{C}\right]\mathrm{z}\left(\mathrm{t}\right)+\left[\mathrm{D}\right]\mathrm{u}\left(\mathrm{t}\right)$$

(6)

in which z(t) is the state of the system, u(t) the control input, T the selection vector, and [A], [B], [C], [D] the SS representation quadruplet. Since the laser sensor could be removed from the system in an eventual next prototype, the adopted control design uses only the first row of [C], using the laser sensor only for limit switch purposes. The initial state array is:

$${\mathrm{z}}_0\stackrel{\mathrm{def}}{=}\mathrm{z}\left(\mathrm{t}=0\right)={\left\{{\stackrel{\cdot }{\mathrm{x}}}_{\mathrm{a}}\left(0\right),\stackrel{\cdot }{\mathrm{x}}\left(0\right),{\stackrel{\cdot }{\mathrm{x}}}_{\mathrm{b}}\left(0\right),{\mathrm{x}}_{\mathrm{a}}\left(0\right),\mathrm{x}\left(0\right),{\mathrm{x}}_{\mathrm{b}}\left(0\right)\right\}}^t\cong {\left\{0,{\mathrm{v}}_0,0,0,0,0\right\}}^t$$

(7)

in which ${\mathrm{v}}_0$ is the approach speed. Such SS system is then reduced to its minimal realization through Mathworks Control System Toolbox™, thus yielding positive results in both controllability and observability tests, as both the controllability and the observability matrices are fully ranked.

The modeled system’s SS representation is then inputted in the MPC, which was adopted due to high constraint adaptiveness. In fact, the peak force of the selected motor (105 N) is too close to the highest force level to be expected (about 100 N) for choosing serendipitous, i.e., simpler and software-saturated, control strategies such as PID control. The MPC bases its computations on the following control subjects: prediction horizon, quadratic program (QP) optimization, and receding horizon (RH) technique. Said control algorithm was adopted due to constraint adaptiveness.

Prediction horizon (Hp) is the number of steps ahead to be computed (and controlled) by the algorithm, forming the macro vector:

$$\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)=\left\{\mathrm{u}\left(\mathrm{k}|\mathrm{k}\right)\mathrm{u}\left(\mathrm{k}+1|\mathrm{k}\right)\dots \mathrm{u}\left(\mathrm{k}+{\mathrm{H}}_{\mathrm{p}}-1|\mathrm{k}\right){\mathsf{\epsilon }}_{\mathrm{k}}\right\}$$

(8)

in which caption “i|j” means that the prediction of the i-th time step is being calculated during the j-th time step and ${\mathsf{\epsilon }}_{\mathrm{k}}$ is the slack variable for soft constraint violation;

QP vector U(k|k) optimization is subject to plant system morphology (SS quadruplet) and constraints (in this case, the actuator’s peak force), through the following cost function:

$$\mathrm{J}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)={\mathrm{J}}_{\mathrm{y}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)+{\mathrm{J}}_{\mathrm{u}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)+{\mathrm{J}}_{\Delta \mathrm{u}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)+{\mathrm{J}}_{\mathsf{\epsilon }}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)\stackrel{\mathrm{def}}{=}{\displaystyle \sum }_{\mathsf{\xi }=1}^4{\mathrm{J}}_{\xi }\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)$$

(9)

$${\mathrm{J}}_{\mathsf{\xi }}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)={\displaystyle \sum }_{i=0}^{H_p-1}{\xi }^t\left(\mathrm{k}+i|\mathrm{k}\right)\mathsf{\Omega }\xi (\mathrm{k}+i|\mathrm{k})$$

(10)

As Equation (10) shows, each cost function part is regulated by a tunable weight matrix $\mathsf{\Omega }$. The higher the weight, the more optimized the respective target $\xi$ in the QP optimization. A description of every addendum of the cost function is shown in

Table 3. QP cost function parameters legend.

Symbol	$\mathbf{Optimization}\mathbf{Target}\mathbf{\xi }$	$\mathbf{Tunable}\mathbf{Weight}\mathbf{\Omega }$
${\mathrm{J}}_{\mathrm{y}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)$	Tracking error	Q
${\mathrm{J}}_{\mathrm{u}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)$	Control input	${\mathrm{R}}_{\mathrm{u}}$
${\mathrm{J}}_{\Delta \mathrm{u}}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)$	Control input rate	${\mathrm{R}}_{\Delta \mathrm{u}}$
${\mathrm{J}}_{\mathsf{\epsilon }}\left(\mathrm{U}\left(\mathrm{k}|\mathrm{k}\right)\right)$	Soft constraint violation	${\mathsf{\rho }}_{\mathsf{\epsilon }}$ 1

¹ In this case equal to zero due to the presence of hard constraints only.

.

RH technique extracts control input u(k|k) from the QP optimized macro vector U(k|k) and inputs it to the plant.

Since the state of the system is not totally measured, the controller needs the implementation of a state observer. The implementation used a Kalman Filter observer, whose matrices L and M are set by default as optimal static gains for the plant model by the MathWorks Model Predictive Control Toolbox™ and Control System Toolbox™. The Kalman Filter asymptotic stability is thus checked positive through the following eigenvalues argument investigation:

$$\max (\mathrm{arguments}\left(\mathrm{eigenvalues}\left(\mathrm{A}-\mathrm{LC}\right)\right))<1$$

(11)

A comprehensive scheme showing the MPC algorithm loop is portrayed in Figure 3.

The resulting control input signal, which was originally computed as a force (N), had to be converted into a voltage signal (V) in order to be correctly processed by the driver of the electric motor. The whole system was modeled through the use of Simulink^® (The MathWorks Inc., Natick, MA, USA) environment.

2.3. Simulations and Experimentation

To investigate the behavior of the control system, preliminary analyses were performed by Model-In-the-Loop (MIL) testing of the strike phase. Then, the control system was implemented in the test bench shown in Figure 1a,b for HIL testing. The latter revolved around the trial and error optimization of Hp and of the QP tunable weights shown in Table 3, in order to minimize the tracking error and to improve the constrained control input behavior as much as possible. To have the most fitting profiles, the control input gain performing the V/N conversion, found through static measurements and equal to 0.125 V/N, was also modified by tuning a dimensionless corrective factor (CF).

Before starting the HIL test, it was necessary to electrically isolate the load cell from the motor to block any parasitic current formed by the alimented phases of the motor. In addition, the load cell signal was filtered with a spike recognition and removal digital filter. The system also implemented a filter (FN 2010-20-06, Schaffner EMC Srl, Monza, Italy) at the driver’s power supply to limit noise affecting the instrumentation.

To evaluate the postural response to the stimulus, the perturbation must be large enough to elicit a clear detectable postural response but not evoking a step strategy response. Moreover, the test should not be painful for the patient. Based on preliminary experiments, the perturbation’s amplitude and time interval (TI), namely and respectively, should be set between 20 N and 100 N and between 75 ms and 250 ms. Both MIL and HIL tasks focused on two rectangular reference signals:

• R1: rectangular reference signal of 50 N with 75 ms duration (reference FI = 3.75 Ns);
• R2: rectangular reference signal of 50 N with 250 ms duration (reference FI = 12.5 Ns).

For each trial, the TI of each perturbation was defined by the intersections of the force signal with a defined threshold (equal to 3.5 N), while the FI was calculated through the trapezoidal rule. Performance indices were the TI and FI deviations, respectively, TID and FID. Each deviation was computed as:

$$\mathrm{Deviation}=\left(1-\frac{\mathrm{measured}\mathrm{value}}{\mathrm{reference}\mathrm{value}}\right)\ast 100$$

(12)

For each set of control parameters and test condition, five equal perturbations were delivered during the HIL test. The coefficient of variation (CV) of the FI, i.e., the standard deviation divided by the absolute mean value, was estimated to measure the repeatability of the perturbations. The CV was multiplied by 100 to ease the readability of the results.

3. Results

3.1. MIL Results

MIL MPC performance results are presented in Figure 4 and

Table 4. MIL MPC tuned parameters and the resulting performance indices. The labels and colors refer to the corresponding curves presented in Figure 4.

Profile	${\mathbf{H}}_{\mathbf{p}}\left(\mathbf{Steps}\right)$	Q	${\mathbf{R}}_{\mathbf{u}}$	${\mathbf{R}}_{\mathbf{\Delta }\mathbf{u}}$	Ref	TI (s)	TID (%)	FI (Ns)	FID (%)
P1	10	3.8	0	10	R1	0.264	−5.6	12.56	−0.47
(cyan)					R2	0.089	−18.67	3.82	−1.77
P2	10	6	0	10	R1	0.263	−5.2	12.57	−0.56
(brown)					R2	0.088	−17.33	3.83	−1.98
P3	10	3.8	1	10	R1	0.269	−5.6	11.00	11.96
(pink)					R2	0.089	−18.67	3.37	10.13

. As shown in Table 4, the MPC algorithm gives a good performance in terms of TID, FID, and FI CV.

3.2. HIL Results

HIL MPC performance results are shown in Figure 5 and Figure 6. The performance indices obtained for the HIL are reported in

Table 5. HIL MPC tuned parameters and the resulting performance indices. The colors refer to the corresponding curves presented in Figure 5 and Figure 6.

Profile	${\mathbf{H}}_{\mathbf{p}}\left(\mathbf{Steps}\right)$	Q	${\mathbf{R}}_{\mathbf{u}}$	${\mathbf{R}}_{\mathbf{\Delta }\mathbf{u}}$	CF	Ref	TI (s)	TID (%)	FI (Ns)	FID (%)	FI CV * 100
P1	10	3.8	0	10	1.14	R1	0.256	−2.4	8.35	33.17	7.14
(cyan)					1.45	R2	0.077	−2.7	3.08	17.92	2.66
P4	5	3.8	0	10	1	R1	0.252	−0.8	10.59	15.29	0.46
(black)					1.45	R2	0.077	−2.7	3.83	−2.19	0.50
P5	5	1.5	0	10	1.23	R1	0.254	−1.6	9.36	25.10	0.34
(orange)					2.29	R2	0.078	−4	2.93	21.87	0.63
P6	5	6	0	10	1	R1	0.255	−2	11.01	11.93	0.88
(violet)					1	R2	0.077	−2.7	3.03	19.12	1.47
P7	5	3.8	1	10	1	R1	0.252	−0.8	6.69	46.47	0.55
(blue)					1.45	R2	0.076	−1.3	2.72	27.35	1.25

. As highlighted by comparison between Figure 4 and Figure 5, the MIL optimized tunes (P1 profile) show a significant performance loss if applied to the real system. As reported in Table 5, the optimized black profile (P4) shows the best compromise results in terms of the analyzed output signal performances.

4. Discussion

As shown in Figure 4, the results of the simulations highlighted that the proposed architecture is able to provide sufficiently accurate impulse (FI) values for different selections of the control parameters. In particular, the model showed low deviation values of force impulse and duration in two of the configurations studied, i.e., the ones corresponding to P1 and P2 profiles. For these profiles, the larger error observed was related to the duration of the force profiles in the R2 condition, when an impulsive force lasting only 75 ms had to be tracked by the control system. As expected, this condition is particularly demanding since it requires high dynamics and accuracy. However, the TI observed (about 89 ms) is still appropriate for the application. On the other hand, the R1 reference profile, lasting 250 ms, was more efficiently tracked.

Similar behavior was observed in HIL testing, although the tracking accuracy, especially for the R2 profile, was definitely lower than the one observed in MIL simulations. The TI accuracy observed was particularly high, with low deviations for all the selections of control parameters (see Table 5). This result confirms that the electric actuated solution presented in this work can achieve particularly fast responses, which are fundamental to reduce the overlap between the perturbation and the active neuromuscular response of the subject. This condition is necessary to avoid biased analyses of postural control, and it could not be achieved by a former, pneumatic actuated solution presented in previous research [30]. Furthermore, as presented in Figure 6, the HIL testing showed high repeatability of the stimuli, which was confirmed for most of the configurations considered as highlighted by the low FI CV values reported in Table 5. This result is relevant to ease the definition of clinical protocols that, in addition, must guarantee the safety of the experimentation thus ensuring to obtain a measure of the subject’s postural performance comparable over time. As shown in Table 5, better performances were achieved by shortening Hp, i.e., the prediction power of the controller. This is likely because the controller bases its computation on a highly simplified model, whose linearity could serve as a sufficient approximation of the reality only for limited future time intervals.

The mismatch observed between MIL and HIL results, especially for the R2 force profile, can be also motivated by the approximations considered for the modeling of the plant. In particular, as discussed in Section 2.2, the nonlinearities involved in the impact phase were generally left out from the definition of the model for the sake of simplicity. This strong approximation was necessary to evaluate the feasibility of the study, and to test the efficacy of the MPC algorithm, in such a complex scenario where several parameters are unknown and cannot be easily measured. A more accurate model representation of the actuation system, including, for example, the inertia of the piston and the driver response and update times, could also improve the efficiency of the control system as it takes into account the latency between the start of the reference signal and the start of the perturbation. Among the several aspects to be considered, the finite response time of the actuation system has likely the most significant impact on the MPC strategy, since the controller assumes the control input to be instantly issued to the plant and, thus, acts accordingly. Moreover, in a real human–machine interaction, it could be even more difficult to retrieve unknown data related to the physical characteristics of the subject, which could be highly specific and impractical to measure. However, the outcomes of this study signal that, in order to achieve a better matching between the MIL and HIL results, a refinement of the analytical model is necessary to improve the performance of the predictive algorithm. For example, detailed modeling of the actuator’s dynamics could be necessary to overcome part of the limitations observed in the current architecture of the control system, especially for short-lasting perturbations. In addition, a more refined model of the plant could help in the identification of the optimal control parameters, that have been obtained by iterative trial and error approach in the current study.

With respect to the model simulations, experimental trials showed lower tracking accuracy, as confirmed by the significantly higher force impulse deviations collected in Table 5. Among the several selections of control parameters, the one corresponding to the P4 profile (see Figure 5 and Figure 6) showed the best results in terms of impulse accuracy. In such configuration, the FID observed for the short-lasting R2 reference profile was even lower than the one measured for the R1 profile, however, it was not possible to assess the same result regarding the tracking accuracy. As shown in Figure 6, the oscillations in the force profile were particularly substantial for the 75 ms reference profile. A similar result was highlighted for the other selections of control parameters (see Figure 5 and Table 5). A large overshoot in the contact force profile could be potentially dangerous for the subject, therefore it must be avoided in a final embodiment of the perturbation system. The CF was tuned to achieve a reasonable compromise between the accuracy of the force impulse and the steadiness of the control but, still, it was not possible to remove the undesired overshoot without affecting the general performance of the perturbation device.

5. Conclusions

A linear electric motor, in combination with Model Predictive Control, was used to develop an automatic perturbator prototype for postural control analysis. The application requires that high forces are developed in a short time interval and that the perturbation is accurately controlled to avoid patient’s injuries. Whereas traditional systems are based on the shifting and tilting of the BoS below the patient, the device presented in this work is aimed to exert controlled stimuli to the human body, with no limitations in terms of magnitude and point of application or direction of the perturbation.

While both the simulation and the HIL experimental trials have shown that the dynamics of the system is appropriate for the application, some limitations were found in the accuracy of the impulse and of force profile tracking, likely because the controller formulation is based on a highly simplified model of the real experimentation setting.

Future steps will focus on improving the accuracy of the plant model, in combination with a more refined actuator’s transfer function system identification, to enhance the predictiveness of the MPC algorithm and thus the control performances. The plant’s model improvement will also involve a proper modeling of the involved strike nonlinearities, in parallel with exploring the option of a nonlinear MPC action.

6. Patents

Ferraresi, C., Franco, W., Maffiodo, D., De Benedictis, C., Roatta, S., Dvir, Z. Striker unit for postural analysis, Italian patent n. 102018000010030, filed 5 November 2018 and issued 8 October 2020.