UNIVERSIDADE ESTADUAL PAULISTA "JÚLIO DE MESQUITA FILHO" School of Engineering of Ilha Solteira - SP Post-Graduate Program in Electrical Engineering HÉBER HWANG ARCOLEZI A Novel Robust and Intelligent Control Based Approach for Human Lower Limb Rehabilitation via Neuromuscular Electrical Stimulation Nova Abordagem Robusta e Inteligente de Controle para Reabilitação de Membros Inferiores Humanos via Estimulação Elétrica Neuromuscular Ilha Solteira 2019 HÉBER HWANG ARCOLEZI A Novel Robust and Intelligent Control Based Approach for Human Lower Limb Rehabilitation via Neuromuscular Electrical Stimulation Nova Abordagem Robusta e Inteligente de Controle para Reabilitação de Membros Inferiores Humanos via Estimulação Elétrica Neuromuscular Dissertation presented to the São Paulo State University (UNESP) - School of En- gineering - Campus of Ilha Solteira, in ful- fillment of one of the requirements for ob- taining the Master’s degree in Electrical Engineering. Specialty: Automation. Prof. Dr. Aparecido Augusto de Carvalho Advisor Ilha Solteira 2019 ARCOLEZI A Novel Robust and Intelligent Control Based Approach for Human Lower Limb Rehabilitation via Neuromuscular Electrical StimulationIlha Solteira2019 99 Sim Dissertação (mestrado)Engenharia ElétricaAutomaçãoNão . . FICHA CATALOGRÁFICA Desenvolvido pelo Serviço Técnico de Biblioteca e Documentação Arcolezi, Héber Hwang. A novel robust and intelligent control based approach for human lower limb rehabilitation via neuromuscular electrical stimulation / Héber Hwang Arcolezi. -- Ilha Solteira: [s.n.], 2019 99 f. : il. Dissertação (mestrado) - Universidade Estadual Paulista. Faculdade de Engenharia de Ilha Solteira. Área de conhecimento: Automação, 2019 Orientador: Aparecido Augusto de Carvalho Inclui bibliografia 1. Estimulação elétrica neuromuscular. 2. Controle da articulação do joelho. 3. Controlador RISE. 4. Algoritmo genético melhorado. 5. Redes neurais recorrentes. A675n ACKNOWLEDGMENTS Primarily, I am thankful to the unique God, that by love, compassion and a lot of mercy, provided me with life, health, a family and people to help me during this journey. For my family, my beloved parents and brothers, my greatest thank for each of you, that have been supporting and taking care of me during my entire life. From each of you, a different type of love has been demonstrated as the years pass by, and with pleasure, I consider and retribute it with all the love I can provide for all of you. I thank my advisor, Aparecido Augusto de Carvalho, for accepting me at the Laboratório de Instrumentação e Engenharia Biomédica. I thank him for trusting and guiding me through this work; even more for supporting me on decision-making and opportunities that I searched for. I also thank professor Marcelo Teixeira, for his trust and a lot of advice on professional and research subjects, whose the last he does it with passion. I am grateful for meeting the most special person to me, thank you, Selene, for being this person. I admire Selene for her great willing for helping and sharing with others. I learned a lot with her, both technically and life experience, which essentially helped me during this research. I am thankful for all colleagues at the research laboratory, especially Rafael and Willian who directly contributed to this work in both theoretical and practical questions. To Willian, I have a lot to thank about and he knows it, it is inspiring his willing to research and to contribute, as well as his patience for helping. I am also grateful for each person that accepted to volunteer to this research, even more for those with spinal cord injury, that with locomotion limitation disposed themselves to come at the research laboratory to contribute with this study. I thank professors Rubén Lázaro, Rodrigo Cardim, Jean Marcos, and Sérgio Kurokawa, in no strict order, for supporting and helping me in the decisions taken during this master’s de- gree. I am also grateful to Christophe Guyeux, Raphaël Couturier, and Jean-François Couchot, professors at the University of Franche-Comté, for their welcome and support in the AND team. Lastly, but not less important, I thank to UNESP, which provided the opportunity for pursu- ing this study with very qualified personnel and equipment. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001". RESUMO Nos últimos anos, vários estudos foram realizados mostrando que a estimulação elétrica neuro- muscular (EENM) pode produzir bons resultados terapêuticos em pacientes com lesão medular (LM). Esta pesquisa introduz uma nova metodologia robusta e inteligente baseada em con- trole para a reabilitação de membros inferiores humanos via EENM usando uma técnica de controle de tempo contínuo chamada robust integral of the sign of the error (RISE). Embora na literatura o controlador RISE tem demonstrado bons resultados sem qualquer método de ajuste fino, uma abordagem de tentativa e erro poderia levar rapidamente à fadiga muscular em pacientes com LM. Portanto, foi mostrado nesse estudo que o desempenho do controle para ras- trear com robustez um sinal de referência pode ser melhorado através da abordagem proposta, fornecendo um ajuste inteligente para cada voluntário. Resultados de simulação com um mod- elo matemático e oito sujeitos identificados da literatura são fornecidos, e experimentos reais são feitos com sete indivíduos saudáveis e dois paraplégicos. Além disso, esta pesquisa introduz a aplicação de redes neurais profundas e dinâmicas, especificamente o perceptron multicamadas, uma rede neural recorrente simples e a arquitetura Long Short-Term Memory, para identificar a relação não-linear e variante no tempo entre a EENM fornecida e a posição angular alcançada. Os resultados de identificação indicam boa adaptação aos dados e erro quadrático médio muito baixo usando poucos dados para treinamento, provando ser métodos muito prospectivos para propor modelos orientados ao controle. Palavras-chave: Estimulação elétrica neuromuscular. Controle da articulação do joelho. Con- trolador RISE. Algoritmo genético melhorado. Redes neurais recorrentes. ABSTRACT In the last few years, several studies have been carried out showing that neuromuscular elec- trical stimulation (NMES) can produce good therapeutic results in patients with spinal cord injury (SCI). This research introduces a new robust and intelligent control-based methodol- ogy for human lower limb rehabilitation via NMES using a continuous-time control technique named robust integral of the sign of the error (RISE). Although in the literature the RISE con- troller has shown good results without any fine-tuning method, a trial and error approach would quickly lead to muscle fatigue in SCI patients. Therefore, it was shown in this study that the control performance for robustly tracking a reference signal can be improved through the pro- posed approach by providing an intelligent tuning for each voluntary. Simulation results with a mathematical model and eight identified subjects from the literature are provided, and real ex- periments are performed with seven healthy and two paraplegic subjects. Besides, this research introduces the application of deep and dynamic neural networks namely the multilayer per- ceptron, a simple recurrent neural network, and the Long Short-Term memory architecture, to identify the nonlinear and time-varying relationship between the supplied NMES and achieved angular position. Identification results indicate good fitting to data and very low mean square error using few data for training, proving to be very prospective methods for proposing control- oriented models. Keywords: Neuromuscular electrical stimulation. Knee joint control. RISE controller. Im- proved genetic algorithm. Recurrent neural networks. LIST OF FIGURES Figure 1 The proposed robust and intelligent control-based methodology. . . . . 21 Figure 2 Genetic operations to advance one generation to the next. . . . . . . . 23 Figure 3 Real and positive coded genes of a chromosome. . . . . . . . . . . . . 24 Figure 4 Tournament selection. . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Figure 5 Fast genetic algorithm operators to perform one generation and select the best individual. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Figure 6 Parallel identification scheme of nonlinear dynamic system. . . . . . . 29 Figure 7 Series-parallel identification model of nonlinear dynamic system. . . . 29 Figure 8 Schematic representation of the lower limb with surface electrodes. . . 31 Figure 9 Normalized nominal torque with smooth, moderate and critical non- idealities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 10 Analysis of tracking responses using empiric and IGA tuning on ideal conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Figure 11 Analysis of system responses using IGA tuning considering smooth non-idealities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 12 Analysis of system responses using IGA tuning considering moderate non-idealities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Figure 13 Analysis of system responses using IGA tuning considering critical non-idealities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Figure 14 The complete test platform for electrical stimulation experiments. . . . 41 Figure 15 Open-loop tests to determine ρmin and ρmax. . . . . . . . . . . . . . . 42 Figure 16 Randomly selected PW during a random time. . . . . . . . . . . . . . 43 Figure 17 Individual P1 in the instrumented chair. . . . . . . . . . . . . . . . . . 47 Figure 18 Experimental results for individual P1. . . . . . . . . . . . . . . . . . 49 Figure 19 IGA comparison of simulation and real experiment for individual P1. . 50 Figure 20 Individual P2 in the instrumented chair. . . . . . . . . . . . . . . . . . 51 Figure 21 Experimental results for individual P2. . . . . . . . . . . . . . . . . . 52 Figure 22 Experimental results for subject H1 on sessions one and two respectively. 54 Figure 23 Experimental results for subject H1 on sessions three, four and five respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 24 IGA comparison of simulation and real experiment for individual H1. . 57 Figure 25 Experimental results for subject H2 on sessions one and two respectively. 59 Figure 26 Experimental results for subject H2 on sessions three, four and five respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Figure 27 Experimental results for subject H3 on sessions one, two and three respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Figure 28 Experimental results for subject H3 on sessions four and five respectively. 63 Figure 29 IGA comparison of simulation and real experiment for individual H3. . 64 Figure 30 Experimental results for subject H4 on sessions one and two respectively. 66 Figure 31 IGA comparison of simulation and real experiment for individual H4. . 67 Figure 32 Experimental results for subject H5 on sessions one and two respectively. 69 Figure 33 IGA comparison of simulation and real experiment for individual H5. . 70 Figure 34 Experimental results for subject H6 on sessions one, two and three respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Figure 35 Experimental results for subject H6 on sessions four and five respectively. 73 Figure 36 IGA comparison of simulation and real experiment for individual H6. . 74 Figure 37 Experimental results for subject H7 on sessions one and two respectively. 75 Figure 38 IGA comparison of simulation and real experiment for individual H7. . 76 Figure 39 Comparison of pulse width in both empirical and IGA tuning ap- proaches for H1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Figure 40 Comparison of pulse width in both empirical and IGA tuning ap- proaches for H6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Figure 41 RMSE analysis for both empiric and IGA tuning approaches. . . . . . 80 Figure 42 Multilayer perceptron. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Figure 43 Recurrent neural network. . . . . . . . . . . . . . . . . . . . . . . . . 84 Figure 44 LSTM cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Figure 45 Identification results for individuals P1, P2 and H1 respectively. . . . . 89 Figure 46 Identification results for individuals H2, H3 and H4 respectively. . . . 90 Figure 47 Identification results for individuals H5, H6 and H7 respectively. . . . 91 LIST OF TABLES Table 1 RISE controller for lower limb tracking control. . . . . . . . . . . . . . 17 Table 2 Empiric ideal response metrics. . . . . . . . . . . . . . . . . . . . . . . 34 Table 3 IGA ideal response metrics. . . . . . . . . . . . . . . . . . . . . . . . . 34 Table 4 IGA smooth response metrics response. . . . . . . . . . . . . . . . . . 34 Table 5 IGA moderate response metrics. . . . . . . . . . . . . . . . . . . . . . 35 Table 6 IGA critical response metrics. . . . . . . . . . . . . . . . . . . . . . . 35 Table 7 Specific data on analyzed individuals. . . . . . . . . . . . . . . . . . . 41 Table 8 Example of how datasets are encoded. . . . . . . . . . . . . . . . . . . 44 Table 9 Technical information on experiments for individual P1. . . . . . . . . 47 Table 10 Identification results for individual P1. . . . . . . . . . . . . . . . . . . 48 Table 11 Metrics on experimental results for individual P1. . . . . . . . . . . . . 48 Table 12 Technical information on experiments for individual P2. . . . . . . . . 50 Table 13 Identification results for individual P2. . . . . . . . . . . . . . . . . . . 51 Table 14 Metrics on experimental results for individual P2. . . . . . . . . . . . . 51 Table 15 Technical information on experiments for individual H1. . . . . . . . . 53 Table 16 Identification results for individual H1. . . . . . . . . . . . . . . . . . . 53 Table 17 Metrics on experimental results for individual H1. . . . . . . . . . . . . 54 Table 18 Technical information on experiments for individual H2. . . . . . . . . 58 Table 19 Identification results for individual H2. . . . . . . . . . . . . . . . . . . 58 Table 20 Metrics on experimental results for individual H2. . . . . . . . . . . . . 58 Table 21 Technical information on experiments for individual H3. . . . . . . . . 61 Table 22 Identification results for individual H3. . . . . . . . . . . . . . . . . . . 61 Table 23 Metrics on experimental results for individual H3. . . . . . . . . . . . . 61 Table 24 Technical information on experiments for individual H4. . . . . . . . . 65 Table 25 Identification results for individual H4. . . . . . . . . . . . . . . . . . . 65 Table 26 Metrics on experimental results for individual H4. . . . . . . . . . . . . 65 Table 27 Technical information on experiments for individual H5. . . . . . . . . 68 Table 28 Identification results for individual H5. . . . . . . . . . . . . . . . . . . 68 Table 29 Metrics on experimental results for individual H5. . . . . . . . . . . . . 68 Table 30 Technical information on experiments for individual H6. . . . . . . . . 71 Table 31 Identification results for individual H6. . . . . . . . . . . . . . . . . . . 71 Table 32 Metrics on experimental results for individual H6. . . . . . . . . . . . . 71 Table 33 Technical information on experiments for individual H7. . . . . . . . . 76 Table 34 Identification results for individual H7. . . . . . . . . . . . . . . . . . . 76 Table 35 Metrics on experimental results for individual H7. . . . . . . . . . . . . 76 Table 36 Identification results for all individuals (P1-P2, H1-H7). . . . . . . . . . 87 LIST OF ABBREVIATIONS & ACRONYMS CGA Complete Genetic Algorithm FES Functional Electrical Stimulation FGA Fast Genetic Algorithm GA Genetic Algorithm GRASP Greedy Randomized Adaptive Search Procedure IGA Improved Genetic Algorithm LSTM Long Short-Term Memory MLP Multilayer Perceptron NMES Neuromuscular Electrical Stimulation NN Neural Network PW Pulse Width RIP Real Initial Population RISE Robust Integral of the Sign of the Error RMSE Root Mean Squared Error RNN Recurrent Neural Network SCI Spinal Cord Injury LIST OF SYMBOLS R+ Set of Real numbers U Input space Y Output space F Functional space ⊗ Element-wise multiplication tanh Hyperbolic tangent sgn(·) Signum function TABLE OF CONTENTS 1 INTRODUCTION 15 1.1 CONTEXT OF THE PROBLEM 15 1.2 MOTIVATIONS 16 1.3 OBJECTIVES AND HYPOTHESES 17 1.4 MASTER THESIS OUTLINE 18 2 PROPOSED METHODOLOGY AND THEORETICAL BACKGROUND 20 2.1 RISE CONTROL DEVELOPMENT 21 2.2 IMPROVED GENETIC ALGORITHM 22 2.3 SYSTEM IDENTIFICATION VIA NEURAL NETWORKS 27 3 SIMULATION RESULTS WITH A MATHEMATICAL MUSCULAR MODEL 30 3.1 HUMAN LOWER LIMB MATHEMATICAL MODEL 30 3.2 MATERIAL AND METHODS 31 3.3 RESULTS AND DISCUSSION 33 4 EXPERIMENTAL RESULTS WITH NEURAL NETWORK BLACK-BOX MODELS 40 4.1 INSTRUMENTATION 40 4.2 ANALYZED INDIVIDUALS 40 4.3 EXPERIMENTAL SETUP 42 4.3.1 First session 43 4.3.2 Two up to five sessions 45 4.4 RESULTS AND DISCUSSION 45 4.4.1 Individual P1 46 4.4.2 Individual P2 49 4.4.3 Individual H1 52 4.4.4 Individual H2 57 4.4.5 Individual H3 61 4.4.6 Individual H4 65 4.4.7 Individual H5 67 4.4.8 Individual H6 70 4.4.9 Individual H7 74 4.5 CONCLUSION 77 5 NONLINEAR IDENTIFICATION OF THE KNEE JOINT ANGULAR PO- SITION UNDER NMES/FES APPLICATION VIA DEEP AND DYNAMIC NEURAL NETWORKS 81 5.1 NEURAL NETWORK METHODS 82 5.1.1 Multilayer Perceptron (MLP) 82 5.1.2 Recurrent neural network (RNN) 83 5.1.3 Long Short-Term Memory (LSTM) 83 5.2 MODEL SELECTION 85 5.3 DATA ENCODING 86 5.4 RESULTS AND DISCUSSION 86 6 GENERAL CONCLUSIONS 92 6.1 FUTURE WORKS 93 6.2 PUBLICATIONS 93 REFERENCES 95 15 1 INTRODUCTION 1.1 CONTEXT OF THE PROBLEM In the last few years, several studies have been carried out showing that neuromuscular electrical stimulation (NMES) and functional electrical stimulation (FES) can produce good results in rehabilitation treatments of spinal cord injured patients. On the one hand, damages to the spinal cord may be occasioned by traumatic causes like a road accident, sports injuries, and violence, or, nontraumatic ones such as diseases that destroy the neurological tissues and tumors. Notwithstanding, the problem itself, spinal cord injury (SCI), is often irreversible and it can cause some issues like the inability to complete daily activities or occupational ones; loss of sensation and function of muscles; difficulties related to sexual functions; total or partial paralysis; and a lot of pain to the individual. However, the main consequences depend on several factors, such as the personal condition of the patient itself, the extent of the damage, availability of time and resources, and socioeconomic factors, e.g., in low-income countries the SCI normally lead to death in contrast with high-income ones where patients enjoy a better and more productive life. The application of NMES/FES for SCI rehabilitation is one of the most frequent methods, which provides many health and social benefits for its patients, such as the maintenance and re- covery of muscles strength; prevention of flaccidity and hypotrophy, which are classic signs of muscle inactivity; offering higher expectation and quality of life; and by allowing social reinser- tion. NMES/FES are techniques based on the use of equipment that generates electrical current for muscle stimulation at the motor level, aiming to generate a muscle contraction via surface or intramuscular electrodes. On the use of electrical stimulation, the strength of muscle con- traction is controlled by manipulating pulse frequency, amplitude or duration, given stimulation being delivered as a waveform of electrical current pulses (HO et al., 2014; LAW; SHIELDS, 2007; LYNCH; POPOVIC, 2008). Efforts have been made motivated by the promising therapeutic treatment and beneficial results of NMES to increase the efficacy of motor rehabilitation on accomplishing functional tasks (e.g., grasping, walking, standing, cycling and sitting pivot transfers) (KAPADIA; ZI- VANOVIC; POPOVIC, 2013; JAIME; MATJACIC; HUNT, 2002; FONSECA et al., 2017) where it is named as FES. However, there are numerous challenges to be faced in designing automatic stimulation strategies, e.g., the system must run in real-time and safely even in face of bodily uncertainties such as the muscle response to the electrical stimulus, pH, body fat, hydration level, body temperature, inter and intra-subject variability in muscle properties. A 1.2 MOTIVATIONS 16 particular problem of variability in muscle properties is the difficulty of predicting the exact contraction force exerted by the muscle, which generates an unknown mapping between the stimulus parameters and the generated muscle force. To overcome some of these problems, bet- ter control procedures using feedback and/or feedforward techniques has been developed (FER- RARIN et al., 2001; PREVIDI; CARPANZANO, 2003; JEZERNIK; WASSINK; KELLER, 2004; LYNCH; GRAHAM; POPOVIC, 2011; NUNES et al., 2018; LYNCH; POPOVIC, 2012; SANCHES et al., 2014; SHARMA et al., 2009; OLIVEIRA et al., 2017; CHENG et al., 2016; STEGATH et al., 2007, 2008) to benefit of acquired information about the controlled system, i.e., the movements measured in real-time with sensors, are used by the control system to adjust the control operation in a closed-loop design. 1.2 MOTIVATIONS Even though there are several investigations on closed-loop control of FES/NMES for lower limb rehabilitation, these systems are hardly put into production. Alternatively, there exist com- mercial stimulators normally available on open-loop designs and with pre-programmed electri- cal stimulus, which are not enough to deal with the nonlinear and time-varying nature of mus- cles. Thus, the main reason to continue investigating on this field is that real-world NMES/FES applications to rehabilitate SCI patients require control strategies that compensate for model- ing errors on the plant, system’s faults, individual’s muscles behavior, external disturbances, nonideal muscle conditions (fatigue, tremor, and spasms) and many other factors (LYNCH; GRAHAM; POPOVIC, 2011; MOHAMMED et al., 2012; LEW et al., 2016; PECKHAM; KNUTSON, 2005). The present study investigated the continuous and robust control technique for uncertain nonlinear systems that have been reported in the literature as robust integral of the sign of the error (RISE) (XIAN et al., 2003; XIAN; QUEIROZ; DAWSON, 2004). In the literature one can find Stegath et al. (2007, 2008) and Sharma et al. (2009) as pioneers authors on develop- ing the RISE controller for the lower limb tracking control. Afterward, Sharma et al. (2012) presented an improvement of the RISE control method for the same application using an NN feedforward term. Downey, Cheng and Dixon (2013) and Downey et al. (2015b), developed a RISE controller for asynchronous stimulation to the lower limb. Kawai et al. (2014) simulated the tracking control performance of a RISE-based controller to a model to the co-contraction control of the human lower limb. Lastly, Kushima et al. (2015) modeled an FES knee bending and stretching system and developed a RISE-based controller to stimulate the quadriceps and hamstrings muscle groups. The RISE method was chosen to perform this study by some intrinsic characteristics, such as it considers disturbances that were not previously modeled and it also has implicit learn- ing characteristics, which are important in performing rehabilitation experiments. However, 1.3 OBJECTIVES AND HYPOTHESES 17 adjusting the parameters of the controller is the main component to guarantee high-quality per- formance, i.e., the method can only guarantee good responses selecting appropriately the gain constants. Considering the stability analysis and the gains sufficient condition initially provided by Sharma et al. (2009) for an uncertain nonlinear muscle model, it is difficult to reach ex- pected control effects in practice just adjusting empirically the gain parameters in the extremely large search space R+. In daily routines of NMES/FES application to SCI patient rehabilitation exist several problems like muscle fatigue, due to incomplete tetanus and even from the electri- cal stimulus application itself, which would be increased by applying a “trial and error” method to tune the controller. Although the literature indicates good lower limb tracking performance using the RISE controller, the motivation of this research emerges from the lack of intelligent techniques to adjust the RISE controller parameters, wherein the aforementioned researches authors did not inform how they tuned the controller or that was used an empiric approach (pretrial tests) to define gain parameters before conducting the real experiments. Besides, those controllers were tested only on healthy patients, and as well-known, muscles of SCI patients do not have the same strength as a healthy one, and the brain may have lost the regular way of communication given certain damage in the spinal cord tissue. Table 1 presents information on the number of patients and the tuning procedure used by the aforementioned researches. Table 1 - RISE controller for lower limb tracking control. Authors and years Validation Tuning Stegath et al. (2007, 2008) 2 healthy subjects Not informed Sharma et al. (2009, 2012) 5 and 9 healthy subjects Not informed Kawai et al. (2014) Simulation Adjusted by simulation Kushima et al. (2015) 7 healthy subjects Not informed Downey et al. (2015b) 4 healthy subjects Pretrial tests Source: Developed by the author 1.3 OBJECTIVES AND HYPOTHESES To overcome the aforementioned problems, this paper proposes a novel robust and intelli- gent control-based methodology to human lower limb rehabilitation via NMES/FES, aiming to overcome the empiric tuning procedure to the RISE controller observed in the literature. More- over, this study proposes to extend the analysis of the RISE controller to paraplegic individuals who present not ideal conditions as healthy ones. This methodology includes an identification step based on neural networks black-box models where the novelty is the use of past identi- fication and control data for each patient, the RISE control method (or fundamentally similar 1.4 MASTER THESIS OUTLINE 18 control laws) to guarantee the system’s stability and an offline controller optimizer. Specific objectives of this research are: • Propose an improved genetic algorithm (IGA) as an offline RISE controller optimizer to the problem of nonlinear lower limb tracking via NMES/FES; • Validate the proposed methodology with simulation results, considering healthy and para- plegic identified subjects from Ferrarin and Pedotti (2000) on ideal and nonideal muscle conditions (including fatigue, spasms, and tremor) using a mathematical muscle model; • Substantiate the proposed methodology with experimental results, considering healthy and paraplegic subjects identified via neural networks black-box models; • Analyze the effectiveness of using past identification and control data on the construction of NNs black-box models for each subject; • Investigate the performance of deep and dynamic neural networks namely the multilayer perceptron, a simple recurrent neural network, and the Long Short-Term Memory archi- tecture in the specific task of identifying the knee joint angular position under electrical stimulation. On the one hand, the first hypothesis for this research is that by using an empirical approach to clinical procedures, it would present a large number of poor performances, while a more adequate tuning with a more representative identified model can provide better tracking control of the lower limb. And, the second hypothesis is that by using past rehabilitation data for identifying a patient, this model will improve the description of the relationship between angular position and the delivered electrical stimulus, where fatigue and other problems as tremors are already implicit in the data. 1.4 MASTER THESIS OUTLINE The following chapters are structured as follow. • Chapter 2 introduces the proposed control-based methodology for the human lower limb rehabilitation with NMES/FES and its theoretical background. It proposes the improved genetic algorithm (IGA), which aims to better-acquire the RISE controller parameters. It provides the basic principles of the RISE control design. Lastly, it describes system identification and modeling via neural networks. 1.4 MASTER THESIS OUTLINE 19 • Chapter 3 presents simulation results for five healthy and three paraplegic subjects iden- tified by Ferrarin and Pedotti (2000) using a mathematical model. It mathematically de- scribes the dynamic lower limb model and the reproduction of the “non-idealities” block (modeling of fatigue, tremors, and spams) from Lynch, Graham and Popovic (2011). It defines the metrics used to evaluate the control performance and presents the materials and methods. Simulation results are presented and discussed to isotonic and isometric contractions on ideal and nonideal conditions using an empiric and an IGA tuning ap- proaches. • Chapter 4 presents experimental results made with seven healthy subjects and two para- plegic ones. It describes the instrumentation apparatus, analyzed individuals, and the experimental procedure used in stimulation sessions. It defines the metrics used to eval- uate the control and identification performances. Results are presented and discussed to both isotonic and isometric contractions using an empiric and an IGA tuning approach. • Chapter 5 introduces the investigation of deep and dynamic neural networks for identify- ing the knee angular position under NMES/FES. It briefly describes the MLP, a simple RNN, and the LSTM architecture. It presents how data were structured and the selec- tion of each model and its hyperparameters. Results are presented using data acquired in Chapter 4 comparing the effectiveness of each NN model. • Finally, in Chapter 6 the general conclusions of this research and future works are pre- sented. 20 2 PROPOSED METHODOLOGY AND THEORETICAL BACKGROUND In the literature, it is shown that an empirical tuning approach to the RISE controller can also guarantee stability if gains inequalities are respected given the stability proof analysis pro- posed by Sharma et al. (2009). However, once gains selections are immense it is likely to one choose combinations that would not guarantee the best performance and accuracy during the rehabilitation of SCI patients. Hence, Figure 1 illustrates the proposed methodology, where at the first session a patient has no data and a stimulation session is made on the attempt to acquire information on the relationship between delivered electric stimulus and reached angular position. The acquired data is appropriately treated to pass through an identification step via neural network black-box models. Therefore, assuming that this relationship is efficiently mapped, a simulation process making use of a global and fast search optimization algorithm that tests numerous gain com- binations evaluating a well-defined control task is made to efficiently tune the RISE controller for this patient. At final, the rehabilitation procedure is retaken with fine-tuned gains for a bet- ter control-stimulation session, aiming to prevent premature fatigue and other problems of SCI patients that would be present by not choosing an appropriate gains combination. In future sessions, all data (control and identification) from previous rehabilitation sessions are used for training a neural network model in an offline scheme, where before each next sec- tion, all data from a patient is combined to a single dataset and used to better map its relationship with electrical stimulus. Thus, the same optimization process using the trained model provides fine-tuned gain parameters to be afterward applied to the rehabilitation procedure. In Section 4.3 the methodology is experimentally implemented, where it is explained more comprehensively. The use of neural networks is motivated by the advantages of these methods for the non- linear system identification problem and by a high power for computation and storage of data encountered nowadays. However, as one can see in Chapter 3 and its preliminary results in Ar- colezi et al. (2019), our proposal is not limited to black-box modeling, where a mathematical model and identified parameters of five subjects and three SCI patients are used to provide simulation results. To the identification step, the novelty of the proposal is motivated by the use of past reha- bilitation data. The primary purpose is to build up a dataset for each patient, where the number of data will increase during rehabilitation sessions, and the identified model will improve with more data and details about the nonlinear muscular behavior. As highlighted in the literature, muscular behavior is very susceptible to parametric variation between one day to another, and 2.1 RISE CONTROL DEVELOPMENT 21 Figure 1 - The proposed robust and intelligent control-based methodology. Source: Developed by the author for instance, evolution and gain of strength due to previous rehabilitation sessions. Moreover, one of the primary advantages of performing simulations for an NMES-base knee extension is the liberty of studying this problem from different perspectives and divergent levels of abstraction with the acquired data. Different to quite of human limitations to NMES ap- plications, that restrict the number of experiments, simulation provides uncountable executions to better study feasibility and practicality of the designed system, and it supplies continue feed- back to maintain improving the system (JEZERNIK; WASSINK; KELLER, 2004; DURIEZ; BRUNTON; NOACK, 2017; LYNCH; GRAHAM; POPOVIC, 2011). The following sections of this chapter describe each component of the methodology, the RISE control method, an improved genetic algorithm (IGA) for the optimization procedure, and system identification via neural networks black-box models. 2.1 RISE CONTROL DEVELOPMENT The RISE control technique has been proposed as a continuous-time and high gain feedback control approach for uncertain nonlinear systems, which even in spite of bounded smooth exter- nal disturbances and bounded modeling uncertainties, the control law can guarantee asymptotic 2.2 IMPROVED GENETIC ALGORITHM 22 tracking (XIAN et al., 2003; XIAN; QUEIROZ; DAWSON, 2004). To achieve the stated control objective, i.e., to enable the lower limb to track a desired angular trajectory even in spite of ex- ternal disturbances and modeling uncertainties, a position tracking error denoted by e1(t) ∈ R, is defined as e1(t) = θd(t)−θ(t), (1) where θd(t) is the desired angular trajectory assumed to have bounded continuous time deriva- tives, and θ(t) the actual position. Additionally, filtered tracking errors e2(t) ∈ R and r(t) ∈ R are defined as e2(t) = ė1(t)+α1e1(t), (2) r(t) = ė2(t)+α2e2(t), (3) where α1,α2 ∈R denote positive and adjustable control gains. Authors in Sharma et al. (2009), Stegath et al. (2008), proved semi-global asymptotic stability for an uncertain nonlinear muscle model with the RISE control law defined as u(t) = (ks +1)e2(t)− (ks +1)e2(0)+ ∫ [(ks +1)α2e2(τ)+β sgn(e2(τ))]dτ , (4) where ks,β ∈ R are also positive and adjustable control gains, and sgn(·) denotes the standard signum function. As highlighted in the literature, a recommended differentiator for simulations, and used in this work, is mathematically described as H(s) = Y (s) U(s) = s τs+1 , (5) being composed by the low pass filter 1/(τs+1); where Y (s) and U(s) are the input signal and its filtered derivative, respectively; and τ is a time constant. 2.2 IMPROVED GENETIC ALGORITHM Genetic algorithms (GAs) are global search algorithms, which solves mathematical opti- mization problems based on mechanisms of natural selection. GAs are a class of evolutionary algorithms used to identify and optimize parameters of input-output relationships, in which evo- lutionary algorithms are an important category of machine learning techniques for optimization. Every generation of a GA, a set of individuals named population, compete at a well-defined 2.2 IMPROVED GENETIC ALGORITHM 23 task, and a new set of individuals are created using certain genetic operations. Generally, the initial population of a GA is randomly generated, and every individual is evaluated concerning an objective function. Individuals with lower costs (minimization problem) are more likely to advance to the next generation when passing by a selection method. Afterward, the whole pro- cedure evolves certain genetic operations, such as elitism where the best individuals advance directly to the next generation; replication where individuals also advance directly to the next generation, but with a probability rate related to the cost function; crossover, in which two in- dividuals are selected to interchange parameters; and mutation where a small random change is made in an individual. Fig 2 illustrates an iteration of a standard GA with the respective genetic operations (DURIEZ; BRUNTON; NOACK, 2017; FLEMING; PURSHOUSE, 2002). Figure 2 - Genetic operations to advance one generation to the next. Source: (DURIEZ; BRUNTON; NOACK, 2017) In this research an improved genetic algorithm (IGA) is proposed to better optimize gains parameters of the RISE controller, which is partially inspired on the standard methodology of the greedy randomized adaptive search procedure (GRASP) proposed by Feo and Resende (1995), running with GA instead of using hybridization of a simple heuristic algorithm and a simple local search procedure. We took advantage of the basic idea of GRASP, which runs in a multistart framework and has the subsequent steps. Firstly, it is implemented a preprocessing step to identify convergence patterns and initiate the process efficiently bounding gain limits. In the construction phase, we used a simple fast genetic algorithm (FGA) to generate a good set of solutions named the real initial population (RIP) to finally run a local search procedure based on a complete genetic algorithm (CGA). Both fast and complete GAs are very similar and related to general approaches of GAs they are different in some aspects as encoding, mutation algorithm to fit our problem and replacement of individuals in the current population. Description of the whole improved genetic algorithm procedures is detailed below. 2.2 IMPROVED GENETIC ALGORITHM 24 • A solution to our problem is a chromosome consisting of four real and positive num- bers representing the gains parameters of the RISE controller. As shown in Fig 3, every solution representing α1,α2,ks,β ∈ R+ respectively, concerning gain conditions stated in Sharma et al. (2009) is valid. Figure 3 - Real and positive coded genes of a chromosome. Source: Developed by the author • At this stage, a selection procedure is made to determine the best chromosomes, according to fitness values. The tournament selection was chosen for many benefits that it provides, such as it is easy to code, efficient and fast (FLEMING; PURSHOUSE, 2002). The idea of this process is described below and illustrated in Figure 4. i. Provide n (size of population) games, in which j (2 for this study) individuals are chosen randomly from the present population and the individual with best fitness win each game; ii. Finally, the two best individuals, i.e., that won more games during the tournament, are chosen as parents. Figure 4 - Tournament selection. Source: Developed by the author • The two individuals chosen in the tournament selection process as parents are submitted to the crossover algorithm. This operator generally combines features of two parents to 2.2 IMPROVED GENETIC ALGORITHM 25 form two offspring. In this research, the commonly single-point crossover operator is used in the middle of the chromosome. • The mutation algorithm is designed to expand the search space to regions that may not be close to the actual population, i.e., to increase diversity and avoid premature conver- gence (DURIEZ; BRUNTON; NOACK, 2017; FLEMING; PURSHOUSE, 2002). It has the following characteristic: Select one point of the chromosome randomly each stage (concerning to a mutation rate), and add or subtract small or medium value to this point according to the preprocessing step, and take its absolute value (e.g., α1,α2 accept small variations and ks,β accept medium variations). • Further, as the purpose of our nonlinear lower limb tracking problem via NMES/FES is to optimally stimulate the knee angle to track a desired reference, a minimization problem is defined as min : J(α1,α2,ks,β ) = RMSE + penalty, (6) RMSE = ∫ T 0 √ E((θd −θ)2), (7) penalty = ∫ T R 0 √ E((θd −θ)2). (8) where T is the whole period and TR is the transient response. In other words, the main objective is to minimize the metric root mean squared error (RMSE) of the actual and desired knee angle, penalizing poor transient response aiming to obtain fast responses with low overshoot. • Construction phase: the RIP set is generated in an efficient and simple way using the FGA, described in the following steps i. Preprocess the problem to identify convergence patterns; ii. Define the following initial parameters: size of population (small), number of gen- erations, rate of crossover and mutation operators; iii. Generate the initial population randomly and evaluate the fitness function; iv. Select two parents using tournament selection; v. Perform the single-point crossover and the mutation algorithm; vi. Substitute only these two offspring in place of the two worst individual in the current population; 2.2 IMPROVED GENETIC ALGORITHM 26 vii. Evaluate and select the best individual of this small population to be a member of the RIP; viii. Repeat steps (iii) to (vii) until the algorithm reach the predefined number of gener- ations. Note that, the number of generations is equal to the RIP set size and when the FGA stops, provides the RIP set of individuals also evaluated about their fitness function. Figure 5 shows a generation of the FGA, where the best solution of this new population is selected to be part of the RIP. Figure 5 - Fast genetic algorithm operators to perform one generation and select the best individual. Source: Developed by the author • Local search phase: the local search procedure is performed with the CGA describe below i. Evaluate if there are any repeated individuals in the RIP, and if positive excludes them from population; ii. Define the following initial parameters: number of generations, rate of crossover and mutation operators; iii. Select two parents using tournament selection; iv. Perform the single-point crossover and the mutation algorithm; v. Evaluate if the two new individuals are able to replace others from population. The individual is accepted if, and only if it is not equal to another one in the population, and if it performs better than the worst individual in the current population; vi. Repeat steps (iii) to (v) until the algorithm reach the predefined number of genera- tions. Note that the CGA is very similar to the FGA, and when the algorithm reaches the pre- defined number of generations, it is expected to have a whole final set of different great 2.3 SYSTEM IDENTIFICATION VIA NEURAL NETWORKS 27 solutions. If one chooses a multistart framework of k iterations, the whole process starts again with the FGA to construct the RIP set, passing to the local search with CGA until complete the predetermined number of iterations, presenting the best final set as solu- tion. One is instructed to balance between the multistart framework and the total time of execution. 2.3 SYSTEM IDENTIFICATION VIA NEURAL NETWORKS Nonlinear systems identification and modeling have been applied in most areas of science to predict the future behavior of dynamic systems. It has been fomented fields in control theory, and it is an essentially important way to explore, study and understand the world by a formal description of events as a model. On the one hand, there are mathematical models, which are based on the direct description of systems using its physics and mechanics resulting in linear or nonlinear differential equations depending on several parameters. Notwithstanding, there are black-box models, which use no physical insight about the system and will be the approach of this research via NNs. Generally, black-box models are simpler than physical modeling and more appropriate where there is a lack of knowledge of the underlying physiology, or in the case when the physical knowledge is too complex, as it is in this case, where muscular behavior presents time-varying and high nonlinearities properties. The construction of black-box models is essentially based on the quality of measured data about the system, where the collection of as much relevant data as possible, which cover the whole operating range of interest, is considered the bottom line to achieve success with the developed models. The fundamental concept of this approach is to model the direct input-output relationship, i.e., identifying and modeling just with data, in which the main objective is to find the weights and other coefficients (known as hyperparameters) of the NN (ABLAMEYKO, 2003; WANG, 2017; GONZALEZ; YU, 2018; HAYKIN, 2009). Moreover, NNs are based on a collection of inter-connected units named neurons. These neurons are structured into three or more layers, input, hidden(s), and output. NNs are in the core of deep learning (several neurons and hidden layers) field and has become a progressively and very popular research topic. Generally, NNs can be divided into two large classes, feedfor- ward NNs and recurrent neural networks (RNNs), where the first has been used for identification and control of simple nonlinear systems and the later for sequences modeling and time-series forecasting (OGUNMOLU et al., 2016) Fundamentally, an operator F from an input space U to an output space Y expresses the model of the system to be identified, where the goal is to find a function F̂ that approximates F to a specific requirement. By the Stone-Weierstrass theorem, there exists a continuous and 2.3 SYSTEM IDENTIFICATION VIA NEURAL NETWORKS 28 bounded function F , that can be uniformly approximated as closely as desired by a polyno- mial function F̂ . Further, according to the universal approximation theorem, there exists a combination of hyperparameters of an NN that allows it to identify and learn any nonlinear function (HAYKIN, 2009; NARENDRA; PARTHASARATHY, 1990; WANG, 2017). In the following, a short description of using NNs for identification of discrete dynamic system is provided. Thus, consider a single-input and single-output discrete system structure with only the input and output data available as y(k) = f [y(k−1), ...,y(k−n);u(k−1), ...,u(k−m)], (9) where f (·) is an unknown nonlinear difference equation that represents the plant dynamics; u and y are measurable scalar input and output respectively; and m and n are the maximum lags for the system output and input, i.e., they are the last values of the input and output respectively. In short, the next value of the dependent output signal y(k) is regressed on previous values of the output and input signals. Extended versions for multi-input and multi-output cases are also possible, but it is not in the scope of this research. The identification for the discrete-time system in Eq. 9 can be done by following two major types of identification structures presented in the literature as the parallel, and the series-parallel identification model (WANG, 2017; GONZALEZ; YU, 2018; NARENDRA; PARTHASARATHY, 1990). The first structure is formed using past inputs of the system and the outputs of the NN model, which is a more stringent test, and it is only recommended for sta- ble plants. In the second arrangement, both past inputs and system outputs are used (the actual output of the system rather than that of the identification model), which is always recommended for stability reason and motivated by data availability of past system output. Mathematically these models, parallel, and series-parallel are described respectively as ŷ(k) = F [ŷ(k−1), ..., ŷ(k−n);u(k−1), ...,u(k−m)], (10) ŷ(k) = F [y(k−1), ...,y(k−n);u(k−1), ...,u(k−m)], (11) where ŷ is the model output, y is the real system output, F is the model structure, m and n are the regression order for the input and output. These last two parameters are chosen before the identification process, which n is the output memory to indicate how many past steps of output will be used in the system identification, and m refers to time-step of input values and it is the longest memory that a model can store. Both parallel and series-parallel identification schemes are shown in Figures 6 and 7 respectively. 2.3 SYSTEM IDENTIFICATION VIA NEURAL NETWORKS 29 Figure 6 - Parallel identification scheme of nonlinear dynamic system. Source: Adapted from Wang (2017) Figure 7 - Series-parallel identification model of nonlinear dynamic system. Source: Adapted from Wang (2017) 30 3 SIMULATION RESULTS WITH A MATHEMATICAL MUSCULAR MODEL On the attempt to validate the novel approach for lower limb rehabilitation via NMES/FES, simulations results are provided using a mathematical muscle model and parameters for three SCI subjects and five healthy identified from Ferrarin and Pedotti (2000). Besides the new ap- proach, another novelty is reproducing non-idealities (fatigue, spasms, and tremor) encountered on the SCI population based on the work from Lynch, Graham and Popovic (2011) to provide a more realistic scenario for simulating real-world experiments. The following sections of this chapter describe the mathematical muscle model; materials and methods (developed platform, metrics, methodology); results and its discussion. 3.1 HUMAN LOWER LIMB MATHEMATICAL MODEL Stegath et al. (2007, 2008) and Sharma et al. (2009) developed an uncertain nonlinear muscle model to design and implement the RISE controller to yield an asymptotic stability result even in spite of uncertain nonlinear muscle model and in the presence of additive bounded disturbances. Nevertheless, the muscle model described in this chapter was based on Ferrarin and Pedotti (2000) and on Lynch, Graham and Popovic (2011), which includes muscle fatigue, spasms and tremors being the control of the knee joint made by varying the width of the current pulse. The lower limb model considers that the individual is seated and the lower leg is freely suspended, as illustrated in Figure 8. Knee flexion is provided by the gravity force and knee ex- tension is given by external electrical stimulus. The knee angular acceleration θ̈(t) is expressed as θ̈(t) = 1 J (−mgl sinθ(t)− τsti f ness −Bθ̇(t)+υ(t)), (12) where J is the moment of inertia of combined shank and foot, θ(t), θ̇(t), θ̈(t) are the knee angular position, velocity and acceleration respectively, B is the viscous damping coefficient, m is the combined mass of shank and foot, g is the gravity acceleration, and l is the distance between the knee joint and center of mass of shank and foot. The stiffness torque is τsti f f ness = λe−E(θ+ π 2 )(θ + π 2 −ω), (13) 3.2 MATERIAL AND METHODS 31 where λ , E are exponential term coefficients and ω is the resting elastic knee angle. Addition- ally, this research simulated muscle fatigue, tremors, and spasms that modify muscle nominal torque in three levels (smooth, moderate and critical). For more information, Lynch, Graham and Popovic (2011) is highly recommended reading and consulting. The mathematical model of this block is υ(t) = (1+ spm(t)+ tr(t))τquad(t) f at(t), (14) where υ(t) is the output of the “non-idealities” block, spm(t) is a scaled torque resulted from a muscle spasm, tr(t) represents a scaled torque resulted from a muscle tremor, f at(t) indicates a fatigue waveform and τquad(t) is the input to the “non-idealities” block, which relates the stimulated Pulse Width (PW) delivered to quadriceps muscle and resulting torque exerted about the knee joint, expressed in the frequency domain as τquad(s) = G 1+ηs PWquad(s), (15) where PWquad(s) is the pulse width and G, η are constants of muscle activation function. Figure 8 - Schematic representation of the lower limb with surface electrodes. Source: Adapted from Ferrarin and Pedotti (2000) 3.2 MATERIAL AND METHODS Initially, a Matlab/Simulink R⃝ system model that maps stimulation PW and quadriceps torque, using the nonlinear second-order dynamics of the knee and lower leg described in Section 3.1, was developed. Additionally, to perform a more realistic simulation process in- cluding muscle fatigue, tremors and spasms, the “non-idealities” block from Lynch, Graham and Popovic (2011) was reproduced and implemented. Figure 9 illustrates waveforms of all non-idealities included at the same time t = 15s as smooth, moderate and critical to the normal- ized nominal torque. A saturation block is attached to the system model to comply with real 3.2 MATERIAL AND METHODS 32 applications of NMES to individuals, bounding control signal from 0 µs to 500 µs. Figure 9 - Normalized nominal torque with smooth, moderate and critical non-idealities. Source: Developed by the author In this research, the performance of the RISE controller was evaluated considering both time transient and stationary responses, owing to the fact that it will be future tested on real subjects, is important to take into account information about transient response. Three paraplegics (P1- P3) and five healthy (H1-H5) identified individuals from Ferrarin and Pedotti (2000), were used to simulations considering two references trajectories. The first one is a sinusoidal trajectory ranging from 12.5◦ to 45◦ to simulate an isotonic contraction with a repetitive pattern that starts on 20◦, and the second one is a 45◦ step trajectory replicating an isometric contraction. Initial parameters of IGA used in simulations were size of population equal 10, size of RIP equal to 50, crossover rate equal 1, mutation rate equal 0.3, number of generations equal 50 and k equal 1 iteration. Algorithm performed 30 trials given probabilistic properties of GA to ensure convergence pattern to global minimal. Initially, simulations will try to tune the controller to the worst case scenario, i.e., with all non-idealities (fatigue, tremors and spasms) included to the model as critical to find the best gain parameters. Afterwards, with these ‘best’ parameters it will also be tested to the other cases, which combines non-idealities as smooth, moderate and the ideal case with none of them. In addition, simulation results for empirical tuning are presented to compare effectiveness, and all results are presented without penalizing transient response as it is during the optimization process. In view that SCI population generally presents diagnostics with some aforementioned non- idealities, it motivated us to follow this methodology, which would allow examining if and how the RISE controller could compensate in real-world experiments, in contrary to just ideal cases. Thus, our hypothesis is that by previously including critical external disturbances in simulation procedures, which performs several combinations of gains to better fit on a particular individual, will imply on a more robust controller prepared for extreme cases, and good performances in 3.3 RESULTS AND DISCUSSION 33 smoother situations will be less difficult to achieve. Further, different metrics for ideal and nonideal conditions were considered to both tracking and regulation problems. To the sinusoidal trajectory and with ideal conditions, metrics were the lag to the desired and actual knee angles lag(t); and the RMSE Ermse(deg) between the desired and actual knee angles. Considering nonideal conditions, these metrics were the exact time where the controller could not compensate tracking anymore for at least 5 seconds uncT (t), and the Erms. For the step trajectory along ideal conditions, these metrics were calculated as 10%90% rise time τrise(t); the steady-state error ess(deg) between the actual and desired knee angles; the percent overshoot M.O(%) past the steady-state knee angle; and the 2% settling time τsettling(t). For non-ideal conditions, these metrics were calculated as the time where the controller could not compensate regulation anymore for at least 5 seconds uncR, Erms, peak angular position PA (upper or lower the desired knee angle), and the time to peak angular position PT respectively. 3.3 RESULTS AND DISCUSSION Tables 2-6 present metrics performances for control systems, considering an empiric and an IGA tuning approaches for both regulation and tracking situations respectively. Further, Figure 10 presents results on ideal conditions using empiric and IGA tuning for all patients related to both trajectories. Similarly, Figures 11-13 present results for all patients in both trajectories under nonideal muscle conditions (fatigue, tremor, and spasms) as smooth, moderate and critical respectively, using the IGA tuning approach. Results for an empirical tuning approach to the RISE controller could also guarantee stabil- ity on ideal conditions, but once gains selections are immense, one likely choose combinations that would not guarantee best performance and accuracy during rehabilitation. For instance, see Figure 10 presenting commonly found results on our simulations, which by using empiric gains on clinical procedures, the rehabilitation of SCI patients would not result in good tracking performances. On the other hand, the RISE controller can provide high-quality performance with the right selection of gains, motivating the use of some intelligent tuning approach as the one proposed in this research. Considering IGA fine-tuned responses on ideal conditions, very good results for tracking control were achieved for paraplegic patients. However, results adding non-idealities to the model showed different muscular behavior. For smooth non-idealities included to the model, healthy patients well compensated the whole period with a low increment of RMSE; P1 pre- sented small variations and compensated much better than P2 and P3, where P2 in the final cycle and P3 in the two final cycles failed on tracking the reference angle. Similarly, in moderate and critical scenarios healthy patients well compensated the whole period with a low increase of 3.3 RESULTS AND DISCUSSION 34 RMSE; P1 full failed on compensating in the last cycle, and patients P2 and P3 full failed in last cycles. Further, results for the step signal trajectory considering ideal conditions lead to zero steady-state error with good performance. Similarly to tracking situations, responses for smooth non-idealities are well compensated by healthy patients and by P1, while P2 and P3 failed on compensating after 40 and 30 seconds respectively. For moderate and critical scenarios, healthy patients could still compensate for almost the whole period; P1 presented good compensation on moderate and poor one on the critical situation; P2 presented problems on compensating after 30 seconds approximately, and P3 failed compensation just after introducing the non-idealities at 15 seconds. Table 2 - Empiric ideal response metrics. Metric H1 H2 H3 H4 H5 P1 P2 P3 τrise(t) 0.6138 0.5962 0.3565 0.3294 0.3638 0.7995 3.2453 3.0694 ess(deg) 0 0 0 0 0 0 0 0 M.O(%) 16.687 16.902 24.256 25.799 26.232 24.380 0.0289 0.0206 τsettling(t) 6.2046 6.1488 6.7681 6.6522 6.7461 8.1102 6.5566 6.3417 lag(t) 5.4 5.5 5.3 5.5 5.4 6.0 0.8 0.6 Ermse(deg) 2.4956 2.7943 3.9087 3.4864 3.4943 2.8409 0.7007 0.6913 Source: Developed by the author Table 3 - IGA ideal response metrics. Metric H1 H2 H3 H4 H5 P1 P2 P3 τrise(t) 0.3352 0.58487 0.3349 0.5737 0.6076 0.7628 0.6209 0.2638 ess(deg) 0 0 0 0 0 0 0 0 M.O(%) 23.894 18.151 26.610 19.313 17.277 30.414 9.0326 11.678 τsettling(t) 4.97 5.8259 6.6995 7.5461 8.7923 7.4267 6.2543 4.9745 lag(t) 4.87 5 4.5 4.69 4.65 5.3 0.55 1.2 Ermse(deg) 2.4413 2.5868 3.7808 3.5304 3.3886 3.2284 0.6755 0.7308 Source: Developed by the author Table 4 - IGA smooth response metrics response. Metric H1 H2 H3 H4 H5 P1 P2 P3 uncR(s) 60.0 60.0 60.0 60.0 60.0 60.0 44.4 31.4 Erms(deg) 3.8453 3.1109 3.7296 3.4173 3.2838 4.7022 5.0076 7.5408 PA(deg) 56.8 54.168 57.983 54.720 53.784 59.386 18.208 20.714 PT (s) 1.587 1.690 1.682 1.745 1.997 2.455 58.602 58.501 uncT (s) 60.0 60.0 60.0 60.0 60.0 57.1 52.5 52.1 Erms(deg) 2.4531 2.7576 3.9682 3.6422 3.4596 3.6173 1.6189 2.9452 Source: Developed by the author 3.3 RESULTS AND DISCUSSION 35 Table 5 - IGA moderate response metrics. Metric H1 H2 H3 H4 H5 P1 P2 P3 uncR(s) 60.0 54.2 54.4 58.2 56.6 50.1 40.3 26.9 Erms(deg) 4.1903 3.8210 4.1228 3.7129 3.6487 6.0676 8.3207 10.0580 PA(deg) 29.178 25.202 30.432 32.911 31.595 16.555 6.3058 12.741 PT (s) 58.71 58.67 58.66 58.73 58.68 58.54 58.5 58.5 uncT (s) 60.0 60.0 60.0 60.0 60.0 53.2 52.2 40.0 Erms(deg) 2.7362 3.0609 3.9622 3.6388 3.5308 4.4485 3.6215 4.3497 Source: Developed by the author Table 6 - IGA critical response metrics. Metric H1 H2 H3 H4 H5 P1 P2 P3 uncR(s) 54.3 54.0 54.5 54.8 34.2 15.45 33.6 15.2 Erms(deg) 4.8935 4.5646 4.7113 4.2663 4.8730 7.7100 10.3250 11.8600 PA(deg) 24.863 22.284 26.444 29.970 27.652 16.365 3.652 12.049 PT (s) 58.69 58.67 58.68 58.75 58.71 56.7 58.5 58.52 uncT (s) 60.0 60.0 60.0 60.0 60.0 52.3 40.3 27.25 Erms(deg) 3.3695 3.5921 4.3305 3.9430 3.9376 5.6544 4.9096 5.8021 Source: Developed by the author In all cases, transient response presented interesting results, where it seems that stronger muscles present bigger overshoot. However, strong muscles demonstrate less sensitivity to external disturbances modeled in this research. Also, responses are according to reality where a healthy patient even in spite of non-idealities could track and regulate angular position very well and an SCI patient with strong muscles also compensates (P1), but not as well as a healthy one, and SCI patients with weak muscles do not reach well tracking and regulating results with non-idealities in the model. As commented before, it can be due to numerous factors that some papers did not take into account the existing problems with the SCI population assuming that results would be the same for healthy patients. It is noteworthy that responses to the step and sine trajectories are very similar to real exper- iments on SCI patients and healthy ones related in the literature (TEODORO, 2018; RIENER; QUINTERN; SCHMIDT, 1996; SHARMA et al., 2009; OLIVEIRA et al., 2017; STEGATH et al., 2007, 2008; FERRARIN; PEDOTTI, 2000; DOWNEY et al., 2015a) and in Chapter 4 of this research, proving effectiveness of simulating knee-joint and non-idealities that are reg- ularly encountered on real world. Thus, our methodology is based on simulation procedures optimizing a control task by an intelligent technique before implementing a controller on real tests. Moreover, this method will provide a lot of information about human identified system behavior to NMES/FES, permitting to save time and resources. 3.3 RESULTS AND DISCUSSION 36 Figure 10 - Analysis of tracking responses using empiric and IGA tuning on ideal conditions. (a) Step trajectory for patient H1. (b) Sine trajectory for patient H2. (c) Sine trajectory for patient H3. (d) Step trajectory for patient H4. (e) Sine trajectory for patient H5. (f) Step trajectory for patient P1. (g) Sine trajectory for patient P2. (h) Step trajectory for patient P3. Source: Developed by the author 3.3 RESULTS AND DISCUSSION 37 Figure 11 - Analysis of system responses using IGA tuning considering smooth non-idealities. (a) Tracking of both trajectories for patient H1. (b) Tracking of both trajectories for patient H2. (c) Tracking of both trajectories for patient H3. (d) Tracking of both trajectories for patient H4. (e) Tracking of both trajectories for patient H5. (f) Tracking of both trajectories for patient P1. (g) Tracking of both trajectories for patient P2. (h) Tracking of both trajectories for patient P3. Source: Developed by the author 3.3 RESULTS AND DISCUSSION 38 Figure 12 - Analysis of system responses using IGA tuning considering moderate non-idealities. (a) Tracking of both trajectories for patient H1. (b) Tracking of both trajectories for patient H2. (c) Tracking of both trajectories for patient H3. (d) Tracking of both trajectories for patient H4. (e) Tracking of both trajectories for patient H5. (f) Tracking of both trajectories for patient P1. (g) Tracking of both trajectories for patient P2. (h) Tracking of both trajectories for patient P3. Source: Developed by the author 3.3 RESULTS AND DISCUSSION 39 Figure 13 - Analysis of system responses using IGA tuning considering critical non-idealities. (a) Tracking of both trajectories for patient H1. (b) Tracking of both trajectories for patient H2. (c) Tracking of both trajectories for patient H3. (d) Tracking of both trajectories for patient H4. (e) Tracking of both trajectories for patient H5. (f) Tracking of both trajectories for patient P1. (g) Tracking of both trajectories for patient P2. (h) Tracking of both trajectories for patient P3. Source: Developed by the author 40 4 EXPERIMENTAL RESULTS WITH NEURAL NETWORK BLACK-BOX MODELS On the attempt to substantiate and validate the proposed methodology, experimental results are provided for seven healthy individuals and two paraplegic patients. An identification step based on neural networks black-box models is made using previous rehabilitation data as ses- sions passed by, aiming to have a better description on the relationship of delivered PW and achieved angular position, where fatigue and other problems as tremors are already implicit in the data. The following sections of this chapter describe the instrumentation, analyzed individuals, and the procedures used during the real experiments; the results and its analysis; and finally, a conclusion for this chapter is presented. 4.1 INSTRUMENTATION Figure 14 illustrates the test platform used for conducting the experiments at the Instrumen- tation and Biomedical Engineering Laboratory (“Laboratório de Instrumentação e Engenharia Biomédica - LIEB”) at UNESP - Ilha Solteira. This platform is composed of a NI myRIO con- troller (National Instruments R⃝) to operate in real time; a current-based neuromuscular electrical stimulator developed by Sanches (2013); an instrumented chair composed by an electrogo- niometer (Lynx R⃝), a gyroscope LPR510AL (ST Microelectronics R⃝), two triaxial accelerom- eters MMA7341 (Freescale R⃝); and two user interfaces developed in the LabVIEW R⃝ student version edition, one for identification and the other for controlling. The neuromuscular electrical stimulator delivers rectangular, biphasic, symmetrical pulses to the individual’s muscle, allowing a control adjustment of the pulse width in a range of 0− 400µs. The stimulation intensity is controlled by setting the pulse amplitude to the quadriceps and controlling the PW. The stimulus frequency was fixed in 50 Hz and the pulse amplitude in 80 mA (healthy subjects) or 120 mA (paraplegic subjects). Further, surface electrodes with rectangular self-adhesive CARCI 50 mm x 90 mm settings are used in this study. 4.2 ANALYZED INDIVIDUALS The study with volunteers was authorized through a research ethics committee involving human beings (CAAE: 79219317.2.1001.5402) at São Paulo State University (UNESP), and before the participation, written informed consent was obtained from all participants. Table 7 4.2 ANALYZED INDIVIDUALS 41 Figure 14 - The complete test platform for electrical stimulation experiments. Source: Developed by the author presents sufficient data on all male individuals that participated in this study. Table 7 - Specific data on analyzed individuals. H1 H2 H3 H4 H5 H6 H7 P1 P2 Age (years) 24 28 27 22 22 28 25 32 43 Weight (kg) 74.1 70.4 75 94.3 73 68.8 78.3 70.0 96.0 Height (cm) 174 167 180 186 175 170 165 170 183 Injury level - - - - - - - L4, L5 C5, C6 Injury time - - - - - - - 9 years 17 years Source: Developed by the author The inclusion criteria for individuals were male or female people older or equal than 18 years old. Healthy subjects were expected to have no lesion or congenital disease that could compromise the lower limbs movements, and paraplegic subjects were considered the ones with SCI. As exclusion criteria, pregnant women and people with cardiovascular diseases were not considered to this study. Furthermore, healthy individuals who presented fear or discomfort 4.3 EXPERIMENTAL SETUP 42 during the tests by voluntarily moving or holding the movement of the leg (excessively) were excluded from the analysis. By excessively, it means that as healthy individuals, small voluntary movements is probably present during stimulation and in this level, it is acceptable. 4.3 EXPERIMENTAL SETUP Initially, the chair backrest and the knee joint position are adjusted for each volunteer to ensure their comfort. A muscle analysis is made to determine the motor point following the methodology described in Teodoro (2018), in order to guarantee proper positioning of the sur- face electrodes. Therefore, scraping and cleansing procedures are made at the motor point to correctly place the electrodes. After the motor-point identification, a few open loop tests are per- formed applying a step-type signal during four seconds, in the interest of determining a bounded pulse width band ρmin to ρmax concerning to a desired range of angular position θmin = 10◦ and θmax = 40◦. Figure 15 illustrates this open-loop stimulation for subject H1. Subsequently, a small time interval for muscle rest is provided. Figure 15 - Open-loop tests to determine ρmin and ρmax. Source: Developed by the author During the experiments, individuals were instructed to relax, to not influence the leg mo- tion voluntarily and allow the stimulation to control it. Further, during electrical stimulation sessions, the individual could deactivate the stimulation pulses using a stop button, under any displeasure situations. In the following two subsections, the experimental setup is detailed, where, first, it is as- sumed a subject for the first time using the proposed methodology, i.e., without any previous data, and then, subjects who participated in two up to five sessions of experimental tests. 4.3 EXPERIMENTAL SETUP 43 4.3.1 First session Thereafter, a stimulation test is carried out consisting of one minute (less for paraplegic subjects due to premature fatigue problem) of randomly selected PW in the predetermined range (ρmin, ρmax) per individual, being each value applied during minimum four and maximum seven seconds (also random), as Figure 16 illustrates the stimulation test from individual H1. The motivation to follow this methodology is the attempt to map a tracking situation and to recognize the completely nonlinear and time-varying nature of muscles during a high time of electrical stimulation. The pulse width, position, velocity, and angular acceleration data are automatically recorded with 0.02 of time sampling, resulting in a dataset with approximately 3000 samples (60 seconds) at most. Figure 16 - Randomly selected PW during a random time. Source: Developed by the author Afterward, the identification data is read and manipulated for feeding up an MLP feedfor- ward NN with one hidden layer, aiming to map the relationship between the electric stimulus and the achieved angular position. In the literature, one layer has been proved to be sufficient to approximate any continuous function on a compact domain (PREVIDI, 2002; OGUNMOLU et al., 2016). The number of neurons varied from 10 to 1000; hyperbolic tangent activation is used in each neuron from the hidden layer and the output layer is composed of one neuron with linear activation, which gives the estimated output ŷ(k). The time-step value m and the past-step value n were chosen as 1 (for more details, refers to Section 2.3), where Table 8 illustrates an example on 10 samples of how the dataset is ar- ranged, containing the last input value “Pulse_Width(k− 1)” and the last output value “Angu- lar_Position (k− 1)” as features, and the actual output value “Angular_Position (k)” as target. 4.3 EXPERIMENTAL SETUP 44 The MLP NN model requires a normal input arranged as [samples, f eatures], where the obser- vations at previous timesteps are inputted as features to the model, considering that a feedfor- ward NN, as the MLP, can only take a fixed-size input. Table 8 - Example of how datasets are encoded. Features Target Angular_Position (k−1) Pulse_Width (k−1) Angular_Position (k) 13.348546◦ 215µs 13.399383◦ 13.399383◦ 215µs 13.382377◦ 13.382377◦ 215µs 13.325167◦ 13.325167◦ 215µs 13.306247◦ 13.306247◦ 215µs 13.347835◦ 13.347835◦ 248µs 13.387653◦ 13.387653◦ 248µs 13.357460◦ 13.357460◦ 248µs 13.256510◦ 13.256510◦ 248µs 13.131691◦ 13.131691◦ 248µs 13.016152◦ Source: Developed by the author Therefore, using the estimated model, an optimization procedure is performed based on the proposed IGA, introduced in Section 2.2, to find the best gains combination for two reference trajectories. The first trajectory is a sinusoidal wave ranging from 10◦ to 40◦ as an isotonic contraction with a repetitive pattern, and the second one is a 40◦ step wave as an isometric contraction. The motivation to use a smooth range of motion at 40◦ and a small-time period (sine wave) is the attempt to avoid premature fatigue by diminishing muscles effort. Considering real-world application of the proposed methodology and by assuming a limited time for a rehabilitation session, the initial parameters of IGA used for simulations were: size of population equal 8, size of RIP equal to 6, crossover rate equal 1, mutation rate equal 0.5, number of generations equal 6 and k equal 1 iteration. The algorithm ran a single time, and when the set of solutions was provided and gains combination selected for performing the real experiment, previous simulations to both trajectories were made to check the system’s response. Generally, the running time did not exceed 10 minutes of execution. The simulation system was developed using the Matlab/Simulink R⃝ platform, which con- tains both sine and step trajectories, a saturation block bounding control signal from 0 µs to ρmax µs for each subject, the RISE controller block, and the identified neural network block for each individual. Lastly, using empiric gains and then the gains encountered in the previous step, the con- trolling procedure is implemented for both trajectories. Data are recorded with 0.005 of time sampling, generally resulting in a dataset with approximately 12,000 samples (60 seconds) at most. 4.4 RESULTS AND DISCUSSION 45 4.3.2 Two up to five sessions Differently, for individuals who participated in more sessions with at least 48 hours of difference between two consecutive sessions, the identification step is not considered as it is done when an individual participates for the first time. The data from previous rehabilitation sessions are used for training a neural network model in an offline scheme, where before each next section, all data from a individual is combined to a single dataset and used to better map its relationship with electrical stimulus. Thus, using each trained identified model, an optimization procedure using the developed IGA is made to find the best gains combination for both sine and step reference trajectories. As this optimization is made in an offline scheme and before the next session, time and computa- tional costs are not too strict as they are for the first session. Therefore, the initial parameters of IGA used for simulations were: size of population equal 10, size of RIP equal to 30, crossover rate equal 1, mutation rate equal 0.3, number of generations equal 30 and k equal 1 iteration. The algorithm ran a single time, and several gains combinations from the set of solutions were simulated to check the system’s response and select the best one for both trajectories. Generally, this procedure took at most 1 hour. For the experimental part, the electrodes are positioned at the motor-point identified in the first session, and similarly, a few open-loop tests applying a step-type signal during four seconds are performed, to determine a bounded pulse width band concerning to θmin = 10◦ and θmax = 40◦. Afterward, a small-time interval for muscle rest is provided. Therefore, knowing the fine-tuned gains parameters for each individual, the controlling pro- cedure is applied for both references, and then with an empiric gains combination for comparing results. 4.4 RESULTS AND DISCUSSION Results and its discussions will be presented in this section divided by nine subsections, one for each individual that participated in the experiments during one up to five sessions at most. For all individuals will be presented: • a table containing additional technical information about each stimulation session. This table contains data related to initial identification, i.e., the range (ρmin, ρmax) in µs to acquire (θmin, θmax) in degrees. Further, it comprises data for the controlling procedure (empiric and IGA), where ρ0 in µs is the initial pulse width for starting the control operation; sat in µs is the value to bound control signal preventing possible pain to the individual (as tests are also made with healthy subjects) and quick fatigue due to the 4.4 RESULTS AND DISCUSSION 46 electric stimulus; and empiric and fine-tuned gains used to control-stimulate the lower limb. The two ρ0 and sat values were initially chosen according to the initial identification, i.e., concerning (ρmin, ρmax). As the tests advanced and fatigue was detected, both values were increased but under-estimated to prevent the problems aforementioned (pain, fatigue). Further, the same empiric gains were applied for both sine and step trajectories, and in this table, all unit of measurements is omitted due to page length limitation; • a table with information on the identified NN model such as the training time (TT), the correlation between input and output data (Corr), and the mean squared error (MSE); • a table presenting control metrics results for the experiments made. The metrics for evaluating the control system performance on experiments are the RMSE between the desired and actual knee angles considering the whole period of control-stimulation; and the time of effective control (TEC) representing how much time in seconds the lower limb is control-stimulated to track the reference angle. When the lower limb does not track the reference angle, the RMSE metric is represented by NC, meaning “not calculated”; • and to conclude, experimental results are presented in figures comparing both empiric and IGA tuning for all individuals. 4.4.1 Individual P1 Figure 17 illustrates subject P1 who participated in one session. As related by P1, during experiments no discomfort or pain was provoked by stimulation, there is no sensibility under the level of injury (L4, L5), no use of electrical stimulation was made before this session, and currently, none rehabilitation or treatment is being carried out. Hence, the initial procedure was performed for finding the subject’s P1 motor point and the bounded pulse width band (ρmin, ρmax) (just as Figure 15 illustrates), as described in Section 4.3. Thereafter, an stimulation of random PW during 40 seconds was made (similar to Figure 16) and the acquired data were treated and fed into a NN model to map the relationship of PW and angular position. Using the identified model, the IGA optimization was performed to find the best gains combination for both sine and step trajectories. During this time (about 10 minutes), muscle rested from previous stimulations. Therefore, the controlling procedure using gains encountered by the IGA was made for the sine and step waves, respectively during 30 seconds of stimulation at most. The operation point for the step wave was decreased to 30◦, aiming to avoid premature fatigue. To conclude the session, the same controlling procedure was made using empiric gains. 4.4 RESULTS AND DISCUSSION 47 Figure 17 - Individual P1 in the instrumented chair. Source: Developed by the author Therefore, Table 9 describes additional technical information about the stimulation session for subject P1; Table 10 presents information for the identified model, and metrics results on simulations; Table 11 presents control metrics results on experimental results and Figure 18 exhibits control results for subject P1 from the LabVIEW R⃝ platform (using IGA tuning only) and one with zoom to the tracking results (empiric and IGA tunings). As results in Table 11 and Figure 18 show, the proposed methodology could be effectively applied to clinical procedures for treating spinal cord injured patients via NMES/FES. As one will see in the following subsections, the RMSE for the sine wave from P1, was the best result during all experiments made during this research. However, it is notable in the final seconds (about 28 seconds) that the lower limb would start to have more tremors due to the fatigue factor. Similar results were acquired in Chapter 3, where paraplegic subjects presented superior tracking results under ideal conditions, i.e., without non-idealities included such as fatigue, tremors, and spasms. Using an empirical tuning, the lower limb tracked the sine wave with a lag and presented a slow response to the step wave, but still, tracking results are acceptable with the RISE controller. On the other hand, as previously mentioned, very good tracking results were achieved for both trajectories when the RISE controller was tuned via the proposed methodology. Table 9 - Technical information on experiments for individual P1. Identification Empiric Sine IGA Step IGA ρmin;ρmax θmin;θmax ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β 200;250 18;50 200;350 1;2;30;5 180;300 2.61;3.34;48.94;1.78 180;300 2.72;3.57;47.12;1.54 Source: Developed by the author 4.4 RESULTS AND DISCUSSION 48 Table 10 - Identification results for individual P1. Session TT Corr MSE 1st 28(s) 0.836 0.0019 Source: Developed by the author Table 11 - Metrics on experimental results for individual P1. Sine Step Empiric IGA Empiric IGA Session RMSE TEC RMSE TEC RMSE TEC RMSE TEC 1st 9.1471◦ 30(s) 2.9842◦ 30(s) 10.9950◦ 30(s) 5.9786◦ 25(s) Source: Developed by the author Results from P1 validate and substantiate the first hypothesis made in Section 1.3, whereby using an empirical tuning, the rehabilitation of SCI patients presents inadequate performances while by using the proposed approach control results are improved and fatigue or other prob- lems can be diminished. As SCI population have weak muscles or even atrophied ones due to a total or partial paralysis, the application of empiric gains would not provide good tracking, which would accentuate the fatigue factor, motivating the proposal of this research in the way that when the RISE controller is well-tuned, control performance is improved. Moreover, it is noteworthy that even using a higher pulse amplitude of 120 mA (in contrast with 80 mA for healthy subjects) both ρ0 and sat values are higher for paraplegic subjects due to aforemen- tioned problems (weak muscles or even atrophied ones), which allows the control law to switch in a higher range of pulse width to obtain an adequate muscle contraction and control it. From Table 10 it is notable a high correlation between input and output data. When the data present a high correlation between features and target, the generalization and learning are less difficult, which consequently generated a low MSE for the model. Figure 19 illustrates a comparison between simulation and real experiment for the sine wave. As the correlation between data was high, the model could identify quite well the result in comparison with the real experiment. Both results presented similar tremors in the upper and lower peak value, while the step wave was modeled with a fast response with no stationary error. However, as an approximate model using only the identification data of 40 seconds, this model was not prepared for modeling fatigue and tremors that were present in the real experiment. As we hypothesized in Section 1.3, the modeling of knee joint position under NMES/FES could be improved by taking into advantage the use of past rehabilitation data. However, as P1 did not continue participating in more sessions, we could not validate this second hypothesis. 4.4 RESULTS AND DISCUSSION 49 Figure 18 - Experimental results for individual P1. Source: Developed by the author 4.4.2 Individual P2 Just as subject P1, subject P2 illustrated in Figure 20 participated in one session. As related by P2, during experiments no discomfort or pain was provoked by stimulation, there exist pains 4.4 RESULTS AND DISCUSSION 50 Figure 19 - IGA comparison of simulation and real experiment for individual P1. Source: Developed by the author due to spasms contractions during the day, he feels small sensibility in the lower limbs, he has participated in experiments using Russian current, and currently, none rehabilitation or treatment is being carried out. Similarly to P1, the same procedure was performed for the stimulation session including identification and control. For more information, the reader is invited to check the previous Subsection 4.4.1. Therefore, Table 12 describes additional technical information about the stimulation session for subject P2; Table 13 presents information for the identified model, and metrics results on simulations; Table 14 presents control metrics results on experimental results and Figure 21 exhibits control results for subject P2 from the LabVIEW R⃝ platform (using IGA tuning only) and one with zoom to the tracking results (empiric and IGA tunings). As one can notice in Figure 21 and Table 14, tracking results for subject P2 were not as good as they were for P1. However, the poor sine wave tracking under IGA tuning could be due to an under-estimation for saturation sat that bounds the control signal value, where it seems that a possible good tracking would be acquired if the saturation value had been selected higher. This idea is substantiated by results of the step wave of 30◦ performed with three minutes interval for muscle rest after the sine wave using the IGA tuning where a very good regulation around the operation point was achieved during 21 seconds approximately. Table 12 - Technical information on experiments for individual P2. Identification Empiric Sine IGA Step IGA ρmin;ρmax θmin;θmax ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β 150;250 8;38 200;370 1;2;30;5 185;310 2.22;3.54;39.50;1.40 190;360 3.01;1.91;48.34;2.65 Source: Developed by the author Additionally, similar to results from P1, an empirical tuning provided very poor perfor- mance for both sine and step trajectories, by not tracking the sine wave and by regulating with 4.4 RESULTS AND DISCUSSION 51 Figure 20 - Individual P2 in the instrumented chair. Source: Developed by the author Table 13 - Identification results for individual P2. Session TT Corr MSE 1st 33(s) 0.796 0.0032 Source: Developed by the author Table 14 - Metrics on experimental results for individual P2. Sine Step Empiric IGA Empiric IGA Session RMSE TEC RMSE TEC RMSE TEC RMSE TEC 1st 11.2966◦ 30(s) 10.7306◦ 30(s) 10.1067◦ 23(s) 6.6134◦ 21(s) Source: Developed by the author a stationary error around the operation point from the step wave. In contrast, when using the proposed methodology, even the one with poor performance due to under-estimation value for the saturation value, results are better than using empiric gains. For the step wave, the RMSE metric decreased almost 50%, representing very good tracking results. From Table 13 it is notable a high correlation between input and output data, which gener- ated a low MSE for the model. As the correlation between data was high, the model identified quite well the real experiment responses. However, due to the under-estimated value of sat, fatigue damping and tremors from real-world muscle conditions, the results in real experiments varied from simulations. 4.4 RESULTS AND DISCUSSION 52 Figure 21 - Experimental results for individual P2. Source: Developed by the author 4.4.3 Individual H1 Individual H1 participated in five sessions, wherein its first session, the initial procedure was carried out for finding its motor point and the bounded pulse width band, as described in 4.4 RESULTS AND DISCUSSION 53 Section 4.3. Thereafter, an stimulation test was implemented (as Figure 16 illustrates) and the acquired data were fed into a NN model to map the relationship of PW and angular position. While the NN was being trained and the IGA optimization later on performed, the controlling procedure with empiric gains was made for both sine and step-signal trajectories respectively. To conclude the session, the same controlling procedure was made with the gains found by the IGA for the identified NN model. As described in Section 4.3.2, to the second until the fifth session, the initial procedure was carried out and there was no identification step given that the NN model was offline trained with both control and identification previous data. In the second, third and fifth session, at first, the controlling procedure was performed with the gains found by the IGA optimizer, and then with empiric ones, while in the first and fourth sessions, empiric gains were applied first and then the IGA ones. The first session took more time and an additional stimulation than the following others due to the motor-point identification, stimulation test of one minute, and the training/optimization time made to find the best gains combination. Therefore, Table 15 describes additional technical information about each session, being each line of the table a session day for individual H1; Table 16 presents information for the identified models, and metrics results on simulations; Table 17 presents control metrics results on experimental results and Figures 22 and 23 display control results for subject H1 in all stimulation sessions. Table 15 - Technical information on experiments for individual H1. Identification Empiric Sine IGA Step IGA ρmin;ρmax θmin;θmax ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β 100;150 10;42 100;170 1;2;30;5 100;170 3.23;1.08;24.74;5.50 100;170 1.37;1.63;54.03;2.36 80;175 10;40 100;190 0.5;1;30;1.5 100;190 1.76;2.28;32.30;2.39 100;190 0.64;1.66;52.26;4.00 100;150 12;45 100;185 0.8;1.2;20.5;2.5 100;150 3.23;2.52;27.33;2.29 100;165 2.30;4.24;59.26;3.49 120;160 10;40 100;180 5;2;15;3 100;180 2.40;4.10;27.05;2.18 100;180 3.12;5.80;43.162;1.35 140;190 10;44 120;220 4;7;25;8 120;210 3.07;4.37;21.73;1.56 120;210 2.61;3.54;39.50;1.30 Source: Developed by the author Table 16 - Identification results for individual H1. Session TT Corr MSE 1st 34(s) 0.726 0.0022 2nd 48(s) 0.157 0.0386 3rd 106(s) 0.101 0.0428 4th 144(s) 0.159 0.0411 5th 350(s) 0.113 0.0402 Source: Developed by the author As demonstrated in Figures 22 and 23, an empirical tuning could also guarantee stability, i.e., the lower limb tracked the step wave effectively in sessions one and four, and in most of the sine trajectories, the leg poorly, but tried to follow the reference angle. However, as previously 4.4 RESULTS AND DISCUSSION 54 Table 17 - Metrics on experimental results for individual H1. Sine Step Empiric IGA Empiric IGA Session RMSE TEC RMSE TEC RMSE TEC RMSE TEC 1st 7.4946◦ 60(s) 5.8303◦ 60(s) 5.9209◦ 60(s) 6.1677◦ 60(s) 2nd 8.7529◦ 60(s) 5.9333◦ 60(s) 12.2910◦ 60(s) 8.2012◦ 60(s) 3rd 14.0928◦ 60(s) 7.3373◦ 60(s) 7.2664◦ 60(s) 4.1644◦ 60(s) 4th 6.3778◦ 60(s) 3.6292◦ 60(s) 6.7410◦ 60(s) 4.4040◦ 60(s) 5th 6.3835◦ 60(s) 3.5624◦ 60(s) 6.8875◦ 60(s) 4.4251◦ 60(s) Source: Developed by the author Figure 22 - Experimental results for subject H1 on sessions one and two respectively. Source: Developed by the author 4.4 RESULTS AND DISCUSSION 55 discussed in Section 3.3, once gains selections are immense, one likely choose combinations that would not guarantee best performance and accuracy during rehabilitation, e.g., compare simulation results from Figure 10 with the result from session five in Figure 23. As one can notice in Table 17 and in the two previous figures, results for individual H1 using an empirical tuning provided poor tracking performances in most of the sessions, while with an appropriate choice of gains, using the proposed methodology, the lower limb presented better tracking results. Except for session two, results to the step wave presented very efficient regu- lation around the desired point of operation, and for the sine wave, sessions two and three did not present very positive results, but still, the control law tried to compensate the error between desired and real angle. Furthermore, the RMSE tended to decrease using fine-tuned gains as sessions passed by, showing that the proposed methodology could provide good rehabilitation treatments for SCI patients that continuously go on sessions of NMES/FES. Additionally, in all stimulation sessions, even the ones with poor performances, the control- stimulated lower limb tried robustly to track the reference angle during 60 seconds, while the RISE controller in the literature presented by Stegath et al. (2007, 2008) demonstrated track- ing control on 8 seconds for the stepping trajectory and 20 seconds for a sine wave; Sharma et al. (2009, 2012) presented tracking control on 30 seconds for a step and sine-type signal; and Kushima et al. (2015), Downey et al. (2015b) presented tracking control on 30 and 45 seconds respectively for a sine trajectory. From Table 16, it is noticed a fast decrement from the first session to the remaining ones in the correlation between input and output data as sessions passed by (given the addition of control data from every new session), which resulted in increments to the MSE metrics for each trained NN model in comparison with the first one. Due to less correlation between data, the generalization and learning procedure of a NN is much harder, however, this model better describes what happens in real experiments, where nonlinearities and the time-varying behavior of muscles are explained by data. 4.4 RESULTS AND DISCUSSION 56 Figure 23 - Experimental results for subject H1 on sessions three, four and five respectively. Source: Developed by the author 4.4 RESULTS AND DISCUSSION 57 In this context, Figure 24 illustrates a comparison between results from simulation and real experiments from the individual H1 for sessions one and four using the IGA tuning approach. In the first session, the data were much more correlated and the model could predict quite well what would happen in the real experiment, wherein responses for the step wave are very comparable. Similarly, for the fourth session, the model did also provide a good description of the real experiment, predicting a sine trajectory with tremors in both lower and upper values of the wave and a quick response to the step-signal trajectory with no stationary error. However, as an approximate model, neither of them would be prepared for an exact description of the real system, for instance, considering that voluntary movements, fear or other issues related to thoughts of the individual itself (e.g., social or personal life) that can affect results, could not be predicted by the model. Figure 24 - IGA comparison of simulation and real experiment for individual H1. Source: Developed by the author 4.4.4 Individual H2 Individual H2 also participated at five sessions, where the set-up protocol for control- stimulation was similar as it was for H1 (for more details, refers to the previous Subsection 4.4.3 or Section 4.3). Therefore, Table 18 describes additional technical information about each ses- sion, being each line of the table a session day for individual H2; Table 19 presents information for the identified models, and metrics results on simulations; Table 20 presents control metrics results on experimental results; and Figures 25 and 26 exhibits control results for subject H2 in all stimulation sessions. 4.4 RESULTS AND DISCUSSION 58 Table 18 - Technical information on experiments for individual H2. Identification Empiric Sine IGA Step IGA ρmin;ρmax θmin;θmax ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β ρ0;sat α1;α2;ks;β 120;150 10;40 100;200 1;2;30;5 100;220 1.90;3.50;48.00;3.00 100;220 1.90;3.50;48.00;3.00 180;220 10;40 200;290 0.8;1.2;20.5;2.5 180;240 1.57;2.37;48.45;1.05 200;260 1.38;1.34;64.41;3.72 150;180 14;41 130;200 0,5;1;30;1,5 150;230 1.47;3.31;30.01;1.87 170;250 3.54;3.83;54.88;1.92 130;150 12;42 150;220 5;2;15;3 110;170 1.42;3.68;35.09;1.99 130;190 2.07;1.75;36.32;1.96 130;160 8;40 120;230 4;7;25;8 120;180 2.03;3.02;36.07;2.23 120;200 2.17;1.76;38.53;1.85 Source: Developed by the author Table 19 - Identification results for individual H2. Session TT Corr MSE 1st 94(s) 0.869 0.0064 2nd 132(s) 0.416 0.0397 3rd 148(s) 0.308 0.0390 4th 249(s) 0.282 0.0370 5th 327(s) 0.292 0.0382 Source: Developed by the author Table 20 - Metrics on experimental results for individual H2. Sine Step Empiric IGA Empiric IGA Session RMSE TEC RMSE TEC RMSE TEC RMSE TEC 1st 5.2124◦ 60(s) 5.0555◦ 44.9(s) 9.7644◦ 57(s) 6.2122◦ 37(s) 2nd 8.3179◦ 60(s) 3.8853◦ 60(s) NC 30(s) 7.8566◦ 25(s) 3rd 11.7414◦ 35(s) 3.6331◦ 60(s) 11.8222◦ 57(s) 5.4572◦ 37(s) 4th 4.8877◦ 60(s) 3.5629◦ 40(s) 6.4240◦ 34(s) 4.8901◦ 45(s) 5th 10.7136◦ 33(s) 4.8582◦ 23(s) 6.2262◦ 35(s) 7.2336◦ 38(s) Source: Developed by the author As one can notice in Table 20 and in Figures 25 and 26, the lower limb barely tracked the reference trajectories when control-stimulated using empiric gains constants in most of sessions. In sessions one and five, for the sine trajectory, the lower limb oscillated the whole period of stimulation, being more accentuated in the fifth session, which is very similar to simulation results in Figure 10 and results obtained for individual H1 in Figure 23. Additionally, for the step-signal trajectory, in the first and fourth session, the lower limb tracked the reference angle very well during 30 seconds approximately, and in the last one, it took a while to pass through a very large (20 seconds) transitory period but still presented a good tracking result. On the other side, by using the proposed methodology to tune the RISE controller, results to the sine wave in all sessions were less or equal to 5◦ with minimum 23 seconds and at most 60 seconds (two of five) of stimulation, representing very satisfactory tracking results. Moreover, for the step-signal trajectory, the control-stimulated lower limb was regulated around 4.4 RESULTS AND DISCUSSION 59 the desired angle, with a poor result in session three where the oscillatory period took 22 se