%===============================================================================
% $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.

\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.}

\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}

\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 at 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 modelling 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), (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 behaviour 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, 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). SMC also these days is very useful in many studies such a pneumatic cylinders as actuators for robot manipulators, (Paul et al, 1994) and  the hydraulic dynamics of the manipulator, (Guo et al, 2008). 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. This is the first time that a 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).  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 is defined for anybody, 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 also one alternative estimation procedure to compare with our results.
\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}
\begin{aligned}
\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} 
\end{aligned}
\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  hyperaemia  in  active  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. $ 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  $  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  $  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, and 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  an $ n $ dimensional  potential  solution. 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 using the PSO algorithm.

 Each particle evaluates its fitness (2) 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. 
In this case, the underlying stopping condition takes the form. 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 \\
 x_{id} (k+1) &=& x_{id} (k) + v_{id} (k+1)
  \nonumber 
 %\end{aligned}
 \end{eqnarray}
In fact, according to Eqs.(6), 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). Although, the most common approach to restrict the particle in the search space is to set the violated components of the particle equal to the value of the violated boundary :
 \begin{eqnarray}
 { x }_{ id }=\begin{cases} 0\quad ,\quad { x }_{ id }<0 \nonumber \\ { 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 and repeated by the equations, so that the maximum number of repetitions can be reached or the speed upgraded approaches zero and the performance of each particle is estimated by the merit criterion. The  parameters  $c_1$ and  $c_2$ are the so-called cognitive  and social  parameters,  respectively, and satisfy  $0<c_1,c_2<1$, [Rini et Al, 2014]. Finally,  $w$  is  the  inertia weight selected as $0.4\le w\le 0.9$, range where the algorithm provides the best results, (Bansal et al, 2011).  
 
  Typically, the number of subsequent iterations without improvement of the best solution and/or the dispersion of the particles current (or best) positions in the search space have been used as indicators of search stagnation.
Frequently, the aforementioned termination criteria are combined in forms such as:
 \begin{equation}
IF(\left| { I }_{ k+1 }-{ I }_{ k } \right| \le \varepsilon )\quad OR(\quad k\ge { k }_{ max })).ThenStop
 \end{equation}
where ${ I }_{ k }$ is the targets in the search space and function values and $k$ stands for the iteration number, respectively, and $ \varepsilon$ is the corresponding user-defined tolerance. However, the search stagnation criterion can prematurely stop the algorithm even if the computational budget is not exceeded. Successful application of this criterion is based on the existence of a proper stagnation measure.
 
 Figure 1 displays the pseudocode of the PSO algorithm developed in our study. 
\begin{figure}[h]
\begin{center}
\includegraphics[height=12cm,width = 8.5cm]{PSOPsucode.page1.eps}
\caption{PSO pseudocode employed to solve the optimization problem.}
\label{fig:1}
\end{center}
\end{figure}
  The PSO search is carried out by the speed of the particle. During the development of several generations, only the most optimistic particles can transmit information to the other ones. One of the advantages of the PSO method is that it can be applied to optimization problems of large dimensions, often producing quality solutions more rapidly than alternative methods. Also, the algorithm is terminated after a given number of iterations, or once the fitness values of the particles (or the particles themselves) are close enough in some sense. 
 
 In the results section (Section 4), we describe the M-estimator as an alternative model to compare with our estimation results, showing that the proposed PSO method outperforms the M-estimator. Once the parameter estimation procedure has been described, the next step is to design a controller to make the heart rate track a pre-specified  profile.  The design of such a controller is carried out in the following section. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Super-twisting sliding mode control}
The design of control strategies for nonlinear systems has attracted considerable research interest in the recent past, (Shtessel and Yuri, 2017). Sliding
mode control (SMC), as an effective robust control scheme, has been successfully applied to a wide variety of systems, (GAO and Xuehui, 2016),(Abu-Rmileh et al, 2010). This section contains the design of the sliding mode control for the system (1). Thus, define the tracking error as:
 \begin{equation}
e=R-y 
\label{eq.trackerror}
\end{equation}
 where $R$ denotes the reference signal (that is, the HR profile to be tracked) and $y$ is the output of our system. The role of the controller is to ensure that system's output accurately tracks the reference signal $R$. When the system is perturbed or uncertain, the finite time stabilization is not ensured. Hence, a reaching law based discontinuous control is developed which rejects the uncertainties of the system and ensures that the control objectives are fulfilled. The uncertainties in our system can be modelled as:
 \begin{eqnarray}
