%=============================================================================== % $Id: ifacconf.tex 19 2011-10-27 09:32:13Z jpuente $ % Template for IFAC meeting papers % Copyright (c) 2007-2008 International Federation of Automatic Control %=============================================================================== \documentclass{ifacconf} % Add packages and commands here. The following packages are loaded in our class file: fontenc, calc, indentfirst, fancyhdr, graphicx, lastpage, ifthen, lineno, float, amsmath, setspace, enumitem, mathpazo, booktabs, titlesec, etoolbox, amsthm, hyphenat, natbib, hyperref, footmisc, geometry, caption, url, mdframed, tabto, soul, multirow, microtype, tikz \usepackage{graphics} % for pdf, bitmapped graphics files %\usepackage{epsfig} % for postscript graphics files \usepackage{mathptmx} % assumes new font selection scheme installed %\usepackage{times} % assumes new font selection scheme installed \usepackage{amsmath} % assumes amsmath package installed \usepackage{amssymb} % assumes amsmath package installed \usepackage[pdftex]{graphicx} \graphicspath{{figure/}} \DeclareGraphicsExtensions{.pdf,.jpeg,.png,.eps} %=============================================================================== \begin{document} \begin{frontmatter} \title{Identification and Robust Control of Heart Rate During Treadmill Exercise at Large Speed Ranges} % Title, preferably not more than 10 words. \author{Ali Esmaeili $^ {(1) }$, Asier Ibeas$^{ (1,2) }$, Jorge Herrera$^{(2)}$ and Nazila Esmaeili$^{ (3) }$} \address[First]{\quad Department of Telecommunications and Systems Engineering \\ Universitat Aut\`onoma de Barcelona, 08193 Bellaterra, Barcelona, Spain. Email: Ali.Esmaeili@e-campus.uab.cat \\} \address[Second]{\quad Dep. de Ingenier\'ia, Facultad de Ciencias Naturales e Ingenier\'ia ,\\ Universidad de Bogot\'a Jorge Tadeo Lozano, 110311, Bogot\'a DC, Colombia\\} \address[Third]{Faculty of Electrical Engineering and Information Technology, Otto-von-Guericke-Universitat, 39106, Magdeburg, Germany} %\thanks[{This work was partially supported by the Spanish Ministry of Economy and Competitiveness through grant DPI2016-77271-R and by the University of the Basque Country (UPV/EHU) through grant PPG17/33.} \begin{abstract} % Abstract of not more than 250 words. The objective of this paper is to design a heart rate (HR) controller for a treadmill so that the HR of an individual running on it tracks a pre-specified, potentially time-varying profile specified by doctors for the cardiac recovery of the person. Initially, a parameter estimation algorithm is presented with the aim of estimating the values of the parameters of a model relating the speed of the treadmill with the HR of an individual. The parameter estimation problem is formulated as an optimization one and solved by using Particle Swarm Optimization (PSO). Afterwards, a super-twisting sliding mode controller is designed to perform the robust control of treadmill's speed in the presence of potential unmodelled dynamics or parametric uncertainties. Numerical examples show that the estimation procedure is able to obtain accurate values for the system's parameters while the proposed control approach is able to obtain zero tracking error without chattering, definitely achieving the control objectives. In both cases the range of treadmill's speed goes from 2 to 14 km/h, range that is not usually employed in previous studies. \end{abstract} \begin{keyword} {Heart rate control, PSO Identification of nonlinear systems, Super-twisting sliding mode control} \end{keyword} \end{frontmatter} %=============================================================================== \section{Introduction} The objective of this paper is to design a heart rate (HR) controller for a treadmill so that the HR of an individual running on it tracks a pre-specified, potentially time-varying profile specified by doctors for the cardiac recovery of the person. The controller modifies the speed of the treadmill so as to make individual's HR track such a profile. The ability to control the Heart Rate (HR) in treadmill exercise is of the great importance in design of exercise routines and it is one of the most significant vital signs to reflect cardiovascular events, (Hunt et al, 2016). The understanding of HR response with exercise may also lead to an improvement in developing training protocols for athletics, more efficient weight loss protocols for the overweight people, and in facilitating assessment of physical fitness and health of individuals, (Weippert et al, 2014). In addition, treadmill exercise has an important role for people recovering from cardiac disease or surgery as well as for people involved in weight loss programs. The design of the controller follows two steps. The first one is the modeling of the HR response to the treadmill exercise while the second one designs the controller based on the previous model by means of a sliding mode control. Several models have been proposed in the literature to model heart rate response to treadmill velocity, (Jang and Dae-Geun, 2016), (Cheng et al, 2008), (Peter and Huber, 1964). For instance, (Shtessel and Yuri, 2010) proposes a Hammerstein model composed of a static non-linearity defined by a look-up table followed by a linear dynamical system, while (Zhang et al, 2011) and (Scalzi et al, 2012) propose a nonlinear dynamical model. Anyway, nonlinearity must be present in the model due to the nonlinear response of the heart rate to the exercise. In this study, we use the nonlinear dynamical system proposed by, (Zhang et al, 2011), which is able to capture the dynamic behavior in a compact way by means of a reduced number of parameters, fact that is convenient for control purposes. Heart rate treadmill models are parametrized by a number of parameters that capture the individual HR response to exercise. It is vital, thus, to design parameter estimation procedures that allow us to have a personalized model for each individual from measured data. In this work, parameter estimation is set up as an optimization problem whose solution leads to the estimated model parameters. The optimization problem is solved by using a Particle Swarm Optimization (PSO) algorithm. Particle swarm optimization was introduced by Kennedy and Eberhart in 1995, (Yan et al, 2013), (Shtessel and Yuri, 2010). At that time, the wide success of Evolutionary Algorithms (EAs) motivated researchers worldwide to develop and experiment with novel nature-inspired methods. Thus, besides the interest in evolutionary procedures that governed EAs, new paradigms from nature were subjected to investigation. The first PSO models introduced the novelty of using difference vectors among population members to sample new points in the search space. This novelty diverged from the established procedures of EAs, which were mostly based on sampling new points from explicit probability distributions. Additional advantages of PSO were its potential for easy adaptation of operators and procedures to match the specific requirements of a given problem, as well as its inherent decentralized structure that promoted parallelization, (Storn and Price, 1997). PSO has gained much attention nowadays and has wide applications in different fields such a fitness distance ratio, (Yan et al, 2013), adaptive mutation and inertia weight, (Zhang et al, 2011) and parameter identification in magnet synchronous motors, (Liu et al, 2008), hybrid neural network and the level of seismic inversion, (Yang et al,2017). In this study, the parameter estimation problem is solved by using the Particle Swarm Optimization (PSO) method, which is the first time where this method is used for the problem at hand. The parameter estimation will be focused in the range of treadmill speeds from 2 up to 14 km/h (2-14 km/h). In many studies, (Bansal et al, 2011), (Ibeas et al, 2016), (Shtessel and Yuri, 2010), the usual range of speed is (2-8 km/h), or (2-10 km/h). Therefore, the proposed methodology is applicable in the large range of speeds from 2 to 14 km/h. Despite this range is popular in many rehabilitation and training exercises, it is the first time considered in this problem. On the other hand, a Sliding Mode Controller (SMC) approach is adopted to design the controller. The SMC has revealed very useful in the robust control of multiple systems, such as pneumatic cylinder as actuators for robot manipulators, (Paul et al, 1994) and the hydraulic dynamics of the manipulator, (Guo et al, 2008).However, this is the first time that SMC is used to control the HR during treadmill exercise since previous works used different control techniques such $H_{\infty}$ robust control approaches or model predictive controllers, (Scalzi et al, 2012), (Shtessel and Yuri, 2017), (GAO and Xuehui, 2016). Especial attention will be devoted to the chattering effect since this is a crucial aspect in biomedical applications, (Lu et al, 2016). To this end, a super-twisting based sliding mode controller will be designed instead of a traditional SMC, (Shtessel and Yuri, 2010), (GAO and Xuehui, 2016). The super-twisting approach will allow avoiding the undesired oscillations that a classical sliding mode controller may cause. Simulation results will show that our approach definitely improves the accuracy of the model parameter estimation for these treadmill speeds, (up to 14 km/h) and the SMC also improves the heart rate control which has not been considered in the past. The rest of paper is organized as follows. In section (2) we introduce the problem formulation. In this part, we provide the model description. In the next section, the PSO algorithm proposed for the parameter estimation is introduced, while the advantages of this method are commented. In addition, Section (3) also contains the controller design procedure based on SMC. Finally, the last section presents the simulation examples and its comparison with alternative estimation and control methods. \section{Problem Formulation} The following nonlinear model describes the relationship between speed and heart rate during treadmill exercise, (Jang and Dae-Geun, 2016), (Su et al, 2010): \begin{eqnarray} \dot x_1(t) &=& - a_1 x_1(t) + a_2 x_2(t) + a_3 u^2(t)\\ \dot x_2(t) &=& - a_4 x_2 (t) +\phi(x_1(t))\\ \phi (x_1(t)) &=& \frac{{ a }_{ 5 }{ x }_{ 1 }(t) }{1 + \exp{(-( x_1(t) - a_6))}}\\ y(t) &=& x_1(t) + HR_{rest} \label{eq.mode} \end{eqnarray} Where $x(0) = [x_1(t) , x_2(t)] = [0 , 0]$ is the usual initial condition and ${a_1}$, ..., ${a_6}$ are positive scalars that are adjusted from real data to describe the particular response of each individual to exercise. The output $y(t)$ relates to the change of HR of the person, and $HR_{rest}$ is the value of the Heart Rate at rest. The control input $ u(t)$ describes the speed of the treadmill. The component $x_1(t)$ describes the change of HR from the heart rate at rest mainly due to the central response to exercise, whereas the component $ x_2(t)$ describes the slower and more complex local peripheral effects. The positive feedback signal $x_2$, or a dynamic disturbance input to the $x_1$ subsystem, may be treated as a reaction of HR to the effects from the peripheral local responses or factors. In this case, the metabolites from the peripheral local metabolism further accelerate the HR during exercise. For instance, in the case of the peripheral local metabolism, the accumulated metabolic by-products, such as adenosine, K+, H+, lactic acid and other metabolites, cause vasodilatation and hyperemia inactive muscles, (Su SW et al, 2010). Vasodilatation in the active muscles causes a reduction in total peripheral resistance which in turn causes a decrease in mean arterial blood pressure. In order to regulate the blood pressure, the cardiac output needs to be increased, meaning that stroke volume and HR are increased via the baroreceptor reflex, (Zhang et al, 2011). The nonnegative nonlinear function $\phi(x_1 )$ has the property that $\phi(x_1) << 1$ when $x_1$ is small, whereas when $x_1$ is much larger than $ a_6$, $\phi(x_1(t))$ approaches the linear function $x_1(t)$. If $x_1$ is small and $ a_6$ is large, the variable $x_1$ is multiplied by a small factor(i.e. $\frac { { a }_{ 4 } }{ 1+exp-{ x }_{ 1 }(t)-{ a }_{ 6 } } \approx 0$ in the second equation of (1), so $x_2$ becomes nearly independent of $x_1$ . If $x_1(0) = x_2(0) = 0$ and the input $u(t)$ is small, the state $x_1(t)$ may not be large enough to make the factor $ \frac { { a }_{ 4 } }{ 1+exp-{ x }_{ 1 }(t)-{ a }_{ 6 }}$ significant, and $x_2(t)$ will remain close to zero. As a result, system (1) can be approximated by the system $x_1(t) = - a_1x_1(t) + a_3u^2(t)$ with $x_2(t) = 0$. On the other hand, if the input $u(t)$ is sufficiently large, the state $x_1(t)$ will be driven to a level that the factor $\frac { { a }_{ 4 } }{ 1+exp-{ x }_{ 1 }(t)-{ a }_{ 6 }$ is significant, and $x_2(t)$ is no longer independent of $x_1(t)$. It is important to bear in mind that the objective of the paper is twofold. On the one hand to propose a parameter estimation method based on PSO, while, on the other hand, to design a SMC controller, both things used for the first time in this problem and also for speeds ranging from 2 to 14 km/h. In this way, the estimation of the model's parameters $X=(x_i)=[ a_1..... a_6]$ is formulated as an optimization problem. Hence, the optimization of a cost function will provide an estimation of the parameters. The PSO algorithm will be used to solve the so-obtained optimization problem. \label{ss.PSO} \section{Parameter Estimation and Controller design} This section contains the description of the PSO algorithm along with the derivation of the sliding mode controller. \label{s.PSO} \subsection{Parameter estimation of the model} PSO algorithm was utilized for understanding the regulations dominating the swarms of birds and their sudden changes. PSO consists of a population (or swarm) of $M$ particles, each of which represents a $ n $ dimensional potential solution of the optimization problem. In our approach, $n=6$ is the number of parameters to be estimated. Particles are assigned random initial positions and they change their positions iteratively to reach the global optimal solution. It is desired to minimize the fitness function as the PSO iterations progress. The parameter estimation problem is cast into an optimization one so that the minimum of the fitness cost function will provide an estimation of the parameters of the system. The Squared Error Loss (SEL) is the most common cost function to be optimized for speed estimation problems and it is also the easiest to work with from a mathematical point of view. The SEL is linked with variance and bias of an estimator, so that the cost function is formulated for the estimated parameter vector ${ { \hat { X } }_{ k } }$ at iteration step $k$ as : \begin{equation} { I }_{ k }= Var({ \hat { y } }_{ k })+ bias(\hat { X } _{ k }) \end{equation} where $\hat{y}_{k}$ is the estimated output. Both terms in Eq. (2) are nonnegative i.e. $(Var(δ\hat { y } _{ k })>0 , biasδ(\hat { { X }_{ k } }) \ge 0)$, so that the minimum of the cost function $I_k$ is given by $I_k =0$, $argmin{ I }_{ k }=\left( 0,0 \right) $. Therefore, when the cost function $I_k$ vanishes then $(Var δ \hat{y}_{ k } = 0)$ and $biasδ(\hat { { X }_{ k } })=0$ implying that the estimation of the parameters is performed adequately. The minimization of such cost function is done by using the PSO algorithm. Each particle evaluates its fitness (given by Eq. (5)) and every particle $i=[1,..., M]$ has a memory to store the value of its best own position $Pbest$ $_i$$_d$, which is defined as the position where the particle has minimum fitness. Besides, the best of $Pbest$ $_i$$_d$ of all particles, called $Gbest_d$, is stored too. At each iteration $k$, the PSO modifies each dimension of the position $x $$_i$$_d$ in a particle by adding a velocity $v_i$$_d$ and moves the particle towards the linear combination of $Pbest$ $_i$$_d$ and $Gbest_d$ according to: \begin{eqnarray} %\begin{aligned} v _{id} (k+1) &=& w\; {v}_{id} (k) + c_1 \; rand_1 (p_{id} - x_{id})\;+c_2 \; rand_2 (p_{gd} -x_{id}) \nonumber \\ \label{eq.6} \\ x_{id} (k+1) &=& x_{id} (k) + v_{id} (k+1) \label{eq.7} %\end{aligned} \end{eqnarray} In fact, according to Eqs. (\ref{eq.6})-(\ref{eq.7}), some new particles may be out of the search space so that a projection to the boundaries of such space is included in the algorithm, (Yan et al, 2013). Moreover, the most common approach to restrict the particle position in the search space is to set the violated components of the particle equal to the value of the violated boundary. In the problem at hand, the constraint violation appears when the algorithm provides a negative value for the parameters. Consequently, the projection algorithm takes the form: \begin{eqnarray} { x }_{ id }=\begin{cases} 0\quad ,\quad { x }_{ id }<0 \\ { x }_{ id },\quad otherwise \end{cases} \end{eqnarray} In this way, we can guarantee that the estimated parameters are nonnegative and the velocity and position of each particle are updated by the equations until a termination condition is met and the algorithm finally stops. The parameters $c_1$ and $c_2$ are the so-called cognitive and social parameters, respectively, and satisfy $0 {a_2}| \tilde{x}_2(t) - f_{uncer1}(x)| \label{eq.conK1} \end{equation} Condition (\ref{eq.conK1}) can be further elaborated in the following way. The dynamics of the observation error is obtained by Eq. (27), whose result is: \begin{equation} \dot{\tilde{x}}_2(t) = -a_4 \tilde{x}_2(t) - f_{under2}(x) \label{eq.obserror} \end{equation} The solution to this equation is given by: \begin{equation} \tilde{x}_2(t) = e^{-a_4 t} \tilde{x}_2(0) - \int_0^t e^{-a_4(t - \tau)} f_{uncer2}(x) d\tau \label{eq.solobs} \end{equation} If the uncertain function $f_{uncer2}$ is upper-bounded, i.e. $\sup |f_{uncer2}|<\infty$, fact that holds since according to Assumption 1 the uncertainty terms are bounded, then (\ref{eq.solobs}) can be upper-bounded accordingly as: \begin{eqnarray} |\tilde{x}_2(t)| &\leq& e^{-a_4 t} |\tilde{x}_2(0)| + \frac{1}{a_4} \sup |f_{uncer2}(x)| \left( 1 - e^{-a_4 t} \right) \label{eq.upperobs} \\ &\leq& |\tilde{x}_2(0)| + \frac{1}{a_4} \sup |f_{uncer2}(x)| \label{eq.upperobs2} \end{eqnarray} for all $t \geq 0$. In this way, (\ref{eq.conK1}) is satisfied if the following condition holds: \begin{equation} K a_3 > \left( a_2 |\tilde{x}_2(0)| + \frac{a_2}{a_4} \sup |f_{uncer2(x)}| + \sup |f_{uncer1}(x)| \right) \label{eq.condK2} \end{equation} since \begin{eqnarray} && a_2 |\tilde{x}_2(0)| + \frac{a_2}{a_4} \sup |f_{uncer2(x)}| + \sup |f_{uncer1}(x)| \nonumber \\ &\geq& a_2 |\tilde{x}_2(t)| + \sup |f_{uncer1}(x)| \nonumber \\ &\geq& |a_2 \tilde{x}_2(t) - f_{uncer1}(x)| \label{eq.justi} \end{eqnarray} Thus, the switching gain $K$ must be selected to fulfill (\ref{eq.condK2}), a condition that depends on an upper-bound of the observation error and upper bounds of the uncertainties. In the end, we must bear in mind that we are working with the square of the control signal so that the actual speed command is given by: \begin{equation} u_{actual} = \sqrt { max(0,u) } \label{eq.condiU} \end{equation} In the case of recovery and training programs, the controller is able to make the heart rate follow the predefined profile set up. In the next section, the simulation results showing the performance obtained by the estimation and control algorithms are presented. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Simulation results} This Section is composed of three subsections. First, the estimation results obtained by means of the PSO algorithm are discussed in section {4.1}. Secondly, the control results obtained by using the super-twisting control law are presented in section {4.2}, while the comparison between the PSO and other parameter estimation procedures and SMC comparison with PID is presented in Section {4.3}. \subsection{ PSO parameter estimation results} In this part, we will apply the PSO parameter estimation algorithm to the data described in tables below corresponding to ten subjects. These data are numerical data used for testing the algorithm and they do not correspond to real subjects. The parameters of the PSO algorithm are given by $c_1= 0.87$, $c_2=0.67$, $w=58$, $rand_1=0.1$, $rand_2=0.5$ and $\varepsilon \ge 0.15$ then tolerance is small. The actual parameters of the ten subjects are given in Tables 1 and 2 and the estimated ones obtained from the PSO proposed approach are in Tables 3 and 4. \begin{table}[htbp] \centering \caption{Actual parameters of subjects 1 - 5.} \begin{tabular}{|c|c||c|c|c||c|c|} Parameters & Subject 1 & Subject 2 & Subject 3 & Subject 4 & Subject 5 \\ a1 & 2.512 & 2.791 & 2.683 & 2.592 & 2.611 \\ a2 & 25.92 & 25.74 & 25.25 & 25.41 & 25.84 \\ a3 & 0.81 & 0.85 & 0.79 & 0.8 & 0.81 \\ a4 & 0.9021 & 0.9087 & 0.909 & 0.9011 & 0.9108 \\ a5 & 0.038 & 0.041 & 0.035 & 0.039 & 0.042 \\ a6 & 5.37 & 5.51 & 5.43 & 5.65 & 5.29 \\ HRrest & 64 & 69 & 66 & 68 & 69 \\ \end{tabular}% \label{tab:addlabel}% \end{table}% \begin{table}[htbp] \centering \caption{Actual parameters of subjects 6 - 10.} \begin{tabular}{|c|c||c|c|c||c|c|} Parameters & Subject 6 & Subject 7 & Subject 8 & Subject 9 & Subject 10 \\ a1 & 3.12 & 3.7 & 3.05 & 3.25 & 3.9 \\ a2 & 21.25 & 21.95 & 21.1 & 21.55 & 21.8 \\ a3 & 1.5 & 1.6 & 1.3 & 1.4 & 1.7 \\ a4 & 1.9 & 1.8 & 1.8 & 1.25 & 1.9 \\ a5 & 1.01 & 1.21 & 1.16 & 1.09 & 1.25 \\ a6 & 8.15 & 8.35 & 8.5 & 8.6 & 8.3 \\ HRrest & 62 & 67 & 64 & 68 & 61 \\ \end{tabular}% \label{tab:addlabel}% \end{table}% \begin{table}[htbp] \centering \caption{Estimated parameters obtained by running the PSO algorithm of subjects 1 - 5.} \begin{tabular}{|c|c||c|c|c||c|c|} Parameters & Subject 1 & Subject 2 & Subject 3 & Subject 4 & Subject 5 \\ a1 & 2.508 & 2.789 & 2.681 & 2.588 & 2.615 \\ a2 & 25.88 & 25.69 & 25.2 & 25.48 & 25.71 \\ a3 & 0.849 & 0.855 & 0.788 & 0.798 & 0.809 \\ a4 & 0.9011 & 0.9079 & 0.9088 & 0.9019 & 0.9098 \\ a5 & 0.04 & 0.043 & 0.039 & 0.032 & 0.04 \\ a6 & 5.41 & 5.49 & 5.47 & 5.59 & 5.59 \\ \end{tabular}% \label{tab:addlabel}% \end{table}% \begin{table}[htbp] \centering \caption{Estimated parameters obtained by running the PSO algorithm of subjects 6 - 10.} \begin{tabular}{|c|c||c|c|c||c|c|} Parameters & Subject 6 & Subject 7 & Subject 8 & Subject 9 & Subject 10 \\ a1 & 2.95 & 3.15 & 2.8 & 2.95 & 3.5 \\ a2 & 20.8 & 21.5 & 20.5 & 20.85 & 21.25 \\ a3 & 1.1 & 1.2 & 0.95 & 1.15 & 1.2 \\ a4 & 1.5 & 1.3 & 1.2 & 0.95 & 1.5 \\ a5 & 0.8 & 0.9 & 1.0 & 0.85 & 1.0 \\ a6 & 7.75 & 7.9 & 8.2 & 8.1 & 7.9 \\ \end{tabular}% \label{tab:addlabel}% \end{table}% As it can be noticed from Tables 1, 2, 3 and 4, the estimated parameters obtained by the proposed PSO approach are close to the actual ones. These results mean that the PSO algorithm performs very well. Figure 2 shows the actual heart rate corresponding to one of our subjects (subject No.3) along with the output of the estimated model obtained by using the PSO algorithm described in Section {3.1}. As it can be observed in this figure, the estimated model captures the dynamics of the heart rate of each person. It means that the PSO algorithm has a good response. Moreover, PSO is able to provide accurate parameter values and able to reproduce the behavior of the heart rate. Figure 3 displays the HR generated by the original model parametrized by the parameters in Table 1 along with the values obtained when the estimated parameters values are used to obtain the HR response according to (1)-(4). We can observe some mistmatch between the actual output and the estimated one in Figure 3, which is due to the uncertain dynamics that can not be adequately captured by the model parameters. This unmodeled dynamics will be counteracted by means of the sliding mode control. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j53.page1.eps} \caption{Heart rate response with using PSO at speed of (2-14 km/h)} \label{fig:2} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j54.page1.eps} \caption{Heart rate response with using PSO at speed of (2-14 km/h)- for 4 random subjects} \label{fig:3} \end{center} \end{figure} The following figures (Fig. 4 and Fig. 5) display the evaluation of the estimated parameters for subjects No.1 and No.8 in Table 1 and the estimated onse. These figures show that after a small number of iterations the estimated parameters are close to the actual ones, the fact that is displayed numerically in Tables 2, 3 and 4. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j68.page1.eps} \caption{Evaluation of estimated parameter subject No.1} \label{fig:4} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j69.page1.eps} \caption{Evaluation of estimated parameter subject No.8} \label{fig:5} \end{center} \end{figure} The following Figure 6 displays the actual value of parameters for the ten subjects. Despite the value of the actual parameters with exhibit a large variability, showing that the algorithm works in a variety of situations the parameters. Thus, we have Figure 7, which is the relative error of the parameters with high variability. As it can be observed, the proposed PSO algorithm is able to achieve a superb estimation since the relative error is given by $\frac { { a }-{ a }_{ estimation } }{ a } \times 100$ which is very low and algorithm with this high variabilities has a superb response. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{Jr(4).eps} \caption{Actual value of parameter.} \label{fig:6} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{Jr(7).eps} \caption{Relative Error of parameter.} \label{fig:7} \end{center} \end{figure} \subsection{ Control Results} In this part, we choose one person (subject No.3) to show the result. We want to highlight that similar good results are obtained for all the parameters. The controller parameters are $\alpha =0.5$ and $K=10$. In Fig. 8, the actual Heart rate (HR) and the reference signal are shown, In this figure (Fig. 8), the output and the reference signal are super-impressed implying that the control objective has been achieved. The zoom of first 150 seconds in the previopus figure (Fig. 8) is shown in Figure 9. On the other hand, Figure 10 shows the speed calculated from the SMC given by Eq. (10). In this case, the tracking error after the reaching phase is very small, despite the changes in the reference signal and it shows that how the SMC had a great response regards to the absence of chattering in the output and in the control command. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j75.page1.eps} \caption{Heart Rate provided by the SMC controller.} \label{fig:8} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j76.page1.eps} \caption{The first 150 seconds HR provided by the SMC controller.} \label{fig:9} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j74.page1.eps} \caption{Speed provided by the SMC controller.} \label{fig:10} \end{center} \end{figure} \subsection{Estimation and control comparisons} Previously in some studies, researchers used different methods such as the M-estimator to solve the parameter estimation problem, (Peter and Huber, 1964). This method is effective, from the statistical point of view, to obtain an adequate estimation of the parameters. The comparison with the M-estimator procedure shows that PSO has better behavior than the previous approach, while the accuracy of the estimation parameter is increased by using PSO. In this subsection, the theory of M-estimation is introduced for comparison purposes. The previous researchers used robust performance of estimation for two main reasons, namely: 1)there may be outliers in the data, that are sample values considered very different from the majority of the sample and 2)the data may depart from the underlying distribution assumptions, (Cheng et al, 2008), (Peter and Huber, 1964). This method is good for estimating the parameters, which is the reason why we compare our approach with this one. The class of M-estimators contains the maximum likelihood estimator (ML) as a special case. If we assume that the data come from the model distribution $F (\mu,\sigma)$ then the log-likelihood can be written as: \begin{equation} \sum _{ i=0 }^{ n }{ \{ log({ f }_{ 0 } } (\frac { { x }_{ i }-{ \mu } }{ \sigma } )-log\sigma ) \end{equation} The first order condition for the M-estimator of $\mu$ is then given by: \begin{equation} \frac { 1 }{ n } \sum _{ i=0 }^{ n }{ { \psi }_{ M } } (\frac { { x }_{ i }-{ \mu }_{ M } }{ \sigma } )=0 \end{equation} while the M-estimator of scale verifies \begin{equation} \frac { 1 }{ n } \sum _{ i=0 }^{ n }{ { \rho }_{ M } } (\frac { { x }_{ i }-{ \mu } }{ { \sigma }_{ M } } )=1 \end{equation} with $\psi _ M (u)$ being the so-called score function, and ${ \rho }_{ M }(u)={ \psi }_{ M }(u)u$. Under regularity conditions the ML estimators have a 100 percent efficiency, meaning that their asymptotic variance equals the inverse of the Fisher information, the lower bound of the Cramer-Rao inequality, (Tian et al,2014), (Peter and Huber, 1964). The parameters are coming from the particular subject (subject No.3) of Table 1, estimated by solving Eq. (17) with respect to $x_i$. The estimated parameters obtained by using the M-estimator are applied again to parametrize Eq. (1)-(4). Figure 11 demonstrates that in many points the output of PSO intersects with the reference signal and it means that PSO estimation is really close to the reference. On the other hand, the M-estimator displays an output that intercepts at some points with the reference, it has many variations in most of the points and this is not as effective as PSO. Whereas, we are calculating the fit-in error (since it is an error coming from a difference in the output of the models and actual data) in this part, that is only in open loop and it is only for estimation purposes. The fit-in error is calculated as: \begin{eqnarray}\centering %\begin{aligned} { J }_{ 1,k } &=& [{ r }_{ k }-{ p }_{ k }]\\ { J }_{ 2,k } &=& [{ r }_{ k }-{ m }_{ k }] %\end{aligned} \end{eqnarray} where : ${ r }_{ k }$ = $reference$ ${ p }_{ k }$ = $PSO$ $output$ ${ m }_{ k }$ = $M-estimator$ $parameter$ $output$ Fig. 12 shows the fit-in results of the output in the PSO and M-estimator. We want to show that PSO outperforms the M-estimator. As it can be seen at the beginning of the process in Figure 12, PSO after 500 seconds starts definitely better than the M-estimator and it slightly goes better afterward. Both of these models have a good response, but PSO performs better during the fit-in error process. On the other hand, PSO has a lower error against the M-estimator. Finally, in Fig. 13 the value of the cost function (Eq. (5)) in Particle Swarm Optimization and M-estimator are displayed. The cost function (2) in the PSO after 5 iterations reduces faster than the M-estimator. So the fact that can be interpreted as that estimation is performed faster in PSO method against the M-estimator, and this fact is reflected in the quality of estimation. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j61.page1.eps} \caption{Comparing heart rate tracking result with PSO and M-estimator.} \label{fig:11} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j62.page1.eps} \caption{Fit-in error for the PSO and M estimators.} \label{fig:12} \end{center} \end{figure} \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{j58.page1.eps} \caption{Comparing error with PSO and M-estimator.} \label{fig:13} \end{center} \end{figure} On the other hand, PID is a common approach in the control of systems. For this reason, it is used to solve the control problem in many studies. The SMC control will be compared with the PID controller implemented in (Girard et al, 2016). Since PID controllers are widely used in practice we will show the results achieved by the proposed controller in this scenario. In Fig. 14, the comparison of the actual HR obtained by using the SMC and the PID controller is shown. As it is clear, the SMC works much better than PID showing that it is able to obtain a zero tracking error despite the presence of uncertainties in the system's model. Overall, the presented method is able to obtain an appropriate and superb closed-loop behaviour. \begin{figure}[h] \begin{center} \includegraphics[height=5.9cm,width = 8.9cm]{Jr(2).eps} \caption{Heart rate Comparision with SMC and PID.} \label{fig:14} \end{center} \end{figure} \section{CONCLUSIONS} This paper has considered the design of a sliding mode controller for the HR control during treadmill exercise. Initially, a Particle Swarm Optimization algorithm has been proposed to obtain an accurate estimation of the model parameters. Secondly, a super-twisting based sliding control law has been designed for the system in order to counteract the remaining potential unmodelled dynamics or parametric uncertainty in the system. In both situations, the range of treadmill speeds goes from 2 up to 14 km/h, range not usually employed in previous studies. Simulation results show how the proposed PSO algorithm is able to obtain accurate estimations of the model parameters while the super-twisting sliding controller is capable of obtaining zero tracking error without chattering. \begin{ack} {This work was partially supported by the Spanish Ministry of Economy and Competitiveness through grant DPI2016-77271-R and by the University of the Basque Country (UPV/EHU) through grant PPG17/33.} \end{ack} %\bibliography{ifacconf} %bib file to produce the bibliography % with bibtex (preferred) \begin{thebibliography}{50} % you can also add the bibliography by hand \bibitem{c1} Hunt, Kenneth J., and Simon E. Fankhauser. "Heart rate control during treadmill exercise using input-sensitivity shaping for disturbance rejection of very-low-frequency heart rate variability." Biomedical Signal Processing and Control 30 (2016): 31-42. \bibitem{c2} Weippert, M., Behrens, M., Rieger, A., Behrens, K. Sample entropy and traditional measures of heart rate dynamics reveal different modes of cardiovascular control during low intensity exercise,(2014). Entropy, 16(11), 5698-5711. \bibitem{c3} Cheng TM, Savkin AV, Celler BG, Su SW, Wang L. Heart rate regulation during exercise with various loads: identification and nonlinear H infinity control. InIFAC World Congress 2008. 2008 IFAC. \bibitem{c4} Yan, X., Wu, Q., Liu, H., & Huang, W. (2013). An improved particle swarm optimization algorithm and its application. International Journal of Computer Science Issues (IJCSI), 10(1), 316. \bibitem{c5} Jang, D. G., Ko, B. H., Sunoo, S., Nam, S. S., Park, H. Y., & Bae, S. K. (2016, August). A preliminary study of a running speed based heart rate prediction during an incremental treadmill exercise. In Engineering in Medicine and Biology Society (EMBC), 2016 IEEE 38th Annual International Conference of the (pp. 5323-5326). IEEE. \bibitem{c6} Su SW, Huang S, Wang L, Celler BG, Savkin AV, Guo Y, Cheng TM. Optimizing heart rate regulation for safe exercise. Annals of biomedical engineering. 2010 Mar 1;38(3):758-68. \bibitem{c7} Zhang, Yudong, and Lenan Wu. "Crop classification by forward neural network with adaptive chaotic particle swarm optimization." Sensors 11.5 (2011): 4721-4743. \bibitem{c8} Tian, Xiaomin, and Shumin Fei. "Robust control of a class of uncertain fractional-order chaotic systems with input nonlinearity via an adaptive sliding mode technique." Entropy 16.2 (2014): 729-746. \bibitem{c9} Cheng TM, Savkin AV, Celler BG, Su SW, Wang L. Nonlinear modeling and control of human heart rate response during exercise with various work load intensities. Biomedical Engineering, IEEE Transactions on. 2008 Nov;55(11):2499-508. \bibitem{c10} Peter J. Huber, Robust Estimation of a Location Parameter , The Annals of Mathematical Statistics, Vol. 35, No. 1. (Mar., 1964), pp.73-101. \bibitem{c11} Maronna, R.A. (1976). Robust M-estimators of Multivariate Location and Scatter, The Annals o Statistics 4, 51-67. \bibitem{c12} Rini DP, Shamsuddin SM, Yuhaniz SS. Particle swarm optimization: Lari A, Khosravi A, Rajabi F. Controller design based on $\mu$ analysis and PSO algorithm. ISA transactions. 2014 Mar 31;53(2):517-23. \bibitem{c13} Liu, W., Liu, L., Chung, I. Y., and Cartes, D. A. (2011). Real-time particle swarm optimization based parameter identification applied to permanent magnet synchronous machine. Applied Soft Computing, 11(2), 2556-2564. \bibitem{c14} Bansal, J. C., Singh, P. K., Saraswat, M., Verma, A., Jadon, S. S., and Abraham. "Inertia weight strategies in particle swarm optimization." Nature and Biologically Inspired Computing (NaBIC),IEEE, 2011. 2011 Third World Congress on. \bibitem{c14} Eswaran, T., and V. Suresh Kumar. "Particle swarm optimization (PSO)-based tuning technique for PI controller for management of a distributed static synchronous compensator (DSTATCOM) for improved dynamic response and power quality." Journal of Applied Research and Technology 15.2 (2017): 173-189. \bibitem{c15} Hamed, E., Anushya, A., Alzoubi, R., and Vincy, B. A. (2017, August). An analysis of particle swarm optimization for feature selection on medical data. In 2017 International Conference on Energy, Communication, Data Analytics and Soft Computing (ICECDS) (pp. 227-231). IEEE. \bibitem{c16} Ebrahimi, N., Ozgoli, S., and Ramezani, A. (2018). Model-free sliding mode control, theory and application. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 0959651818780597. \bibitem{c17} Meyer, Daniel, and Veit Senner. "Evaluating a heart rate regulation system for human–electric hybrid vehicles." Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology 232, no. 2 (2018): 102-111. \bibitem{c18} Scalzi, S., Tomei, P., and Verrelli, C. M. (2012). Nonlinear control techniques for the heart rate regulation in treadmill exercises. IEEE Transactions on Biomedical Engineering, 59(3), 599-603. \bibitem{c19} Esmaeili, A., and Ibeas, A. (2016, October). Particie Swarm Optimization modelling of the heart rate response in treadmill exercise. In System Theory, Control and Computing (ICSTCC), 2016 20th International Conference on (pp. 613-618). IEEE. \bibitem{c20} Ibeas, A., Esmaeili, A., Herrera, J., and Zouari, F. (2016, October). Discrete-time observer-based state feedback control of heart rate during treadmill exercise. In System Theory, Control and Computing (ICSTCC), 2016 20th International Conference on (pp. 537-542). IEEE. \bibitem{c21} Shtessel, Yuri B. Super-twisting adaptive sliding mode control: A Lyapunov design. In: Decision and Control (CDC), 2010 49th IEEE Conference on. IEEE, 2010. p. 5109-5113. \bibitem{c22} Shtessel, Yuri B., Jaime A. Moreno, and Leonid M. Fridman. "Twisting sliding mode control with adaptation: Lyapunov design, methodology and application." Automatica 75 (2017): 229-235. \bibitem{c23} GAO, Xuehui. Variable gain super-twisting sliding mode control for Hammerstein system with Bouc-Wen hysteresis nonlinearity. In: Control Conference (CCC), 2016 35th Chinese. IEEE, 2016. p. 3369-3372. \bibitem{c24} Kranjec, Beguš, Geršak , Drnovšek. "Non-contact heart rate and heart rate variability measurements: A review." Biomedical Signal Processing and Control 13 (2014): 102-112. \bibitem{c25} Lu, C. H., Wang, W. C., Tai, C. C., and Chen, T. C. (2016). Design of a heart rate controller for treadmill exercise using a recurrent fuzzy neural network. Computer methods and programs in biomedicine, 128, 27-39. \bibitem{c26} Abu-Rmileh, A., Garcia-Gabin, W., and Zambrano, D. (2010). Internal model sliding mode control approach for glucose regulation in type 1 diabetes. Biomedical Signal Processing and Control, 5(2), 94-102. \bibitem{c27} Anselmino, M., Scarsoglio, S., Saglietto, A., Gaita, F., and Ridolfi, L. "A computational study on the relation between resting heart rate and atrial fibrillation hemodynamics under exercise." PloS one 12.1 (2017): e0169967. \bibitem{c28} Yuan. B, Gallagher. M, (2003), Playing in Continuous Spaces: Some Analysis and Extension of Population-Based Incremental Learning, IEEE Evolutionary Computation, vol: 1, pp: 443-450. \bibitem{c29} Swikir, Abdalla, and Vadim Utkin. "Chattering analysis of conventional and super twisting sliding mode control algorithm." Variable Structure Systems (VSS), 2016 14th International Workshop on. IEEE, 2016. \bibitem{c30} Wang, X., Wang, S., and Ma, J. J. (2007). An improved co-evolutionary particle swarm optimization for wireless sensor networks with dynamic deployment. Sensors, 7(3), 354-370. \bibitem{c31} Argha, A., Su, S. W., and Celler, B. G. (2017). Heart rate regulation during cycle-ergometer exercise via event-driven biofeedback. Medical & biological engineering & computing, 55(3), 483-492. \bibitem{c32} Girard, C., Ibeas, A., Vilanova, R., and Esmaeili, A. (2016, May). Robust discrete-time linear control of heart rate during treadmill exercise. In Electrical Engineering (ICEE), 2016 24th Iranian Conference on (pp. 1113-1118). IEEE. \bibitem{c33} Vaidyanathan, Sundarapandian. "Super-Twisting Sliding Mode Control of the Enzymes-Substrates Biological Chaotic System." Applications of Sliding Mode Control in Science and Engineering. Springer International Publishing, 2017. 435-450. \bibitem{c34} Dan, A. M., and Dragomir, T. L. (2012, May). A model of the control function in the case of constant effort of the cardiovascular system. In Applied Computational Intelligence and Informatics (SACI), 2012 7th IEEE International Symposium on (pp. 169-173). IEEE. \bibitem{c35} Belkaid, A., Gaubert, J. P., and Gherbi, A. (2016). An improved sliding mode control for maximum power point tracking in photovoltaic systems. Journal of Control Engineering and Applied Informatics, 18(1), 86-94. \bibitem{c36} Xue, B., Zhang, M., and Browne, W. N. (2013). Particle swarm optimization for feature selection in classification: A multi-objective approach. IEEE transactions on cybernetics, 43(6), 1656-1671. \bibitem{c37} Storn R, Price K (1997) Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces. J Glob Optim 11:341–359 \bibitem{c38} Suganthan PN (1999) Particle swarm optimizer with neighborhood operator. In: Proceedings of the IEEE congress on evolutionary computation, Washington, DC, pp 1958–1961 \bibitem{c39} Souravlias D, Parsopoulos KE (2016) Particle swarm optimization with neighbourhood-based budget allocation. Int J Mach Learn Cybern 7(3):451–477. \bibitem{c40} Yang, H., Xu, Y., Peng, G., Yu, G., Chen, M., Duan, W., and Wang, X. (2017). Particle swarm optimization and its application to seismic inversion of igneous rocks. International Journal of Mining Science and Technology, 27(2), 349-357. \bibitem{c41} Paul, Arun K., J. E. Mishra, and M. G. Radke. "Reduced order sliding mode control for pneumatic actuator." IEEE Transactions on Control Systems Technology 2, no. 3 (1994): 271-276. \bibitem{c42} Guo, H., Liu, Y., Liu, G., and Li, H. (2008). Cascade control of a hydraulically driven 6-DOF parallel robot manipulator based on a sliding mode. Control Engineering Practice, 16(9), 1055-1068. \end{thebibliography} \end{document}