\begin{aligned}
\dot x_1(t)  =& - a_1 x_1(t) +  a_2 x_2(t) + a_3 u(t)^2 + f_{uncer1}(x)\\
\dot x_2(t)  =& - a_4 x_2 (t) +\phi(x_1(t)) + f_{uncer2}(x) \\
\phi (x_1(t)) =& \frac{{ a }_{ 5 }{ x }_{ 1 }(t) }{1 + \exp{(-(  x_1(t) -  a_6))}} \\
\  y(t)  &= x_1(t) + HR_{rest}
\end{aligned}
\label{eq.modeuncertain}
\end{eqnarray}
where $f_{uncer1}(x)$ and $f_{uncer2}(x)$ account for the unmodelled dynamics and parametric uncertainty in each of the model equations. On the other hand we need to consider the following assumptions.

$Assumption 1$: $f_{uncer1}(x)$ and $f_{uncer2}(x)$ are upper-bounded \nonumber \\
$Assumption 2$: One upper-bound for each one of these terms is known \nonumber 
  
These are common assumptions in SMC, (Shtessel and Yuri, 2010). The sliding mode controller is composed of two parts:
\begin{equation}
u={ u }_{ equiv }+{ u }_{ sliding }
\label{eq.modetotal}
\end{equation}
where $u_{equiv}$ is the so-called equivalent control used to remove certain terms in (10) while the sliding term $u_{sliding}$ is the term used to counteract the uncertainties of the system and will be of the super-twisting type, (Shtessel and Yuri, 2017). This approach will also help us avoid the chattering effect, that would be very harmful in the control system. Initially, the equivalent control will be derived while the final control law will be obtained by incorporating the super-twisting sliding term according to (\ref{eq.modetotal}).

The following sliding manifold with the integral term is proposed: 
\begin{equation}
 S(t) =& e{(t)} +{ \lambda  }\int _{ 0 }^{ t }{e(\tau) \quad d\tau  } } 
\label{eq.slidingsurface}
\end{equation}
 where  ${ \lambda  }$ are strictly positive constant. The equivalent control is obtained by derivating (\ref{eq.slidingsurface}) with respect to time and then equating the so-obtained derivative to zero. In this way we have:  
 \begin{eqnarray}\centering
\begin{aligned}
 e(t)= R(t)-y(t) = R(t)-{ x }_{ 1 }(t) \nonumber \\
  \dot { e } (t)&=& \dot { R } (t)-\dot { { x }_{ 1 } } (t) \nonumber\\
   = \dot { R } (t)+{ a }_{ 1 }{ x }_{ 1 }(t)-{ a }_{ 2 }{ x }_{ 2 }(t)-{ a }_{ 3 }{ u }^{ 2 }(t)+f_{ uncer1 }(t)
\end{aligned} 
\end{eqnarray}
In this way, if we substitute the above expressions into (\ref{eq.trackerror}) and simplify we obtain. Now, the derivative of the sliding manifold reads: 
\begin{eqnarray}\centering\raggedleft
\begin{aligned}
\dot { S } (t)=\dot { e } (t)+\lambda e(t)
\end{aligned}
\label{eq.simplifideriS}
\end{eqnarray}
\begin{eqnarray}
\dot { S } (t)=\dot { R } -\dot { y } +\lambda e(t)=\dot { R } -\dot { { x }_{ 1 } } +\lambda e(t) 
\end{eqnarray}
\begin{eqnarray}
\dot { S } (t)=\dot { R } -(-{ a }_{ 1 }{ x }_{ 1 }+{ a }_{ 2 }{ x }_{ 2 }+{ a }_{ 3 }{ u }^{ 2 }+{ f }_{ uncer1 })+\lambda e(t)
\end{eqnarray}

 If $\dot { S } (t)=0$ we have:
\begin{eqnarray}
\dot { R } +{ a }_{ 1 }{ x }_{ 1 }-{ a }_{ 2 }{ x }_{ 2 }-{ a }_{ 3 }{ u }^{ 2 }-{ f }_{ uncer1 }+\lambda e(t)=0
\end{eqnarray}
Now, if we isolate ${ u }^{ 2 }$ we obtain: 
 \begin{eqnarray}\centering\raggedleft
\begin{aligned}
{ a }_{ 3 }{ u }^{ 2 }=\dot { R } (t)+{ a }_{ 1 }{ x }_{ 1 }-{ a }_{ 2 }{ x }_{ 2 }-{ f }_{ uncer1 }+\lambda e(t)
\end{aligned}
\end{eqnarray} 
\begin{eqnarray}
\begin{aligned}
{ u }^{ 2 }(t)=\frac { 1 }{ { a }_{ 3 } } (\dot { R } (t)+{ a }_{ 1 }{ x }_{ 1 }-{ a }_{ 2 }{ x }_{ 2 })+\lambda e(t)
\end{aligned}
\label{eq.mode}
\end{eqnarray}
 The uncertain terms ${ f }_{ uncer1 }$ do not appear in (\ref{eq.mode}) since they are unknown. Therefore, they do not appear in the equivalent control part.
 
 The super-twisting sliding term is given by, \cite{c20}:
\begin{equation}
u_{sliding} &=& K {|S|}^{ \alpha  } sign(S)
\label{eq.slidingterm}
\end{equation}
so that the final control law is obtained from (\ref{eq.modetotal}).
  The state observer is given by:
\begin{equation}
 \dot { \hat { { x }_{ 2 } }  }  \left( t \right) =-{ a }_{ 4 }\hat { { x }_{ 2 } } \left( t \right) +\phi \left( { x }_{ 1 }\left( t \right)  \right) 
\end{equation}
 With arbitrary initial condition \hat{x}_2 (0), since ${ x }_{ 2 }$ is infusible to obtain and ${ f }_{ uncer1 }$ is unknown. The control law reads:
\begin{equation}
u(t) = \frac{1}{a_3} \left( \dot R(t)+a_1 x_1(t)-a_2 \hat{x}_2(t)+ \lambda e(t) \right) + K |S|^\alpha sign(S)
\label{eq.control}
\end{equation}

The fact is that initial condition maybe arbitrary that must useful choice  could be $\tilde{x}_2(0)=0={x}_2(0)$, Because preferably effects are zero at the starting from the rest. Although, when we starting from the zero it implies that the error is zero or close to the zero.

$Assumption 3$: Observation error at initial time is bounded and an upper bound is known \nonumber 

The switching gain $K$ has to be selected so as to guarantee the stability and reference tracking of the closed-loop system. In order to obtain a guideline for its tuning we consider the following Lyapunov function candidate:
\begin{equation}
V(t) = \frac{1}{2} S^2
\label{eq.V}
\end{equation}
Its time-derivative is given by:
\begin{eqnarray}
\dot V(t) &=& S \dot{S} = S \left( \dot R - \dot x_1(t) + \lambda e(t) \right) \nonumber \\
&=& S \left( a_2 \hat{x}_2(t) - a_2 x_2(t) - K a_3 |S|^\alpha sign(S) - f_{uncer1}(x) \right) \nonumber \\
&=& S \left( a_2 \tilde{x}_2(t) - K a_3 |S|^\alpha sign(S) - f_{uncer1}(x) \right) \nonumber \\
&=& - K a_3 |S|^{\alpha+1} + S \left( a_2 \tilde{x}_2(t) - f_{uncer1}(x) \right)
\label{eq.der}
\end{eqnarray}
where $\tilde{x}_2(t) = \hat{x}_2(t) - x_2(t)$, represents the observation error. In order to ensure the appropriate operation of the controller, the time derivative (\ref{eq.der}) should be negative-definite, fact that is achieved if: 
\begin{equation}
K a_3 > |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.(24),  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 fulfil (\ref{eq.condK2}), condition that depends of 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, controller is able to make the heart rate follow the predefined profile set up.
  In the next section, we have all of the results that achieved from the controller. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and Simulation}

This Section is composed of three subsections. First, the parameter 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 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 the table below corresponding to six 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$ then tolerance is small.  The actual parameters of the six subjects are given in Table 1 and 2 and the estimated ones obtained from the PSO proposed approach are in Table 3 and 4. 
\begin{table}[htbp]
  \centering
  \caption{Actual parameters of 1-5 subjects.}
    \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 6-10 subjects.}
    \begin{tabular}{|c|c||c|c|c||c|c|}
   
    Parameters  & Subject#6 & Subject#7 & Subject#8 & Subject#9 & Subject#10 \\
    
    a1     & 2.691 & 2.725 & 2.540 & 2.754 & 2.711 \\
    
    a2     & 25.51 & 25.19 & 25.78 & 25.18 & 25.91 \\
    
    a3     & 0.88 & 0.83 & 0.87 & 0.82 & 0.84 \\
    
    a4     & 0.908 & 0.9071 & 0.903 & 0.9041 & 0.9101 \\
    
    a5     & 0.041 & 0.036 & 0.039 & 0.040 & 0.037 \\
    
    a6     & 5.35 & 5.52 & 5.49 & 5.65 & 5.29 \\

    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 for each person.}
    \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 for each person.}
    \begin{tabular}{|c|c||c|c|c||c|c|}
   
    Parameters  & Subject#6 & Subject#7 & Subject#8 & Subject#9 & Subject#10 \\
    
    a1    & 2.699 & 2.723 & 2.539 & 2.749 & 2.71 \\
    
    a2     & 25.47 & 25.23 & 25.72 & 25.19 & 25.89 \\
    
    a3     & 0.871 & 0.826 & 0.871 & 0.829 & 0.836 \\
    
    a4     & 0.9089 & 0.9084 & 0.9021 & 0.9032 & 0.9087 \\
    
    a5     & 0.04 & 0.0398 & 0.0399 & 0.0403 & 0.0389 \\
    
    a6     & 5.29 & 5.66 & 5.56 & 5.66 & 5.3 \\
    
    \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 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 these figures, the estimated models capture the dynamics of the heart rate of each person. It means that PSO algorithm has a good response. Moreover, PSO is able to provide accurate parameter values and able to reproduce the behaviour of the heart rate. Figure 4 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). We are showing randomly other four subjects to better understanding of the PSO's work. On the other hand, we have some mismatching in this figure 3, that is a part of the actual dynamics, implying that there is a part of the actual dynamics that cannot be captured by the model equation. 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) are display the evaluation of the estimated parameters for subjects  No.1 and No.8 in table 1 and the estimated once. These figures show that after small number of iterations the estimated parameters are close to the actual ones, fact that is the contain numerically in table 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} 
 
\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 people. Parameters of the controller are $\alpha =0.5$ and $K=10$.  In Fig.6, The actual Heart rate (HR) and the reference signal are shown, which is clear to understand the reference and the out put are close to each other and both plots are super-impressed implying that the control objective has been achieved. The tracking error of the SMC is obtained and shown in Figure 7. The error which is provided by SMC after 60 second is going close to the reference signal  and is close to the zero. On the other hand, Figure 8 shows the speed calculated from the SMC given by Eg.(14). 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.
\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:6}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[height=5.9cm,width = 8.9cm]{j76.page1.eps}
\caption{Tracking Error of HR by the SMC controller}
\label{fig:7}
\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:8}
\end{center}
\end{figure}

\subsection{Estimation comparison}
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. In comparison with the M-estimator procedure shows that PSO has better behaviour 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 are: 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 that are come from a 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). Figure 8 demonstrates that in many points the output of PSO intersects with 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 }] \nonumber \\
  { 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.9 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 in the beginning of the process in the Figure 9, PSO  after 500 seconds starts definitely better than the M-estimator and it slightly goes better afterwards. 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.10 the value of the cost function (Eq.(2)) in Particle Swarm Optimization and M-estimator are displayed. The cost function (2) in the PSO after 5 iteration 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 shows the reflects 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:9}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[height=5.9cm,width = 8.9cm]{j62.page1.eps}
\caption{Tracking error in PSO and M-estimator}
\label{fig:10}
\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:11}
\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. Entropy, 16(11), 5698-5711,(2014).
\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, Xuesong,Xuesong Yan, Qinghua Wu, Hanmin Liu and Wenzhi Huang,. "An Improved Particle Swarm Optimization Algorithm and Its Application." International Journal of Computer Science (2013): 316-324.
\bibitem{c5} Jang, Dae-Geun. "A preliminary study of a running speed based heart rate prediction during an incremental treadmill exercise." Engineering in Medicine and Biology Society (EMBC), 2016 IEEE 38th Annual International Conference of the. IEEE, 2016. 
\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 of
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, Li, Wenxin Liu, and David A. Cartes. "Particle swarm optimization-based parameter identification applied to permanent magnet synchronous motors." Engineering Applications of Artificial Intelligence 21.7 (2008): 1092-1100. 
\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), 2011 Third World Congress on. IEEE, 2011.
\bibitem{c15} Scalzi, Stefano, Patrizio Tomei, and Cristiano Maria Verrelli. "Nonlinear control techniques for the heart rate regulation in treadmill exercises." Biomedical Engineering, IEEE Transactions on 59.3 (2012): 599-603.
\bibitem{c16} Ali Esmaeili and Asier Ibeas. "Particle Swarm Optimization modelling of the heart rate response in treadmill exercise." System Theory, Control and Computing (ICSTCC), 2016 20th International Conference on. IEEE, 2016.
\bibitem{c17} Asier Ibeas, Ali Esmaeili, J.Herrera, F. Zouari.  "Discrete-time observer-based state feedback control of heart rate during treadmill exercise." System Theory, Control and Computing (ICSTCC), 2016 20th International Conference on. IEEE, 2016.
\bibitem{c18}18) 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{c19} 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{c20} 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{c21} 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{c22} Lu,Chun-Hao, Wang,Wei-Cheng, Tai,Cheng-Chi, Chen,Tien-Chi. "Design of a heart rate controller for treadmill exercise using a recurrent fuzzy neural network." Computer methods and programs in biomedicine 128 (2016): 27-39.
\bibitem{c23}  Abu-Rmileh, Amjad, Winston Garcia-Gabin, and Darine Zambrano. "Internal model sliding mode control approach for glucose regulation in type 1 diabetes." Biomedical Signal Processing and Control 5.2 (2010): 94-102.
\bibitem{c24} 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{c25} 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{c26} 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{c27} Wang, Xue, Sheng Wang, and Jun-Jie Ma. "An improved co-evolutionary particle swarm optimization for wireless sensor networks with dynamic deployment." Sensors 7.3 (2007): 354-370.
\bibitem{c28} Argha, Ahmadreza; SU, Steven W.; Celler, Branko G. Heart rate regulation during cycle-ergometer exercise via event-driven biofeedback. Medical and biological engineering and computing, 2017, 55.3: 483-492.
\bibitem{c29} Clement Girard, Asier Ibeas, Ramon Vilanova, Ali Esmaeili . Robust Discrete-time Linear Control of Heart Rate During Treadmill Exercise.24th of International Conference on Electrical Engineering,Shiraz,IRAN,ICEE 2016.
\bibitem{c30} 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{c31} Dan, Ana-Maria, and Toma-Leonida Dragomir. "A model of the control function in the case of constant effort of the cardiovascular system." Applied Computational Intelligence and Informatics (SACI), 2012 7th IEEE International Symposium on. IEEE, 2012.
\bibitem{c32} Belkaid, Abdelhakim, Jean Paul Gaubert, and Ahmed Gherbi. "An improved sliding mode control for maximum power point tracking in photovoltaic systems." Journal of Control Engineering and Applied Informatics 18.1 (2016).
\bibitem{c33} Xue, Bing, Mengjie Zhang, and Will N. Browne. "Particle swarm optimization for feature selection in classification: A multi-objective approach." IEEE transactions on cybernetics 43.6 (2013): 1656-1671.
\bibitem{c34} 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{c35} Suganthan PN (1999) Particle swarm optimizer with neighborhood operator. In: Proceedings of the IEEE congress on evolutionary computation, Washington, DC, pp 1958–1961
\bibitem{c36} Souravlias D, Parsopoulos KE (2016) Particle swarm optimization with neighbourhood-based budget allocation. Int J Mach Learn Cybern 7(3):451–477.
\bibitem{c37} Yang, Haijun, Yongzhong Xu, Gengxin Peng, Guiping Yu, Meng Chen, Wensheng Duan, Yongfeng Zhu, Yongfu Cui, and Xingjun Wang. "Particle swarm optimization and its application to seismic inversion of igneous rocks." International Journal of Mining Science and Technology 27, no. 2 (2017).
\bibitem{c38} 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{c39} Guo, Hongbo, YongGuang Liu, GuiRong Liu, and HongRen Li. "Cascade control of a hydraulically driven 6-DOF parallel robot manipulator based on a sliding mode." Control Engineering Practice 16, no. 9 (2008): 1055-1068.
\end{thebibliography


\end{document}