Automation in Anesthesia: A Look at L1 Adaptive andPID ControllersbyKousha TalebianBachelor of Engineering Physics - Electrical Option, The University of BritishColumbia, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Biomedical Engineering)The University Of British Columbia(Vancouver)December 2016c© Kousha Talebian, 2016AbstractControl of anesthesia is one of the many tasks performed by anesthesiologists dur-ing surgery. It involves adjusting drug dosage by monitoring patient’s vital andclinical signs. A control system can replace this tedious and routine task, and al-low the anesthesiologists to concentrate on more life threatening procedures.Because of large intra- and inter-variability in patients Pharmacokinetics andPharmacodynamics responses, an adaptive controller is desirable. This thesis thor-oughly investigates the L1 Adaptive Control by applying it on 44 simulation caseswhich cover a wide range of patient demographics. It is found that the controllerapproaches an implantable non-adaptive LTI controller as the adaptation gain in-creases, echoing the results found by other researches. This loss of adaptivity isshown through examples and mathematical derivations. It is concluded that theL1 Adaptive Control in its current form is not applicable to closed-loop control ofanesthesia.As an alternative to adaptive controller, partial adaptivity in a PID controller isinvestigated. iControl, a PID controller designed by us, can sometimes lead to os-cillation in the control signal. It is desirable to automatically detect the oscillationsand tune the controller in order to remove them. A real-time oscillation detectionalgorithm is discussed. It detects multiple oscillations in real-time and providestheir frequency, amplitude, severity and regularity. A PID auto-tuning algorithmis developed that uses the dominant frequency metrics provided by the oscillationdetection algorithm to retune the controller robustly and to guarantee stability. Thistechnique is simulated and tested on 44 cases; the gain and the phase margin in all44 cases are within < 7% of the optimal tuning parameters of the iControl.iiPrefaceAll of the work of this research was conducted at BC Children’s Hospital in col-laboration with The University of British Columbia’s Electrical and Computer En-gineering in Medicine department. All clinical data used was approved by TheUniversity of British Columbia Children’s and Women’s Research Ethic Board(certificate H10-01174).I was the principal investigator on analysis of the L1 Adaptive Control (L1-AC)of Chatper 3 and its application to closed-loop control of anesthesia. A version ofChapter 3 is published in a journal article:• K. van Heusden, K. Talebian, G. A. Dumont. Analysis of L1 adaptive statefeedback control. Why does it approximate an implementable LTI con-troller? In European Journal of Control, 2015The theory of oscillation detection algorithm is based on Wang et a.l [67] whilethe implementation, parameter tuning, and MATLAb code is based on my researchand I was the lead investigator on its fisibility in clinical settings. The turningrules of Chapter 5 is from A˚stro¨m et al. [6], however the implementation and theMATLAB code was provided by myself.The Appendix Bwas a continuation of the work originally conducted by Soltesz,G. I applied his findings on a series of clinical data to determine and compare theirperformance to the well known Varvel metrics. These findings were presented atthe American Society of Anesthesiologist 2013:• K. Talebian, K. Soltesz, G. A. Dumont, M. Ansermino. Clinical assessmentof control performance in closed-loop anesthesia. In American Society ofAnesthesiologist, 2013.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Objectives and Scopes . . . . . . . . . . . . . . . . . . . . . . . 31.3 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . 42 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Monitoring Depth of Hypnosis . . . . . . . . . . . . . . . . . . . 82.1.1 Bispectral Index . . . . . . . . . . . . . . . . . . . . . . 82.1.2 Wavelet-based Anesthetic Value . . . . . . . . . . . . . . 92.2 Drug Effect Modeling . . . . . . . . . . . . . . . . . . . . . . . . 92.2.1 Pharmacokinetics of Propofol . . . . . . . . . . . . . . . 92.2.2 Pharmacodynamics of Propofol . . . . . . . . . . . . . . 10iv2.3 Automatic Control of Anesthesia . . . . . . . . . . . . . . . . . . 112.3.1 Closed-Loop Control: A Review . . . . . . . . . . . . . . 112.3.2 Closed-Loop Control: Oscillation . . . . . . . . . . . . . 132.3.3 Closed-Loop Control: Adaptive vs PID . . . . . . . . . . 142.4 Performance of Closed-Loop Anesthesia . . . . . . . . . . . . . . 143 L1 Adaptive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2 The L1 Adaptive Control . . . . . . . . . . . . . . . . . . . . . . 243.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . 243.2.2 State Predictor . . . . . . . . . . . . . . . . . . . . . . . 253.2.3 L1-norm Stability Condition . . . . . . . . . . . . . . . . 263.3 Achievable Performance Bound . . . . . . . . . . . . . . . . . . 263.3.1 Reference Controller . . . . . . . . . . . . . . . . . . . . 273.3.2 Control Performance . . . . . . . . . . . . . . . . . . . . 283.4 Case Studies of the L1 Controller . . . . . . . . . . . . . . . . . . 293.5 Loss of Adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . 333.5.1 Simple Example of Loss of Adaptivity . . . . . . . . . . . 353.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444 Real-Time Oscillation Detection . . . . . . . . . . . . . . . . . . . . 454.1 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.2 The Discrete Cosine Transform (DCT) Off-line Oscillation Detection 474.2.1 The DCT Definition . . . . . . . . . . . . . . . . . . . . . 474.2.2 The DCT Algorithm . . . . . . . . . . . . . . . . . . . . . 484.2.3 Summary of the DCT Algorithm . . . . . . . . . . . . . . 564.3 Extension to Real-Time . . . . . . . . . . . . . . . . . . . . . . . 574.3.1 Summary of the Real-Time DCT Algorithm . . . . . . . . 594.4 Oscillation Detection Examples . . . . . . . . . . . . . . . . . . 594.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615 Re-Tuning of a Proportional-Integral-Derivative (PID) Controller . 635.1 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 645.2 Overview of Controller Structure . . . . . . . . . . . . . . . . . . 64v5.2.1 iControl Design Structure . . . . . . . . . . . . . . . . . 665.3 Robustness and Performance Design . . . . . . . . . . . . . . . . 695.4 PID Auto-Tuning Rules . . . . . . . . . . . . . . . . . . . . . . . 715.4.1 Auto-Tuning for Robustness . . . . . . . . . . . . . . . . 725.4.2 Auto-Tuning for Performance . . . . . . . . . . . . . . . 745.4.3 Bumpless Parameter Change . . . . . . . . . . . . . . . . 745.5 Auto-Tuning Implementation . . . . . . . . . . . . . . . . . . . . 755.6 Summary of Auto-Tuning Algorithm . . . . . . . . . . . . . . . . 765.7 Simulation Examples and Results . . . . . . . . . . . . . . . . . 765.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816.1 Summary and Contributions . . . . . . . . . . . . . . . . . . . . 816.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 836.2.1 Imminent Future Direction . . . . . . . . . . . . . . . . . 836.2.2 Distant Future Direction . . . . . . . . . . . . . . . . . . 84Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85A Propofol PKPD Modeling . . . . . . . . . . . . . . . . . . . . . . . . 92A.1 Pharmacokinetics . . . . . . . . . . . . . . . . . . . . . . . . . . 92A.2 Pharmacodynamics . . . . . . . . . . . . . . . . . . . . . . . . . 96A.3 The PKPD Model . . . . . . . . . . . . . . . . . . . . . . . . . . 97B Control Performance in Closed-Loop Anesthesia . . . . . . . . . . . 99B.1 Varvel Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . 99B.2 Proposed Control Performance Measures . . . . . . . . . . . . . 101C Limiting Behavior of L1 Adaptive Control . . . . . . . . . . . . . . . 104C.1 Problem Formulation and The L1 Adaptive Controller . . . . . . . 104C.2 Removal of the Internal Feedback over ηˆ(t) . . . . . . . . . . . . 106C.3 Linearizing the L1 Controller with Generic Adaptation Laws . . . 108C.4 Linearization of the Projection Operator in the L1 Adaptive Control 111D Robustness and Performance of iControl . . . . . . . . . . . . . . . 114viE Oscillation Detection MATLAB . . . . . . . . . . . . . . . . . . . . . 119F PID Tuning Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 129viiList of TablesTable 4.1 High and Low components of Example 4.4.1 . . . . . . . . . . 60Table 4.2 Case example from Example 4.4.2. The magnitude regulatoryindex is 3.785. . . . . . . . . . . . . . . . . . . . . . . . . . . 61Table A.1 Propofol Pharmacokinetics (PK) parameters from [48]. BWstands for body weight, ven= 0 is for arterial sampling, ven= 1is for venous sampling, bol = 0 is for infusion administration,and bol = 1 is for bolus administration. . . . . . . . . . . . . . 95Table A.2 Parameter estimates of the the PK model of Table A.1 from [48]. 95Table A.3 PK and Pharmacodynamics (PD) parameters from the Bibianstudy [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Table B.1 Varvel and proposed measures of the example from Figure B.1. 103Table D.1 Robustness comparison of the iControl vs the auto-tunedalgorithm of Chapter 5 for the 44 Pharmacokinetics/Pharmaco-dynamics (PKPD) models. . . . . . . . . . . . . . . . . . . . . 115Table D.2 Output disturbance rejection comparison of the iControl tuningvs the auto-tuned algorithm of Chapter 5 for the 44 PKPD models.116Table D.3 Set-point response comparison of the iControl tuning vs theauto-tuned algorithm of Chapter 5 for the 44 PKPD models. . . 117Table D.4 PID Parameters of the iControl tuning and the auto-tunedalgorithm of Chapter 5 for the 44 PKPD models. . . . . . . . . 118viiiList of FiguresFigure 2.1 MDAPE vs IAE for a systematic small error. . . . . . . . . . 16Figure 2.2 MDAPE vs IAE for a sporadic error. . . . . . . . . . . . . . . 17Figure 2.3 Depth of Hypnosis (DOH) is clearly negatively biased. . . . . 18Figure 2.4 DOH is less biased. . . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.5 Small Induction Phase Duration (ID) of 3.1 min translates to asmall overshoot of 12.8%. . . . . . . . . . . . . . . . . . . . 20Figure 2.6 The large ID of 5.6 min translates to a larger overshoot of48.9% and a longer DOH settling time to set-point. . . . . . . 21Figure 3.1 The original L1-AC block diagram as it appears in [26]. . . . . 26Figure 3.2 The step response of patient #7 and the predictor model . . . . 29Figure 3.3 Simulated system output for the patient #7 controlled by theL1 controller. The upper plot shows the patient’s output forthe three adaptation gains. The reference output is shown inthe thick green line. The lower plot shows the absolute errorof patient’s output to the reference’s output. . . . . . . . . . . 31Figure 3.4 Simulated system output for the patient #2 controlled by theL1 controller. The upper plot shows the patient’s output forthe three adaptation gains. The reference output is shown inthe thick green line. The lower plot shows the absolute errorof patient’s output to the reference’s output. . . . . . . . . . . 32Figure 3.5 Simple feedback for a nonlinear function f (·) . . . . . . . . . 33Figure 3.6 Final form of the linearized adaptation laws for L1 adaptivecontrol with generic adaptation laws. . . . . . . . . . . . . . . 34ixFigure 3.7 Simulated output for the plant G(s) controlled by the L1-ACfor Law A). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is adaptive for low adaptation gain,however it becomes static for higher gains. . . . . . . . . . . 38Figure 3.8 Simulated output for the plant G(s) controlled by the L1-ACfor Law B). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is adaptive for low adaptation gain,however it becomes static for higher gains. . . . . . . . . . . 39Figure 3.9 Simulated output for the plant G(s) controlled by the L1-ACfor Law C). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is nonadaptive for all gains. . . . . . . 40Figure 3.10 Simulated output for the plant G(s) controlled by the L1-ACfor Law D). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is nonadaptive for all gains. . . . . . . 41Figure 3.11 Simulated output for the plant G(s) controlled by the L1-ACfor Law E). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is nonadaptive for all gains. . . . . . . 42Figure 3.12 Simulated output for the plant G(s) controlled by the L1-ACfor Law F). The top figure is the output of the plant, with thethick green line being the output of the reference. The lowerplot is the absolute difference of model output and referenceoutput. The controller is nonadaptive for all dead-zone intervals. 43xFigure 4.1 The signal blue has two components with frequencies 4 and 12units respectively (labeled red and green signals). The DCTof the blue signal is shown on the bottom graph. The first 4points on the DCT captures one of the frequency components,while the next 5 captures the other. . . . . . . . . . . . . . . . 49Figure 4.2 The reconstructed frequency components of Figure 4.1 isshown. The reconstruction can only preserve the frequency,but does not provide an accurate information on the magnitudeand offset of the signal. . . . . . . . . . . . . . . . . . . . . . 50Figure 4.3 x(t) and the frequency component xi(t) . . . . . . . . . . . . 51Figure 4.4 Signal from Example 4.4.1 contains two oscillations of periods1.3 min and 3.4 min. Dominant oscillation period is detectedat 3.397 min. The signal is shown in black and the dominantoscillation is shown in red. . . . . . . . . . . . . . . . . . . . 60Figure 4.5 On-line oscillation detection shows a detected dominant signalof Tp = 3.63 min, M¯ = 5.96 and F = 89.31%. . . . . . . . . . 61Figure 5.1 A 2-degree-of-freedom PID controller . . . . . . . . . . . . . 65Figure 5.2 The iControl Structure . . . . . . . . . . . . . . . . . . . . . 66Figure 5.3 The Nyquist Plot with region of stability . . . . . . . . . . . . 70Figure 5.4 The model is identified at the unstable point Pu and thecontroller is tuned to take the loop function to the stable point Ps. 72Figure 5.5 Group 1, Case 10: The original tuning has Am = 8.72 andφm = 60.29. The re-tuned system has Am = 8.62 and φm = 53.81. 77Figure 5.6 Group 2, Case 17: The original tuning has Am = 6.27 andφm = 61.20. The re-tuned system has Am = 6.13 and φm = 55.64. 77Figure 5.7 Group 3, Case 33: The original tuning has Am = 7.57 andφm = 65.21. The re-tuned system has Am = 7.43 and φm = 58.96. 78Figure 5.8 Group 4, Case 38: The original tuning has Am = 6.89 andφm = 61.14. The re-tuned system has Am = 6.86 and φm = 55.71. 78xiFigure A.1 The 3-compartment pharmacokinetic model. The rapidlyequilibrating compartment models the muscles and viscera.The central compartment models the blood, brain and liver.The slowly equilibrating compartment models the bones and fat. 93Figure A.2 The full PKPD model introduced in [10]. . . . . . . . . . . . . 97Figure B.1 Representative example from a closed-loop DOH controlsystem. The induction phase is shown in solid black, themaintenance phase is shown in solid blue, the emergencephase is shown in solid magenta, the reference is shown inthick green, and the r±10 bounds are shown in dashed black.The red dot represents the overshoot. . . . . . . . . . . . . . 103Figure C.1 L1-AC as formulated in [26]. . . . . . . . . . . . . . . . . . . 106Figure C.2 Equivalent architecture of the L1-AC with removed internalfeedback over ηˆ(t). . . . . . . . . . . . . . . . . . . . . . . . 107Figure C.3 Simplified architecture of the L1-AC to a more coherentstructure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure C.4 Simplified architecture of the L1-AC with generic adaptive laws. 109Figure C.5 Linearized L1-AC with generic adaptation laws. . . . . . . . . 109Figure C.6 Equivalent form of Figure C.5 of the linearized L1-AC withgeneric adaptation laws. . . . . . . . . . . . . . . . . . . . . 110Figure C.7 Final form of the linearized adaptation laws for L1-AC withgeneric adaptation laws. . . . . . . . . . . . . . . . . . . . . 111xiiGlossaryBIS Bispectral IndexCNS Central Nervous SystemDCT Discrete Cosine TransformDFT Discrete Fourier TransformDOA Depth of AnesthesiaDOH Depth of HypnosisDWT Discrete Wavelet TransformEEG ElectroencephalographyER Emergence Phase Rise TimeGABA Gamma-Aminobutyric AcidGS Global ScoreIAE Integrated Absolute ErrorIDCT Inverse Discrete Cosine TransformID Induction Phase DurationIE Integrated ErrorL1-AC L1 Adaptive ControlxiiiLBM Lean Body MassLTI Linear Time InvariantMAP Mean Arterial PressureMDAPE Median Performance Absolute ErrorMDPE Median Performance ErrorMIMO Multi-Input/Multi-OutputMPC Model-Predictive-ControllerOS Percent OvershootPD PharmacodynamicsPEEG processed ElectroencephalographyPE Percent ErrorPID Proportional-Integral-DerivativePI Proportional-IntegralPKPD Pharmacokinetics/PharmacodynamicsPK PharmacokineticsSISO Single-Input/Single-OutputSNR Signal to Noise RatioTCI Target-Controlled InfusionVI Variability IndexWAV Wavelet-based Anesthetic ValuexivAcknowledgmentsThis research would not have been possible without the support of the followingamazing individuals.First, I would like to express my deepest gratitude to my research supervisor,Dr. Guy Dumont, for his guidance, valuable support and insights. I would alsolike to extend my gratitude to Dr. Mark Ansermino who always provided the timeand support to answer my questions. Lastly, I would like to thank my colleagueDr. Klaske van Heusden who always discussed problems and helped me with thechallenges of the research.Second, and most importantly, I would like to thank my parents, Nora andMano, for their constant support of my professional career, and my life as a whole.I owe everything I have achieved and hope to achieve to their love and support.In addition, I would like to thank my dear sister, Kimia, who always showed hersupport without being asked.Third and finally, I would like to thank all my colleagues at the Pediatric Anes-thesia Research Team, Sara Khosravi, Nick West, Parastoo Dehkordi, Joanne Lim,Dr. Matthias Go¨rges, and Aryannah Umedaly.xvChapter 1Introduction1.1 MotivationThe number of surgeries are increasing in the US and around the globe. Accordingto a 2010 report by the Centers for Disease Control and Presentation, 51.4 millionpatients in the US alone went under surgery that year [17]. These surgeries rangeanywhere from a life saving operation, to a plastic surgery, or to correcting andrestoring the physical appearance [17].Anesthesiologists play a central role in the health and comfort of the patientpre-, intra-, and post-operatively. During the procedure, they monitor and controlpatient’s vital functions, including breathing, blood pressure, heart rate, body tem-perature, and etc. They ensure the safety and comfort of the patient by controllingthe hypnotic (unconsciousness) and the analgesic (sense of pain) states as well asassisting the surgeons by controlling the paralysis (relaxation and immobility of theskeletal muscle) through administration of anesthetic drugs. General anesthesia isthe term given to this complex state induced patient.To achieve general anesthesia, the anesthesiologist administers a variety ofdrugs. The combination of hypnotic, opioid, and neuromuscular drugs achievesthe three functional states of anesthesia: hypnosis, analgesic, and paralysis.Anesthesiologists assess adequate anesthesia and analgesia through monitoringthe patient’s vital and other clinical signs such as heart rate, blood pressure, eyemovement, pupil diameter, respiratory rate, facial grimacing, and lacrimation [18].1Depth of Anesthesia (DOA) is a measure of the effect of the hypnotic and analgesicdrugs that cause the unconsciousness and alleviate pain [52].The paralysis state does not contribute to the DOA [66]. The muscle relaxantdrugs that induce paralysis can however, effect some of the patient’s clinical signs,for example the respiratory rate. Anesthesiologists may therefore need to monitorother physiological signs as the paralysis state is induced.There is currently no single variable that can accurately measure the DOA andthe search for this metric remains an active field of research. Monitors use theElectroencephalography (EEG) - the electrical activity of the brain along the scalp,to assess Depth of Hypnosis (DOH)1. Interpretation of the raw EEG signal is timeconsuming and requires a trained neurologist. Therefore, different processed Elec-troencephalography (PEEG) methods have been applied to create an index thatquantifies the state of hypnosis. The most commonly used index, the BispectralIndex (BIS), uses the bispectral analysis of the EEG wavelengths [50]. Our researchteam at the University of British Columbia (UBC) has worked extensively in thisfield; their efforts have resulted in the introduction of the Wavelet-based AnestheticValue (WAV) index that compares well with the standard BIS, see [12], [69], and[11].Opioid (narcotic) analgesic drugs produce their effect through the interactionwith the Gamma-Aminobutyric Acid (GABA) receptors in the Central Nervous Sys-tem (CNS), though they bind at different sites [37] than the hypnotic drugs. Whilethere are proposed indices such as the Analgesia Nociception Index [36], nocicep-tion and antinociception measurements during anesthesia have not been clinicallyproven, see [24] and [30].Anesthesiologists continuously change the administration rate of the hypnotic,analgesic and relaxant drugs to account for stimuli from surgical incision. In fact,they assume the role of a feedback controller; to achieve a given clinical target, thedoctors monitor the clinical signs of the patient and adjust the drugs accordingly. Inmany instances, a computer controlled automation system can assist the anesthesi-ologist by taking over this tedious and routine task, and thus allowing the doctor toonly be involved with outliers events that are life-threatening. This is similar to the1DOH is a measure of depth of hypnosis only, while DOA is the measure of depth of hypnosisand depth of analgesia.2role of the auto-pilot that takes over the cruising of the plane from the pilot. Thepilot intervenes during unforeseen and critical events only, but otherwise adjustsminor details.Our group has shown that a closed-loop Proportional-Integral-Derivative (PID)controller (known as iControl) can effectively control the DOH, see [10], [52], [54]and [68].To account for the large inter- and intra-variability in the patient’s Pharma-cokinetics (PK) and Pharmacodynamics (PD) response, any controller implementedmust be robust. To be clinically feasible, the controller must perform well despitesurgical stimuli and keep the patient’s DOH at the specified target.A set of performance measures allow to assess the quality of the control ofan anesthetic machine. These measurements need to be of clinical significance toprovide the clinicians with details above the anesthetic state of the patient. If themeasurements could also provide insight of the control architecture to an engineer,then these merits can be used as tuning objectives. Currently, a set of four mea-surements merits, known as Varvel measures [65] have constituted the norm (seeSection 2.4.1.2 Objectives and ScopesThis research started as an assessment of the novel L1 Adaptive Control (L1-AC)[26] which yielded limited feasibility. The focus was then turned to detecting os-cillation caused by a PID controller, and to develop an algorithm to remove theoscillations. A set of metrics are also introduced to quantify the use cases of thecontrol algorithms.The objectives of this thesis is then to 1) assess the application of novel L1-AC[26] as applied to closed-loop control in anesthesia using WAV index as the controlsignal; 2) design an oscillation detection algorithm that can detect multi-periodoscillations in real time; and 3) develop a tuning algorithm that can re-tune the PIDcontroller used in iControl to remove the detected oscillation.Conventional control theory establishes a trade-off between robustness andperformance. For the safety of the patient, the current closed-loop controllers ofanesthesia value robustness over performance. According to its developers, L1-AC3guarantees robustness and its fast adaptation implementation implies a near perfectperformance can be achieved. However, this thesis will show that this structure inits current form does not live up to its promise and therefore cannot be used.The most commonly closed-loop anesthetic controller is a PID system. In [2],it was found that many of these controllers are incorrectly tuned and resulted ininferior performance as well as observed oscillation. The implemented PID con-troller at our facility is tuned to be robust with no oscillation. However in practice,some oscillations have been seen due to patient variability. Partial adaptation thatcan detect and remove oscillation in real-time is discussed. The retuning of thecontroller must consider the following 4 design criteria:1. It must detect the dominant oscillation in real-time.2. It must re-tune the PID controller according to the dominant frequency andensure the new system follows the guidelines of a robust controller.3. It must reject output disturbances and compensate for surgical stimuli.4. It must remove the oscillation and provide adequate performance and set-point response.To sufficiently compare different control schemes, the standard means of mea-suring the performance of closed-loop anesthesia are explored. The Varvel mea-sures were introduced in 1992 for target-controlled-infusion systems but their ad-equacy for closed-loop control is debated. A proposed set of measures introducedby Soltesz et al. [53] is assessed on real clinical data.1.3 Thesis OrganizationThis thesis is organized into 6 Chapters, with this Chapter contributing as one.Chapter 6 provides the closing remarks, conclusions and future works. The sup-porting materials for this thesis are organized into 4 Appendices. The Chapters andAppendices are:4Chapter 2: BackgroundThis Chapter gives a brief overview of current practices in closed-loop anesthesia.It will introduce the currently used monitoring systems for depth of hypnosis, areview of the pharmacokinetics and pharmacodynamics of the propofol drug, aswell as the metrics for performance measures in closed-loop anesthesia.Chapter 3: L1 Adaptive ControlThis Chapter uses an L1-AC to simulate the control of the propofol in patients. Theresults are shown to be in-line with claims that L1-AC fails to provide an adaptivealgorithm and at best behaves as an implementable Linear Time Invariant (LTI)controller. The loss of this adaptivity is mathematically proven.Chapter 4: Real-Time Oscillation DetectionThis Chapter provides an off-line oscillation detection algorithm that is capable ofdetecting multiple oscillation frequencies, along with their fitness 2 and magnitude.An extension to real-time is also introduced where a dominant frequency can bemeasured.Chapter 5: Re-Tuning of a PID ControllerThis Chapter provides a tuning methodology to re-tune a PID controller when anoscillation is detected. The data from the previous Chapter is used to tune thecontroller. Simulation results show the robustness and performance of the re-tunedsystem agree with the current implementation of the iControl system.Appendix A: Propofol PKPD ModelingThis Appendix provides an overview of the propofol Pharmacokinetics/Pharmaco-dynamics (PKPD) model introduced in Bibian [10]. The mathematical model andparameters are also included. These models are used for simulation examples inChapter 3 and 5.2A measure of energy of the oscillation as a percentage of the total energy of the signal5Appendix B: Control Performance in Closed-Loop AnesthesiaThis Appendix provides the mathematical description of Varvel and the proposedalternative measures. The proposed measures are used as tuning objectives inChapter 5.Appendix C: Limiting Behavior of L1 Adaptive ControlThis Appendix provides the mathematical proof for the loss of adaptivity in theL1-AC.Appendix D: Robustness and Performance of iControlThis Appendix provides the complete robustness and performance comparison ofthe iControl and the re-tuned controller from Chapter 5. 44 simulation examplesusing the PKPD models are used in this study.6Chapter 2BackgroundThe field of Biomedical Engineering is on demand. In 2013, this profession wasrated as the #2 with respect to overall satisfaction [16]. Interestingly enough, anes-thesiologists were rated #1 with regard to income [51]. Many of the complex prob-lems faced by clinicians can only be answered through the eyes of an biomedicalengineer. A few examples of engineering solutions for medical diagnostics aregiven.Control engineers traditionally have been focused on the aerospace and processindustries, but have recently applied their knowledge to medical devices. The useof closed-loop control to administer drugs has been shown to improve the qualityand safety as well as reducing the total administrated dosage (see [20] and [2]).More specifically, closed-loop control of drugs delivery has been an active area ofresearch in anesthesia (see [39] and [39] as well as our own research group [52]and references within).A closed-loop anesthesia system measures the DOH from a PEEG signal (suchas BIS or WAV) and controls the infusion rate of the hypnotic drug. The mostcommon hypnotic drug for closed-loop control is intravenous propofol due to itsshort-acting mechanism. The controller can take advantage of the short-acting, fastmetabolic, and fast elimination of the drug and provide a much smoother infusiontitration than an anesthesiologist would be able to do manually. By transferring theresponsibility of these routine tasks to a computer, an anesthesiologist can concen-trate on more vital tasks and the safety of the patient.7In the next few Sections, a review of the required components of a closed-loopcontrol system for anesthesia is provided. First, the methods for measuring theDOH is provided. Second, the modeling of the drug effect in a patient is described.Third, prior attempts at the automatic control of anesthesia are reviewed. Fourth,causes and concerns for oscillations in closed-control are discussed. Finally, areview of the current performance and new proposed measures are provided beforeclosing this Chapter.2.1 Monitoring Depth of HypnosisTo fully control anesthesia, a measurement of DOA1 is required. To this date, nosuch index has been developed. Recent studies have shown that the δ wave inthe EEG signal correlates well with depth of hypnosis. Bispectral, wavelet, timedomain, frequency domain, and evoke potential analysis are a few examples oftechniques applied to the raw EEG signal to extract a single index from it.The most commonly used metric uses the bispectral analysis and is appropri-ately called the Bispectral Index (BIS). A more recent approach uses the waveletanalysis and is called Wavelet-based Anesthetic Value (WAV). In the next two Sub-sections, a brief overview of each method is provided.2.1.1 Bispectral IndexThe BIS monitor was first introduced in 1994 by the Aspect Medical Systems, Incand was marketed as a ”novel measure” of level of consciousness from the EEGsignal [50]. The monitor provides a single index that measures the DOH in thescale of 0 (iso-electric EEG) to 100 (fully awake). The BIS monitor was the firstFDA approved monitoring system [8].BIS is statistically based and empirically derived. A large group of volunteers’EEG were collected and using a proprietary statistical methodology, a model wasfitted to the data. Since the dynamics of the system are unknown, it is difficult todesign an optimal controller with this output signal. Moreover, the BIS Monitor hasa time delay between the changes in the patient’s anesthetic state and the changes in1DOA is a measure of both the hypnotic and the analgesic states. DOH on the other hand whichis a measure of the depth of hypnosis only.8the BIS value [8]. These two problems together are the motivation to have anothermonitor whose dynamics are known, and whose response has no delay.2.1.2 Wavelet-based Anesthetic ValueThe Wavelet-based Anesthetic Value (WAV) is an alternative to BIS that is basedon the wavelet decomposition of the raw EEG signal. Proposed by Bibian et. al[11], this hypnotic monitoring value correlates well with the BIS. The dynamicsof the system are described by a simple transfer function 1/(8s2+1) and respondsmuch faster to the changes in anesthetic state than the BIS. There is minimal to nodelay in the signal response [69]. These two advantages of the WAV index make itappealing to be used as a control signal. It has been shown to lead to an improvedperformance in closed-loop control of anesthesia [10].2.2 Drug Effect ModelingThe physiological effect of an administrated drug on a patient is typically describedwith twomodels: the pharmacokinetic and the pharmacodynamic model. The phar-macokinetic model (PK) relate the administrated drug dosage to the drug plasmaconcentration. The pharmacodynamic model (PD) then relate the drug plasmaconcentration to the physiological effect. These models describe the distribution,metabolism, and the clearance of the drug in the body to the resulting physiologicaleffect.In this literature review, an overview of the PK and the PD as described in [10]is introduced and briefly reviewed. A more detailed discussion can be found in[43] and [10].2.2.1 Pharmacokinetics of PropofolPharmacokinetic model represents the drug uptake, distribution and elimination.The mathematical model then relates the infusion rate to the drug plasma concen-tration. The first significant investigation to study effect of sampling site (venous vsarterial) and the method of drug administration (bolus vs infusion) was conductedin 1998 by Schnider et al. [47]. A more recent study in 2000 was conducted bySchu¨ttler and Ihmsen [48] and is discussed in this Section.9The Schu¨ttler and Ihmsen study was a large collaboration between 5 institutes,where 4,112 samples from 270 individuals (150 men, 120 women) of the ages 2-88years, body weights of 12-100 Kg were studied. The objectives of this study was:1. Estimate the pharmacokinetics of propofol with respect to the covariates age,body weight, and gender.2. Evaluate the inter- and intra-patient variability.3. Study the effect of the mode of administration (bolus vs infusion).4. Study the effect of the sampling site (venous vs arterial).The result showed that the pharmacokinetics of propofol is best described by a 3-compartment model, see Figure A.1 in the Appendix A. Weight was determinedto be the most prominent factor; age, gender and mode of administration werealso positively correlated. The sample site had a little influence. The intra-patientvariability in this study was found to be less than 20%.The mathematics of the modeling can be found in Appendix A. The PK param-eters can be found in Table A.3.2.2.2 Pharmacodynamics of PropofolPharmacodynamic model is the observed effect of the drug as a function of the drugplasma concentration. A single drug interacts with multiple organs in the body andhas multiple pharmacological effects. Here, the model for the depth of hypnosisfrom the EEG is considered.There have been a limited number of studies on quantifying the effect of propo-fol on the EEG. A detailed discussion on these studies can be found in [10]. Manyof the studies show large inter-patient variability in the PD model [10]. Moreover,most of these studies derive thePK model for the BIS.In Bibian [10], the dynamics of propofol vs the WAV was modeled throughthe analysis of 44 patients. Using least-squares identification, Bibian estimateda PD model consisting of a Hill function followed by a first-order time delayedtransfer function. The Hill function models the drug-receptor binding interaction.The first-order transfer function was proposed by Sheiner et al. [49] to model the10temporal aspect of the pharmacodynamics. The time delay was added to representthe arm-to-brain circulation time.The details of the modeling can be found in Appendix A. The PD parameterscan be found in Table A.3.2.3 Automatic Control of AnesthesiaDuring surgery, the anesthetic and opioid titration are constantly adjusted to pre-vent under- and over-shoot of the drug plasma concentration and to keep the anes-thetic state constant. An automated system that regulates the administration ofthese drugs thus seems appealing to the anesthesiologists.The idea of closed-loop control of anesthesia has been investigated for half acentury now. The performance and robustness of these controllers depend stronglyon the mathematical model (PKPD) of the patients, the monitoring devices (BISor WAV) as well as the tuning of the controller itself. The ideal controller shouldmeasure the three functional states hypnosis, analgesic and paralysis and regu-late the administration of hypnotic, opioid and neuromuscular drugs. This Multi-Input/Multi-Output (MIMO) system is currently not available due to the limitationof monitoring systems as well as the mathematical models that govern the drugadministration.In recent years, the number of published studies on this field has increasedsignificantly. A literature review on the current attempts on closed-loop control ofhypnosis is provided next.2.3.1 Closed-Loop Control: A ReviewIn the following reviews, adequate anesthesia is considered as the BIS or WAV inthe range of 40-60 [39].In 1999, Frei et al. [21] used a Model-Predictive-Controller (MPC) to controlthe Mean Arterial Pressure (MAP) using the inhaled drug, isoflurane. The studywas performed on over 100 subjects and proved a better performance than manualcontrol. The authors initially designed a PID-like Fuzzy controller. However, thecontroller was unable to account for respiratory dynamics under low flow condi-tions. The MPC model was implemented due to this inadequacy.11In 2001 Struys et al. [56] compared the performance of an adaptive model-based control guided by the patient’s BIS to manually control anesthesia using in-travenous administration propofol. The study was conducted on 20 female subjectsaged 34-50 years undergoing gynecologic laparotomy. Subjects were randomizedwith half under closed-loop control and the other half were manually controlled.The study found that the manually controlled patient had a shorter induction time.The closed-loop controlled patients had a better maintenance performance as val-idated by Varvel measures, as well as a reduced recovery time. No details wereprovided on the controller structure.In 2002, Absalom et al. [2] used a PID controller guided by the patient’s BISusing intravenous administration of propofol in 10 patients undergoing elective hipor knee surgery. Performance was validated using Varvel measures. The authorsreported clinically adequate anesthesia in 9 out of the 10 patients. Three of thepatients’ BIS oscillated around the set-point, although none of these cases showeda sign of inadequate anesthesia. The controller used was from another study byKenny et al. [31] where a PID controller was guided by the auditory evoked poten-tial.In a follow up study by Absalom et al. [1] in 2003, a revised PID controller wasused. In this study, 20 adult patients (12 female, 8 males) undergoing body surfacesurgery were enrolled. The patients were initially controlled with an open-looptarget-controlled-infusion. Once the anesthesia was clinically adequate, the systemwas switched to the revised PID. All 20 patients reported a clinically adequateanesthesia. There was one patient with oscillation.A more interesting study in 2004 by Locher et al. [42] used a cascade structurewith an outer Proportional-Integral (PI) and inner model-based state feedback con-troller guided by the patient’s BIS using isoflurane. The study was performed on23 patients undergoing decompressive spinal surgery who were randomized intoclosed-loop or manual control. The study had two conclusions: 1) the closed-loopcontrol significantly outperformed the manual mode and, 2) the closed-loop controladministrated less total drug and faster wake-up time.In 2006 study by Liu et al. [39], 164 patients undergoing elective minor ormajor surgery were randomized into closed-loop and manual target control infusiongroups. The closed-loop system was an empirically tuned PID controller. The12patient’s BIS was used as the control signal and propofol and remifentanil wereadministrated intravenously. Propofol consumption was lower in the closed-loopgroup, but the induction time was longer. Adequate anesthesia was significantlybetter in the closed-loop group. Recovery time was also shorter for the closed-loopgroup.Finally, L1- output feedback adaptive control was used in 2011 by Ralph et al.[46] in a simulation study using the BIS as the control signal and isoflurane as thehypnotic drug . Seven PKPD models were reconstructed using clinical trial data.A controller was designed based on one of the identified PKPD models. The samecontroller was then applied to the other six models. The result showed adequatereference tracking.2.3.2 Closed-Loop Control: OscillationOscillation in closed-loop control can occur as a combination of any of the follow-ing: 1) marginally stable control loops (due to aggressive control tuning or changesin process gain/phase/time delay); 2) external disturbances; 3) stiction in controlvalve [9]. If the controller is improperly tuned, the oscillation can cause instability.In Chapter 4, a detailed root cause analysis of oscillation is provided.In [2], [1] and recently in our own work [57], oscillation in the patient’s BIS andWAV was detected. Therefore, it is essential that oscillation be detected in real-timeto both warn the anesthesiologist and to remove it by retuning the controller.One of the first attempts at oscillation detection was by Ha¨gglund [25] in 1995.His method computed the Integrated Absolute Error (IAE) between consecutivezero-crossings of the error. When oscillation occurs, the absolute error and thetime between consecutive zero-crossings increase, leading to a higher IAE. Bycounting the instances of IAE larger than a threshold in a given period of time,oscillation can be detected. This method, however, can fail to detect oscillationswhen multiple frequencies exist. Moreover, it cannot determine all the differentoscillation frequencies in a signal.Wang et. al. [67] review a large set of different algorithms. Auto-correlationfunction, Discrete Wavelet Transform (DWT) method, empirical mode decomposi-tion, and Discrete Fourier Transform (DFT) are to just name a few methods applied13since 1995. Many have limitations. DFT has the disadvantage that the default rect-angular window only provides good energy compaction for frequencies that arewhole fractions of the sampling frequency Fs. DWT has the disadvantage that itcan be computationally expensive. Wang et al. provide a new method based on theDiscrete Cosine Transform (DCT) that overcomes all the shortcomings. It is fast,independent of the sampling frequency, and can decompose oscillation into all ofits period components.An improved representation of Wang’s DCT method is provided in Chapter 4.The period, magnitude and fitness of the dominant oscillation is detected in real-time. In Chapter 5, these information are used to re-tune the PID controller.2.3.3 Closed-Loop Control: Adaptive vs PIDA closed-loop control system can be divided into two generic types: adaptive andnon-adaptive (classical). An adaptive controller is a system whose parameters canadapt continuously to the plant it is trying to control [28]. This adaptation canbe in response to initial uncertainty in the plant or the change in the plant itself(for instance, an aircraft loses weight due to fuel consumption). The non-adaptive(classical) controller is a system in which the controller is not changed once it isimplemented. Adaptive controllers in theory can provide better performance androbustness as they adapt to the particular plant. From the reviews in 2.3.1, adaptivecontrol is still not well understood for use in closed-loop anesthesia.In this paper, the newly introduced L1 Adaptive Control (L1-AC) is reviewed[26]. The Proportional-Integral-Derivative (PID) controller, which accounts forabout 90% of all the controllers used in the industry is also considered [5]. Ourresearch group currently uses a PID closed-loop control system for controlling theDOH (see [10] and [52]). In Chapter 5, a tuning algorithm is introduced to auto-matically re-tune the controller in the presence of an oscillation.2.4 Performance of Closed-Loop AnesthesiaA set of four performance measures (MDAPE, MDPE, Divergence, and Wobble),proposed by Varvel et al. [1], have constituted the standard means of assessing per-formance in closed-loop anesthesia. Varvel measures were developed for Target-14Controlled Infusion (TCI) anesthesia systems; they were not developed for EEG-guided closed-loop controllers. These measures are not accepted within the controlcommunity and cannot be used as control tuning parameters. Moreover, they onlyaccount for the maintenance phase of anesthesia. Varvel measures are based onthe median of the relative error. There is no distinction between artifacts, noise,and momentary large errors. There is also no penalty for outliers when adoptedfor EEG-guided DOH control and the metrics are not normalized with respect toduration of the case.Soltesz et al. [53] proposed an alternate set of measures. The key features ofthese measures are : 1) wide acceptance in control community; 2) consideration ofclinical feasibility; 3) separation of metrics for induction, maintenance and emer-gence phases of anesthesia. For the induction phase, Induction Phase Duration (ID)and Percent Overshoot (OS) are proposed. For the maintenance phase, IntegratedError (IE), Integrated Absolute Error (IAE), Variability Index (VI) and percentageof time outside the adequate range are proposed. For the emergence phase, Emer-gence Phase Rise Time (ER) is proposed. The mathematical details of Varvel andthe proposed measures can be found in Appendix B.We analyzed 63 clinical cases that were collected from a study on closed-loopcontrol DOH using the NeuroSense monitor [57]. The study was approved by UBCChildrens and Womens Research Ethic Board (H10-01174), Vancouver, Canada[61]. The population included 32 women, 31 men between the ages of 6-17 yearsold, body weight of 14.5 - 70 Kg, and height of 106 - 182 cm. The propose mea-sures provided more insight about the control performance, as discussed in Ap-pendix B.There are certain scenarios where Varvel measures can be misleading. DOHvalues in the set-point ±10 range are considered adequate. Maintenance phaseslike the one in Figure 2.1 should be more desirable than the one in Figure 2.2;using the error metric on the median (Varvel) has the opposite effect.The IE punishes outliers linearly while the median-based error metric MedianPerformance Error (MDPE) filters out outliers. The DOH in Figure 2.3 is clearlymore negatively biased than in Figure 2.4. The MDPE metric concludes the oppo-site, while IE reflects this bias. Furthermore, IE is used as minimization criterionin existing controller synthesis strategies.150 5 10 15 20 25 30 35 40 45Time (min)020406080100DOH [WAVcns]aiWaveMaintenanceiWAVSet0 5 10 15 20 25 30 35 40 450100020003000Propofol infusion[mcg/kg/min] iPropoDose0 5 10 15 20 25 30 35 40 45Time (min)02468Plasma concentration PlasmaConciCE - iControliCPlasma - iControlFigure 2.1: MDAPE vs IAE for a systematic small error.Varvel metrics do not provide any measures for the induction phase. Lengthof ID affects the initial performance of the maintenance phase. The long ID resultsin a large total initial drug dosage, and an excessive overshoot of the DOH. ShortID results in low plasma concentration and signals the possibility of rapid rising inDOH. This information is available in the proposed measure as seen in Figure 2.5and 2.6.160 10 20 30 40 50 60 70 80 90Time (min)020406080100DOH [WAVcns]biWaveMaintenanceiWAVSet0 10 20 30 40 50 60 70 80 900500100015002000Propofol infusion[mcg/kg/min] iPropoDose0 10 20 30 40 50 60 70 80 90Time (min)02468Plasma concentrationPlasmaConciCE - iControliCPlasma - iControlFigure 2.2: MDAPE vs IAE for a sporadic error.170 10 20 30 40 50 60Time (min)020406080100DOH [WAVcns]aiWaveMaintenanceiWAVSet0 10 20 30 40 50 60050010001500Propofol infusion[mcg/kg/min] iPropoDose0 10 20 30 40 50 60Time (min)02468Plasma concentrationPlasmaConciCE - iControliCPlasma - iControlFigure 2.3: DOH is clearly negatively biased.180 10 20 30 40 50 60 70 80Time (min)020406080100DOH [WAVcns]biWaveMaintenanceiWAVSet0 10 20 30 40 50 60 70 800500100015002000Propofol infusion[mcg/kg/min] iPropoDose0 10 20 30 40 50 60 70 80Time (min)02468Plasma concentrationPlasmaConciCE - iControliCPlasma - iControlFigure 2.4: DOH is less biased.190 10 20 30 40 50 60 70Time (min)020406080100DOH [WAVcns]aiWaveMaintenanceiWAVSet0 10 20 30 40 50 60 70050010001500Propofol infusion[mcg/kg/min] iPropoDose0 10 20 30 40 50 60 70Time (min)02468Plasma concentrationPlasmaConciCE - iControliCPlasma - iControlFigure 2.5: Small ID of 3.1 min translates to a small overshoot of 12.8%.200 10 20 30 40 50Time (min)020406080100DOH [WAVcns]biWaveMaintenanceiWAVSet0 10 20 30 40 5001000200030004000Propofol infusion[mcg/kg/min] iPropoDose0 10 20 30 40 50Time (min)02468Plasma concentrationPlasmaConciCE - iControliCPlasma - iControlFigure 2.6: The large ID of 5.6 min translates to a larger overshoot of 48.9%and a longer DOH settling time to set-point.21Chapter 3L1 Adaptive ControlThere is a high inter-patient variability in the effect of the hypnotic drug on theirDOH. A closed-loop system that controls the drug administration needs to guaran-tee robustness and performance. The rate at which an adaptive controller adapts tothe patient is called the adaptation gain, Γ. It is a well known fact that high-gain inthe feedback loop of a controller leads to amplified high frequency components inthe control signal, reduction in phase margin, and loss of robustness [28]. Numer-ous authors have tried to introduce the concept of fast adaptivity with robustnessas core to the control design (see for instance [28] and [35]) as classical robustnessconcepts are not applicable.In an adaptive controller, Γ shows up in the adaptation law: a nonlinear dy-namic system that identifies a known parameter related to the uncertainty of theplant. As the adaptation gain increases, the rate at which the unknown parameteris identified, also increases.Classically, there has been a trade-off between robustness and performance:as one increases, the other decreases. In adaptive control the same trade-off ex-ists: increasing the adaptive gain will improve the performance (by increasing theadaptivity) at the cost of reducing the robustness of the system.This Chapter discusses L1 adaptive control as introduced in the book ”L1Adaptive Control Theory: Guaranteed Robustness with Fast Adaptation” [26]. Theauthors suggest that through their unique control structure, the adaptivity is decou-pled from the robustness, i.e. one can increase the adaptive gain to arbitrarily large22values without effecting robustness. Robustness is then guaranteed through clas-sical methodology. The L1-AC defines an unimplementable reference model andguarantees the difference between this reference model’s output and the patientmodel’s output decreases as the adaptation gain increases.At the start of this research, there had been some doubts on the validity andclaims of the theory [29]. In [45] it is shown that high-gain leads to system in-stability. In [13], the L1-AC shows inferior performance as compared to otherwell established adaptive controllers. It has also been seen in simulation that theadaptivity is lost as the gain increases [63]. Recently, the authors of L1-AC haveproposed four different adaptation formulations which all lead to the exact sameperformance bounds [64]. One of these formulations is in fact LTI for all adapta-tion gains. All of these research however, have not proved or disproved the L1-AC;they have shown examples where the stability and the adpativity are lost.The structure of this chapter is as follows: in Section 3.2, the L1 AdaptiveControl structure is introduced. The reference system and performance bounds arediscussed in Section 3.3. Simulation examples are provided in Section 3.4. Lossof adaptivity is discussed in Section 3.5.3.1 ContributionThis chapter will review the claims made about L1-AC. First, it is shown that in-creasing the gain results in a loss of adaptivity. Second, it is shown that the limitingbehavior (the case with Γ going to infinity) of the L1-AC can be achieved throughan implementable, non-adaptive LTI controller. Finally, the loss of adaptivity ismathematically shown to be the direct result of inversion of the estimation loop asthe gain increases. An example at the end of the chapter shows how a series ofadaptive, non-adaptive, dynamic, static, linear and nonlinear laws that all lead tothe exact same limiting controller as the adaptation gain increases.233.2 The L1 Adaptive Control3.2.1 Problem FormulationConsider the following state-feedback dynamic controller G(s) within the L1-ACarchitecture (see Chapter 2.2 of [26]):x˙(t) = Amx(t)+b(ω(t)u(t)+θT (t)x(t)+σ(t)), x(0) = x0,y(t) = cT x(t),(3.1)where x(t) ∈ Rn is the measured state of the system; u(t) ∈ R is the control input;y(t) ∈ R is the output; b,c ∈ Rn are assumed known constant vectors; Am is an× n Hurwitz matrix corresponding to the desired closed-loop dynamics; ω ∈ Ris an unknown constant but with known sign; θT (t) ∈ Rn is a vector of unknownparameters; and σ(t) ∈R models input disturbances. The dynamics of the desiredmodel M(s) are given by:x˙m(t) = Amxm(t)− kgbr(t), xm(0) = x0,ym(t) = cT xm(t),(3.2)where kg ,−1/(cTA−1m b) and r(t) is the reference signal.Assumption 1. Boundedness of the unknown parameters: Let the unknownparameters θ(t) and σ(t) be bounded as:θ(t) ∈Θ, |σ(t)| ≤ Σ,where Θ and Σ are both known bounds of θ(t) and σ(t) respectively. Furthermore,let the lower and upper bound of ω(t) be known:ωlb ≤ ω(t)≤ ωub, ∀t ≥ 0.These bounds need to be chosen from prior knowledge of the inter-variability inthe patients’ models and the expected input disturbances.243.2.2 State PredictorThe state predictor in the L1-AC is given by:˙ˆx(t) = Amxˆ(t)+b(ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ (t)), xˆ(0) = x0yˆ(t) = cT xˆ(t).(3.3)The predictor has the same structure as 3.1; the unknown parameters ω(t),θ(t),and σ(t), are replaced by their estimates ωˆ(t), θˆ (t), and σˆ(t).The adaptation laws for the three unknown parameters are given by the follow-ing projection operator [34]:˙ˆθ(t) =−Γ ·Proj(θˆ (t),−x˜TPbx(t)), θˆ (t) = θ0,˙ˆσ(t) =−Γ ·Proj(σˆ(t),−x˜TPb), σˆ(t) = σ0,˙ˆω(t) =−Γ ·Proj(ωˆ(t),−x˜TPbu(t)), ωˆ(t) = ω0,(3.4)where x˜(t) = xˆ(t)− x(t), Γ ∈ R+ is the adaptation gain, and P = PT > 0 is thesolution of the algebraic Lyapunov equation ATmP+PAm = −Q for arbitrary Q =QT > 0.Finally, the L1-AC signal is defined as:u(s) =−kD(s)(ηˆ(s)− kgr(s)), (3.5)where r(s) and ηˆ(s) are the Laplace transforms of r(t) and ηˆ(t) respectively andηˆ(t), ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ(t). (3.6)k > 0 is a feedback gain and D(s) is a strictly proper transfer function such thatthey lead to a strictly proper stable filterC(s):C(s) =ωkD(s)1+ωkD(s). (3.7)The controller is shown in Figure 3.1.25kg kD(s)ηˆ(t) = ωˆ(t)u(t) + θˆT (t)x(t) + σˆ(t)x˙(t) = Amx(t) + b(ωu(t) + θTx(t) + σ(t))y(t) = cTx(t)x˙(t) = Amxˆ(t) + b(ωˆu(t) + θˆTx(t) + σˆ(t))yˆ(t) = cT xˆ(t)˙ˆθ(t) = −ΓProj(θˆ(t),−x˜TPbx(t))˙ˆσ(t) = −ΓProj(σˆ(t),−x˜TPb)˙ˆω(t) = −ΓProj(ωˆ(t),−x˜TPbu(t))r +−+−Figure 3.1: The original L1-AC block diagram as it appears in [26].3.2.3 L1-norm Stability ConditionThe L1-AC is subject to the following L1-norm condition:‖L(s)‖L1T < 1 (3.8)where L(s) and T are computed as:T ,maxθ∈Θ‖θ‖, H(s) = (sI−Am)−1b, L(s) = H(s)(1−C(s)). (3.9)If the condition 3.8 in presence of Assumption 1 is satisfied, then the L1-AC isguaranteed to be stable. In calculating the L1-norm, C(s) depends on the unknownparameter ω , which should be chosen as the worst expected case.The claim for the L1-AC is as follows: compute a gain k and a filter D(s) suchthat for the worst case ω , theL1-norm stability condition holds. This will guaranteethe robustness of the system. Then increase Γ as high as computationally possibleto increase the performance. The filter kD(s)will act as the decoupler of robustnessand performance trade-off.3.3 Achievable Performance BoundThe controller cannot achieve the desired system dynamics M(s) as a direct resultof the introduction of the low-pass filter kD(s) in the control loop. Instead, a refer-ence system Gre f (s) is introduced and the control performance of the system G(s)is compared to the performance of the reference system Gre f (s). The reference26system is defined as:x˙re f (t) = Amxre f (t)+b(ω(t)ure f (t)+θT (t)xre f (t)+σ(t)),yre f (t) = cT xre f (t),ure f (s) =C(s)ω(kgr(s)−ηre f (s)),(3.10)where ηre f (s) is the Laplace transform of ηre f (t), θT (t)xre f (t)+σ(t). This con-troller is not implementable as it depends on the system unknowns ω(t), θ(t) andσ(t).3.3.1 Reference ControllerAssume an initial condition xre f (0) = x0 = 0. The reference system can be writtenas:xre f = H(s)(ωure f +ηre f ), (3.11)where the Laplace operator s is intentionally excluded from the signals to simplifythe calculation. The above equation cannot be solved for ηre f (s) since H(s) isnot invertible. Multiplying the equation above by (Pb)T makes the (Pb)TH(s)invertible and ηre f (s) can be solved for:ηre f =(Pb)T xre f(Pb)TH(s)−ωure f . (3.12)Substituting ηre f from 3.12 into 3.10 leads to:ure f =C(s)ω(kgr− (Pb)T xre f(Pb)TH(s)+ωure f). (3.13)Isolating for ure f results in:ure f =C(s)ω(1−C(s))(kgr− (Pb)T xre f(Pb)TH(s)). (3.14)27Taking into account the definition of C(s) from 3.7, the equation above can besimplified to:ure f (s) =kD(s)(Pb)TH(s)(Pb)T(H(s)kgr(s)− xre f (s)). (3.15)Even though the reference system 3.10 is not implementable, this control signalis implementable since it does not depend on the system’s unknown parameters(θ(s), σ(s) and ω(s)). The control signal corresponds to an implementable LTIcontroller whose dynamics only depends on the filter kD(s), the desired modeldynamics H(s) and the solution to the Lyapunov equation, P. Still, L1-AC claimsto achieve this non-adaptive behavior as the Γ → ∞ in an adaptive structure.3.3.2 Control PerformanceThe L1-AC structure guarantees the following bounds:Lemma 1 (From [26]). Let the system G(s) be controlled by the L1 adaptive con-troller from Section 3.2. Assume the L1-norm stability condition of 3.8 is satisfiedand the bounds of Assumption 1 are met. Assume the reference system 3.10 isstable, i.e. the system Gre f (s) is stabilized through the LTI reference controllerηre f (s) from equation 3.15. Then, the system state x(t) and control input u(t) areuniformly bounded:‖xre f − x‖L∞ 6γ1√Γ, ‖ure f −u‖L∞ 6γ2√Γ, (3.16)where γ1 and γ2 are constants. The full details of the calculation is provided inChapter 2.2, pages 40-41 of the book [26].The details of the calculation of the Lemma is not important. Rather, the inverserelationship of xre f − x and ure f − u to the adaptation gain Γ is significant. Theseperformance criteria motivate the use of high adaptation gain.280 500 1000 1500 2000 2500 300000.10.20.30.40.50.60.70.80.91x 10−4 Step ResponseTime (seconds)AmplitudeModelPatientFigure 3.2: The step response of patient #7 and the predictor model3.4 Case Studies of the L1 ControllerIn this Section, the PKPD models from Appendix A are used for two case studies.Since the PK model depends on the demographic information of the patient only,the predictor’s PK model is chosen as the patient’s PK. The PD is chosen as theaverage of all 44 cases, with the delay set to zero. The predictor’s model PKPDis also slightly modified to have a faster response than the patient’s. While thispredictor model is not ideal, it suffices for the purpose of showing the limitation ofthe L1-AC. Figure 3.2 shows the step response of the predictor model and a patient(case #7).Let D(s) = 110s2+sand k = 1 so that the filterC(s) is given as:C(s) =ω10s2+ s+ω. (3.17)This satisfies the L1-norm condition of 3.8: the value of ‖L(s)‖L1T is between0.0168 and 0.0673 for all 44 models. The bounds of the unknown parameters areΘ = 1, Σ = 1 and −5 ≤ ω ≤ 5 and are chosen as the maximum of the 44 PKPDmodels.29The system is simulated for 3 different adaptation gains 1, 5×103 and 1×105.The results of the simulation are shown in Figure 3.3. The patient’s output for thethree different gains is shown in the upper plot, while the absolute error betweenthe patient’s output and the reference’s output is shown in the lower plot.For an adaptive algorithm, it is expected that the error between the referenceand the patient decreases throughout the case. For the small gain of Γ = 1 the resultclearly shows a decrease in error as the case progresses. However, for the highergain, specially Γ = 1× 105, while the initial error is lower, it does not improve asthe case progresses, i.e. the adaptivity of the system is lost. An almost identicalbehavior is shown for another case (patient model #2) as shown in Figure 3.4. Inthis case, even for the intermediate gain Γ = 5×103 the adaptivity is almost lost asthe error does not decrease as time continues.The claim of achieving the non-implementable reference model in an adaptivearchitecture is invalid, since adaptivity is lost. The exact same performance canbe achieved using the implementable reference controller defined in 3.15. Thiscontroller is LTI, does not depend on any unknown parameters, and provides thesame performance.30Figure 3.3: Simulated system output for the patient #7 controlled by theL1 controller. The upper plot shows the patient’s output for the threeadaptation gains. The reference output is shown in the thick green line.The lower plot shows the absolute error of patient’s output to the refer-ence’s output.31Figure 3.4: Simulated system output for the patient #2 controlled by theL1 controller. The upper plot shows the patient’s output for the threeadaptation gains. The reference output is shown in the thick green line.The lower plot shows the absolute error of patient’s output to the refer-ence’s output.32vΓf(·)u+−Figure 3.5: Simple feedback for a nonlinear function f (·)3.5 Loss of AdaptivityThe previous section provided two case study examples that showed the loss ofadaptivity of the L1-AC as the gain increased. This is in-line with the previ-ous claims that this controller does not achieve a better performance than an im-plantable LTI controller (see [13], [60] and [29]). The loss of this adaptivity is dueto the high adaptation gain, and will be discussed in this Section [63].In a simple feedback gain system, shown in Figure 3.5, it is easy to show thatincreasing the gain Γ leads to the inversion of the nonlinearity dynamic f (·) [23].Straight forward calculation shows that u= Γ(v− f (·)u). Solving for u gives:u=Γ1+ f (·)Γv. (3.18)When the gain Γ increases to infinity, the system dynamic are inverted, i.e.ulim = limΓ→∞Γ1+ f (·)Γv=Γf (·)Γv=1f (·)v= f (v)−1. (3.19)In the L1 controller Figure 3.1, there is a similar high-gain feedback over thepredictor’s nonlinearity and it is expected that increasing the gain will also invertthe predictor, albeit the loop over the signal η(s) makes this observation challeng-ing. However, the L1 architecture is also inverted as the gain goes to infinity. Thedetails of this calculation is available in Appendix C for reference.Figure 3.6 shows the linearized system around some equilibrium points θTQ ,ωTQ , and σTQ . Here, G(s) and H(s) refer to the dynamics of the patient and thepredictor model. Fx, Fxˆ and Fu are some LTI functions related to the linearizedcomponents of the adaptation laws as shown in C.14.After linearizing the projection adaptive laws of 3.4, the transfer function be-33kg kD(s) G(s)r +−Fx(s)1−Fxˆ(s)H(s)+kD(s)Fu(s)u xFigure 3.6: Final form of the linearized adaptation laws for L1 adaptive con-trol with generic adaptation laws.tween u and x can then be written as:u(s) =−kD(s)[Γs(xTQxQ+u2Q+1)(Pb)T + θˆTQ]1+ Γs(xTQxQ+u2Q+1)(Pb)TH(s)+ kD(s)ωQx(s), (3.20)where xQ and uQ are the equilibrium values of the signal x and u. Here, it isassumed that at equilibrium the states x and xˆ are the same (i.e. x˜= xˆ−x= 0). Thephrase (xTQxQ+u2Q+1) is just a constant. The effect of increasing the gain is nowabundantly clear - the limiting controller is:ulim(s) = limΓ→∞−kD(s)[Γs(xTQxQ+u2Q+1)(Pb)T + θˆTQ]1+ Γs(xTQxQ+u2Q+1)(Pb)TH(s)+ kD(s)ωQx(s)=−kD(s)(Pb)T(Pb)TH(s)x(s).(3.21)This limiting controller has the exact same control signal as the reference signal3.15 that was derived from the reference system. Note that unlike the referencesignal, this result is not derived from the mathematical description of the plantG(s) and holds true for any kind of plant G(s), given the stability condition is met.This limiting controller derivation indicates that for high adaptation gains:• The integral effect in the adaptation law is canceled.• The predictor is inverted.• The effect of nonlinearity of the controller is canceled.• The choice of the equilibrium point xQ, xˆQ and uQ does not affect the con-troller. Moreover, the exact formulation of the adaptive laws is irrelevant.34In a follow up paper by the authors of the L1-AC, it was shown that four dif-ferent adaptations laws resulted in similar performance bounds [64]. Above resultexplains why this is so. In fact, the result above suggests that any law, adaptiveor non-adaptive, linear or nonlinear, dynamic or static, will result in the same per-formance bounds, provided that the stability condition and bounds are satisfied.An example will now follow: Six laws are shown to provide the exact same per-formance bound as the gain increases. The example used is the one used by theauthors’ themselves in [15] and [14].3.5.1 Simple Example of Loss of AdaptivityConsider the system 3.1 with Am = 0 1−1 −1.4, b=0.51, c= [1 0], ω = 1,θ =22 and σ = 1. This corresponds to the system G(s) = 0.5s+1.7s2−1.6s−1.4 and thepredictor (Pb)TH(s):(Pb)TH(s) =N(s)D(s)N(s) = 1.6s7+8.2s6+20.5s5+31.4s4+31.4s3+20.5s2+8.2s+1.6,D(s) = s8+5.6s7+15.8s6+27.8s5+33.4s4+27.8s3+15.8s2+5.6s+1,(3.22)where P solves the Lyapunov equation with Q = I2×2. With k = 60 and D(s) =1s(1+0.1s) , the filter C(s) is defined as:C(s) =600.1s2+ s+60. (3.23)The L1-norm condition with the choice of this filter is 0.0858. The following 6adaptation laws are considered:35A) The dynamic L1 projection adaptation law:˙ˆθ(t) =−Γ ·Proj(θˆ (t),−x˜TPbx(t))˙ˆσ(t) =−Γ ·Proj(σˆ(t),−x˜TPb)˙ˆω(t) =−Γ ·Proj(ωˆ(t),−x˜TPbu(t))ηˆ(t) = ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ (t).B) Another dynamic nonlinear adaptation law:˙ˆθ(t) =−Γ · x˜TPb‖x(t)‖˙ˆσ(t) =−Γ · x˜TPb˙ˆω(t) =−Γ · x˜TPb|u(t)|ηˆ(t) = ωˆ(t)|u(t)|+ θˆT (t)x‖(t)‖+ σˆ (t).C) The static L1 projection adaptation law:θˆ (t) =−Γ ·Proj(θˆ (t),−x˜TPbx(t))σˆ(t) =−Γ ·Proj(σˆ(t),−x˜TPb)ωˆ(t) =−Γ ·Proj(ωˆ(t),−x˜TPbu(t))ηˆ(t) = ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ (t).D) A higher order nonlinear static adaptation law:ηˆ(t) =−Γ · [xˆ3(t)− x3(t)]TPb.E) A linear adaptation law as used in [64]:ηˆ(t) =−Γ · x˜TPb.36F) A switching adaptation law as used in [33]:θˆ (t) =−∆θ · sgn(dzεθ(xˆ(t)− x(t))T ·Pbx(t))σˆ(t) =−∆σ · sgn(dzεσ(xˆ(t)− x(t))T ·Pb)ωˆ(t) =−∆ω · sgn(dzεω(xˆ(t)− x(t))T ·Pbu(t))ηˆ(t) = ω(t)u(t)+ θˆT (t)x(t)+ σˆ(t),where sgn(·) is the sign function; dz(·) is the dead-zone function; εθ ∈ R+,εσ ∈ R+, and εω ∈ R+ are the dead-zone intervals; ∆θ , ∆σ , and ∆ω are thebounds of the unknown parameters from Assumption 1.Laws A), E), and F) are 3 of the 4 laws which the authors of the L1-AC havethemselves introduced as alternatives and have shown that they achieve the sameperformance bounds [64] (the forth law is only applicable for a different class ofL1 controllers). Also, only Laws A) and B) are dynamic and adaptive; the other4 laws are static and non-adaptive. For the Laws A) to E), the gains used areΓ = 1,1e3,1e4. For Law F), the dead-zone intervals used are dz = 1,0.1,0.01.Figures 3.7 to 3.12 show the simulation results for the 6 laws.For Laws A) and B), the error between the reference output and the plant outputdecreases as the case progresses for low adaptation gain, as expected of an adaptivecontroller. However, as the gain increases, while the initial error is lower, it doesnot improve over time, suggesting that the system has lost its adaptivity.In Law C) the integral action of the projection law is removed and the adapta-tion law is static. None of the gains result in an adaptive controller, yet the limitingcase is identical to the original projection law. This law shows that a static, non-adaptive ”adaptation law” provides the same plant output.Law D) is an unnecessary and computationally heavy law that is not practicalfor any application, and is only intended for demonstration purposes. For lowadaptation gain, the system has sustained oscillation. However, for higher gains,again the system approximates the LTI reference model and has the same output asthe original adaptive projection law.Law E) is a very simple error feedback, non-adaptive law that was suggested asan alternative solution by the authors of the L1-AC [64]. For low adaptation gain,37Figure 3.7: Simulated output for the plant G(s) controlled by the L1-ACfor Law A). The top figure is the output of the plant, with the thickgreen line being the output of the reference. The lower plot is the abso-lute difference of model output and reference output. The controller isadaptive for low adaptation gain, however it becomes static for highergains.the system has a steady-state error; since there is no integral action or adaptivity,the error does not reduce as the case progresses. The steady-state error reduces asthe adaptation gain increases.Law F) is another law suggested by the authors of L1-AC [33]. This law isthe most computationally exhaustive algorithm and has no improvements over theother laws. It is again non-adaptive and the plant’s model approximates the refer-ence model as the gain increases.380 5 10 15 20 25 30 35 40 45 50Time (sec)0123456789OutputLaw BReference = 50 = 1E3 = 1E40 5 10 15 20 25 30 35 40 45 50Time (s)10-1010-810-610-410-2100|Error| = 50 = 1E3 = 1E4Figure 3.8: Simulated output for the plant G(s) controlled by the L1-ACfor Law B). The top figure is the output of the plant, with the thickgreen line being the output of the reference. The lower plot is the abso-lute difference of model output and reference output. The controller isadaptive for low adaptation gain, however it becomes static for highergains.390 5 10 15 20 25 30 35 40 45 50Time (sec)0123456789OutputLaw CReference = 1 = 1e3 = 1e40 5 10 15 20 25 30 35 40 45 50Time (s)10-2510-2010-1510-1010-5100|Error| = 1 = 1e3 = 1e4Figure 3.9: Simulated output for the plant G(s) controlled by the L1-AC forLaw C). The top figure is the output of the plant, with the thick greenline being the output of the reference. The lower plot is the absolutedifference of model output and reference output. The controller is non-adaptive for all gains.400 5 10 15 20 25 30 35 40 45 50Time (sec)024681012OutputLaw DReference = 1 = 1e3 = 1e40 5 10 15 20 25 30 35 40 45 50Time (s)10-1510-1010-5100105|Error| = 1 = 1e3 = 1e4Figure 3.10: Simulated output for the plant G(s) controlled by the L1-ACfor Law D). The top figure is the output of the plant, with the thickgreen line being the output of the reference. The lower plot is theabsolute difference of model output and reference output. The con-troller is nonadaptive for all gains.41Figure 3.11: Simulated output for the plant G(s) controlled by the L1-ACfor Law E). The top figure is the output of the plant, with the thickgreen line being the output of the reference. The lower plot is theabsolute difference of model output and reference output. The con-troller is nonadaptive for all gains.42Figure 3.12: Simulated output for the plant G(s) controlled by the L1-ACfor Law F). The top figure is the output of the plant, with the thickgreen line being the output of the reference. The lower plot is theabsolute difference of model output and reference output. The con-troller is nonadaptive for all dead-zone intervals.433.6 ConclusionThe L1-AC claims fast adaptation while maintaining the robustness. The fastadaptation is achieved through the use of high-gain feedback while robustness isachieved through the use of a low-pass filter that filters out any noise amplificationcaused by the high-gain feedback.In Section 3.4, it was shown that a PKPD model’s output approximates a refer-ence system with an implementable LTI controller as the gain increases. It was alsoshown that the system is adaptive for low gains, but non-adaptive for higher gains.In Section 3.5 the loss of adaptivity was mathematically shown. The feedback gainshows up in the loop transfer function, which when taken to infinity, causes theinversion of the predictor model. This inversion of the predictor’s model is a wellknown concept in classical control theory. This Chapter then showed a simple ex-ample for 6 different adaptation laws, some of which were adaptive for low gainsand some of which were non-adaptive for all gains. The output of all cases ap-proximates the implantable LTI controller of the reference system. This rules outthe use of those L1-AC schemes to address the problem of patient variability inclosed-loop control of anesthesia.44Chapter 4Real-Time Oscillation DetectionOscillation is a common problem in control-loop systems. There are three types ofoscillations: damped oscillation, undamped sustained oscillation, unstable oscilla-tion. Unstable oscillation will lead to increased deviation from the set-point andcan compromise safety and stability. While damped oscillation and the sustainedoscillation will not cause instability, they will lead to a lower quality control signal.Other than safety and stability concerns, oscillation will also result in a highercontrol action [58]. This translates into a higher drug dose to the patient, whichcan lead to post-surgical complications and unnecessary increase in cost.There are several factors that can cause oscillation. A few are:• Marginally stable control loops (due to aggressive control tuning or changesin process gain/phase/time delay).• External disturbances.• Dead-band, also known as hysteresis, of the controller valve (hardware/e-quipment issues).• Stiction in the control valve.• Hitting the upper and/or lower bound limit of the controller valve.Studying and analyzing the oscillation can be performed in three stages:451. Identifying sustained oscillations, and where possible, detecting the differentfrequency components of the oscillation(s).2. Detecting and quantifying the root-cause of the oscillation(s).3. Correcting the root-cause of the oscillation(s).It would be beneficiary for the detection to be performed in real-time. A flagcan be raised to alert the anesthesiologist. If the problem is due to poor-tuning,an auto-tuning method (see Chapter 5) can re-tune the controller and remove theoscillation. If the problem is due to a mechanical failure, the system can be takeninto manual mode to prevent the escalation of the problem.There have been numerous techniques to address the problem of detecting os-cillation in a signal that contains multiple oscillation frequencies. In more-or-lesschronological order, these methods include the integrated of absolute error [25], theauto-correlation function methods [32], the spectral peaks-based method [32], thewavelet-based method [44], the modified empirical mode decomposition method[55], the DCT-based method [67] and many more. Among all, the detection methodbased on the DCT proposed by Wang et. al. [67] is one of the most advanced meth-ods. It can detect multiple oscillations in off-line and real-time and it can determinethe frequency, magnitude and fitness (percent energy of the oscillation) of thesecomponents.In this Chapter, an algorithm for detecting oscillations is discussed. In Section4.2 an off-line method based on the DCT analysis is introduced. In Section 4.3, thismethodology is extended to real-time. In Section 4.4, examples are provided.4.1 ContributionThis Chapter discusses a new method for determining the multiple oscillations in acontrol signal, based on the method proposed by Wang et. al. [67]. First, the off-line algorithm is developed. The extension to real-time is discussed subsequently.The algorithm is able to also determine the dominant oscillation signal, character-ized by its frequency, magnitude and fitness.464.2 The DCT Off-line Oscillation DetectionTraditionally, the Fourier Transform and its discrete algorithm, DFT has been ap-plied to frequency-related problems. Therefore, it may seem natural to pick the DFTrather than the DCT, which is related to the complex portion of the DFT signal. Themain advantage of DCT is its strong ”energy compaction” property [4]: most of thesignal information is concentrated in a few coefficients of the low-frequency com-ponents of the transformed signal and it approaches the Karhunen-Loe´ve transform(which optimally decorrelates the frequency components, but is extremely slow tocompute).More importantly, the default DFT rectangular windows only provides good en-ergy compaction for frequencies that are whole fractions of the sampling frequencyFs, i.e. a DFT analysis of a signal with sampling frequency of 240Hz can effectivelydetect frequencies that are exact (or close to) multiples of 240/N (such as 120Hz,80Hz, or 60Hz). Applying a non-rectangular DFT window (or a moving window),will produce less broadband leakage, but will be lossy near the window’s edge.DCT addresses all the issues above, and is therefore used as the basis for oscillationdetection.Given a time series x(t), its associated frequencies will be distributed separatelyin the signal’s DCT counterpart y(k). That is, the different frequency componentsof the signal x(t) can be studied by observing different segments of y(k). In thefollowing Section, the oscillation detection algorithm is discussed. Noise is alsoconsidered in this discussion.4.2.1 The DCT DefinitionGiven a time series discrete sequence x(nT )∣∣Nn=1, with sampling period T , its DCTcounterpart is defined as:y(k) = ω(k)N∑n=1x(n)cos( pi2N(2n−1)(k−1)), k = 1,2,3, ...N, (4.1)47whereω(k) =1√Nk = 1,√2N2≤ k ≤ N.(4.2)Similarly, the Inverse Discrete Cosine Transform (IDCT) is defined as:xi(n) =N∑k=1ω(k)y(k)cos( pi2N(2n−1)(k−1)), n= 1,2,3, ...N. (4.3)The DCT signal y(k) has a convenient inherent property: given a signal of x(t)=sin(2piωt + φ), y(k) will always be of the form y(k) = ...,0, ...,#,0, ... where #stands for some non-zero value. In other words, for each frequency component inthe signal x(t), the counterpart y(k) will start with a zero directly followed by anon-zero value (call this 0− # pattern), then followed by some integer (could bezero or non-zero), and finally finished by a #−0 pattern. This property of the signalwill be used in the next section to extract the segment of y(k) that corresponds to aspecific frequency of the signal x(t).To visualize this DCT pattern, a signal with two frequency components isshown in Figure 4.1. The signal is shown on the top with its DCT shown at thebottom. The first 4 points follow the discussed pattern and contribute to one of thefrequency components. The next 5 points contribute to the second frequency com-ponent. The reconstructed signals of these two frequency components are shownin Figure 4.2. The reconstruction can only provide information on the frequencyof the signal, and not on the magnitude or the offset of it.4.2.2 The DCT AlgorithmThe segments in y(k) that are within each of the pattern 0−# and #−0 contributeto the different frequencies in the signal x(t). A signal contaminated by noise how-ever, may have all of the y(k)’s component non-zero. The case for white noise iscovered below and the case of colored noise is discussed in the subsequent Section.480 2 4 6 8 10 12 14 16 18 2002468Signal0 2 4 6 8 10 12 14 16 18 2002462 Frequency ComponentsRed signal with = 4 unitsGreen signal with = 12 units0 2 4 6 8 10 12 14 16 18 20-1001020DCT of SignalAn oscillation componentAn oscillation componentFigure 4.1: The signal blue has two components with frequencies 4 and 12units respectively (labeled red and green signals). The DCT of the bluesignal is shown on the bottom graph. The first 4 points on the DCTcaptures one of the frequency components, while the next 5 capturesthe other.Define σˆy as the estimated standard deviation of the signal y(k):σˆy =√1N−1N∑k=1(y(k)− 1NN∑k=1y(k))2. (4.4)The white noise can be filtered by suppressing the values of y(k) smaller than3σˆy and preserving the most significant components [67]:yh(k) =y(k) |y(k)| ≥HY,0 |y(k)| 3 (4.13)To understand the rationale behind this inequality, consider the signal Zi to be dueto equally randomly distributed arrivals, i.e. a random exponential distribution:fTi = λe−λ µT . (4.14)For an exponential distribution, the mean and the standard deviation are equal, i.e.µTi = σTi . A null hypothesis H0 : RT = 1 and the alternative hypothesis H1 : RT > 3is formed. If the condition 4.13 holds, the H0 is rejected and H1 is accepted; Ti isthen claimed to be regular and oscillatory.The sample mean and standard deviation 4.10 and 4.11 cannot be reliably cal-culated with less than 4 sample sets. It is suggested to use at least 10 sample sets[59]. A modified regulatory index is now defined that also considers the number of52sample sets.Define the population coefficient of variance as the ratio of population standarddeviation and the mean:Cv =σTiµTi, (4.15)and the sample coefficient of variance as:Cˆv =1RT=sTiT¯i. (4.16)Let α be a small positive integer such that (1−α)100% is the confidence intervalfor Cv: √L−1Cˆv√χ2L−1,1−α/2 3. (4.19)Equation 4.19 forms the first periodic test: if no Ti passes the period regulatorytest, then the signal x(t) is concluded to be non-oscillatory.The case for colored noise will now be discussed. If colored noise is present,the suppressed signal yh from 4.5 will have too many of its coefficient removedand no longer resembles the noise. As a result, oscillation detection may give falseresults. Instead, a low cut-off value LY is defined as:LY = σˆy, (4.20)and the signal y(k) is suppressed similar to yh 4.5, but with LY , to give yl . yl isthen segmented into its j-th DCT component, y j, similar to 4.3. A pair of yi and y j53that have the same maximum value are matched. The modified regulatory test isperformed on y j whose pair yi has passed the period regulatory test 4.19. If no y jpasses the test, the signal x(t) is concluded to be non-oscillatory.The definition of yl preserves the colored noise. The modified regulatory teston yl then filters out colored noise, while the test on yh filters out white noise. Thetwo cut-off values LY and HY are the two most important constants used in thisderivation. Li et al. performed a series of simulations with different colored andwhite noises to determine their respective value of LY = σˆy and HY = 3σˆy [38].4.2.2.2 Fitness TestTo measure the percentage energy of a component, or its fitness, the followingequation from [41] and [67] is used:F(x,xk) = 100(1− ‖ xk− x ‖‖ x ‖), (4.21)where ‖ x ‖ is the Euclidean norm. Of the (xi,x j) pair that have survived the modi-fied regulatory test, the component that gives the largest fitness, contains the mostenergy in the signal and is therefore the dominant frequency. Since x j containsmore coefficients than xi, it is used to determine the fitness of the dominant fre-quency of the signal x:Fd =maxjF(x,x j). (4.22)If this fitness is larger than a predefined threshold, F0, then the signal is concludedto be oscillatory. This is known as the fitness test:Fd > F0. (4.23)The dominant period is determined by comparing the modified regulatory indexRT,α of the ximax-x jmax pair that correspond to Fd:T¯d =T¯imax RTimax,α ≥ RTjmax,α ,T¯jmax otherwise.(4.24)54The signal x has a dominant period T¯d with a fitness of Fd and the correspondingcomponent ximax and x jmax.4.2.2.3 Magnitude TestThe two tests modified regulatory test and the fitness test are considered adequatefor most scenarios. However, sometimes the time series can pass both tests, butthe magnitude of the oscillation may not be periodic; therefore a further test on themagnitude of the oscillation must be performed.This test may seem unintuitive. An unstable oscillation (a sinusoidal that in-creases in magnitude) or a damped oscillation (a sinusoidal that decreases in mag-nitude) are both considered oscillatory, but have ”irregular” magnitude. The ir-regular pattern is a signal that has sudden drops and peaks (for instance, due tomiscommunication of the sensor) or one that has its magnitude follow an irregularpattern of high and low in no particular order. The magnitude test described belowpasses an unstable and damped oscillation and fails a true ”irregular pattern”. Thepaper Wang et al. [67] has an excellent example that illustrates this case.Define the magnitude series M(m) as:A(m) =max(x(t)∣∣1+lT¯d1+(l−1)T¯d)−min(x(t)∣∣1+lT¯d1+(l−1)T¯d),M(m) = A(m)/2,(4.25)where T¯d is the dominant period as determined by 4.24, l = 1,2,3, ...L. In otherwords, scan the time series x(t) in a window period of T¯d and subtract the maximumfrom the minimum of the sequence.Similar to 4.18, the magnitude index is defined as:RM,α =√χ2L−1,α/2√L−1M¯sM, (4.26)where M¯ and sM are the sampled mean and standard deviation of the signal Mrespectively. The magnitude regulatory test is:RM,α > 2.73, (4.27)55where the threshold is determined as follows: the signal M(m) is approximatelyhalf the size of T (l). The ratio of√χ2L−1,α/2√L−1 to√χ2L/2−1,α/2√L/2−1 is approximately 1.1 fordifferent values of α and L. The value 2.73 is 3/1.1 [67].4.2.3 Summary of the DCT AlgorithmThere are 3 tests: modified regulatory test 4.19, fitness test 4.23, and magnituderegulatory test 4.27. The tuning parameters are α and F0.The algorithm is summarized in the following 13 steps (from here on called theoscillation detection algorithm):Step 1. Remove the mean from the signal x(t) and compute the DCT y(k) from thedefinition 4.1.Step 2. Suppress the elements of y(k) that are smaller than the high cut-off valueHY 4.6 to generate yh as per 4.5.Step 3. Compute i-th DCT components of yh using 4.7 to get yi for i= 1,2,3, ..., I.Step 4. Generate the inverse DCT xi(t) for each yi(t) using 4.3.Step 5. Compute the period sequence Ti(n) for each xi(t) and perform the modifiedregulatory test. If no signal passes the test, then x(t) is concluded to be non-oscillatory.Step 6. Suppress the elements of y(k) that are smaller than the low cut-off valueLY 4.20 to generate yl .Step 7. Compute j-th DCT components of yl using 4.7 to get y j for j= 1,2,3, ...,J.Step 8. Select the y j that have the same maximum value as the yi whose xi passedthe regulator test from step 5.Step 9. Generate the inverse DCT x j(t) for each y j(t) that was selected form theprevious step.Step 10. Perform the same modified regulatory test as step 5. If none pass the test,then x(t) is concluded to not oscillatory.56Step 11. Calculate Fd as the dominant of the fitness of x jmax and perform the fitnesstest of 4.23. If the test fails, then x(t) is concluded to be non-oscillatory.Step 12. Determine the dominant period T¯d using 4.24.Step 13. Determine the magnitude sequence from T¯d and perform the magnitudetest of 4.27. If the test fails, then x(t) is concluded to be non-oscillatory.The oscillating signal is characterized by a dominant period T¯d, magnitude ofM¯ and fitness of Fd. This test has two tunable parameters: α and F0. Suggestedvalues are α = 2.7% and F0 = 25%.4.3 Extension to Real-TimeOscillation can lead to instability and excessive actuator action. Detecting oscilla-tion in real-time allows us to raise a flag and notify the anesthesiologist to quicklymodify the setting and to stabilize the patient’s DOH. In Chapter 5 an auto-tuningmethod is proposed that will re-tune the controller to remove the oscillation.The basis for extending the oscillation detection algorithm from Section 4.2.2is as follows: select an adaptive window range and perform the oscillation detec-tion algorithm. This window range is dependent on the predicted dominant periodand will change in real-time to adapt to the case.In Thornhill et al (1997) [58], it was suggested to use a window range of 50times of the presumed oscillation when applying the Ha¨gglund method. This valuewas then disputed by Wang et al (2013) [67] when applied to the DCT method sincetheir method specifically determines the oscillation period whereas the Thornhillmethod only detects large IAE errors and therefore requires a larger window range.The window range best be large enough to produce sufficient sample sets tocompute the sample mean and standard deviation of the period and magnitude se-quences T (l) andM(m) accurately. A window range of 10 times the presumed pe-riod will produce a period sequence of 20 sample sets, and a magnitude sequenceof 10 sample sets, allowing for an accurate measurement of the sample mean andstandard deviation of T¯ and M¯ [67]. The starting and ending positions of the time57window are:ne = t,ns = ne−10Tp,(4.28)where t is the current time.The presumed period, Tp should be adaptive to allow the system to identifyoscillation of any period. If a dominant period Td is found after applying the os-cillation detection algorithm, then it is set to be the presumed period Tp. If nooscillation is detected, the component pair xi and x j with the maximum Fx j is se-lected. There are 3 scenarios:Scenario 1. The pair xi and x j only contain one zero crossing. In this case, noperiod can be determined, and so the previous presumed period is kept.Scenario 2. The pair xi and x j contain exactly 2 zero crossing. In this case, themodified regulatory index cannot be calculated and so the method of 4.24cannot be applied. Since xi is contaminated less by noise, then T¯xi is selectedas the presumed period.Scenario 3. There are enough zero-crossing points to perform the modified regu-latory test. In this case, the method of 4.24 is used.In Wang et al (2013), it is suggested to allow the presumed period to changefreely. This can cause an issue: assume the system is initially oscillating with asmall period Ts. The oscillation then stops and at a later time an oscillation withmuch higher period Th > 10Ts is formed. The presumed frequency of the algorithmis now stuck at Ts and the system only scans a period of 10Ts and may never beable to capture this new oscillation.Instead, it is suggested to create multiple parallel instances of the algorithmto run simultaneously. Each instance has a predefined minimum and maximumallowed period, [Tlower,Tupper] and Tp is allowed to freely adapt in this period range.The range of each instance and the number of these instances will be the tuningparameter and is related to the problem at hand. This approach also allows usto ignore certain periods that are expected to exists in the system, for instance aknown background noise that might be present in the signal.584.3.1 Summary of the Real-Time DCT AlgorithmThe real-time algorithm can now be summarized in the following 5 steps. Thefollowing steps should be executed for the number of instances that have beenselected to run.Step 1. Specify an initial presumed oscillation to create the starting and endingpositions of the time window 4.28. Here Tp would be a priori knowledgeof the presumed oscillation. If this information is not available, then thebandwidth of the system can be used.Step 2. Wait until sufficient time has passed and there are sufficient data to performthe oscillation detection algorithm.Step 3. Perform the off-line oscillation detection algorithm on the segment ofx(ns) to x(ne).Step 4. Update the Tp according to the 3 scenarios 4.3 and the [Tlower,Tupper] limitsof the instance.Step 5. Repeat Step 2-4 for all the instances. You may need to wait for more dataif the new Tp is larger than the old one.4.4 Oscillation Detection ExamplesTwo examples are provided to highlight the oscillation detection algorithm. Thefirst example will be a simplified simulation example. The second example will bethe Depth of Hypnosis from a surgical case performed by iControl system.Example 4.4.1. Consider the following signal with an Signal to Noise Ratio (SNR)value of 10−1 shown in Figure 4.4:x(t) = sin2pi1.3t+2sin2pi3.4tThe off-line algorithm on the system determines two oscillation periods of 1.3 minand 3.397 min. Based on the fitness of the two signals, the component with periodof 3.397 min is chosen as the dominant table. The magnitude of this signal isdetermined to be 2.002. Table 4.1 summarizes the result.590 10 20 30 40 50 60Time-6-4-20246Figure 4.4: Signal from Example 4.4.1 contains two oscillations of peri-ods 1.3 min and 3.4 min. Dominant oscillation period is detected at3.397 min. The signal is shown in black and the dominant oscillationis shown in red.Table 4.1: High and Low components of Example 4.4.1Example 4.4.2. The following case is taken from one of the 61 surgical cases con-ducted at Royal Columbian Hospital in New Westminster using iControl system.Written consent was taken before the surgery from the patient. This particularpatient (case 6 from the database) underwent a Laparoscopic hemicolectomy.The oscillation starts at time 63min, and lasts until time 84min. The surgeryhad started at time 37min; there was a stimulation at time 63min, as recordedby the anesthesiologist. At time 80min, the patient moved. Immediately after,Rocuronium was administrated and the oscillation was damped out. The dominant60Figure 4.5: On-line oscillation detection shows a detected dominant signalof Tp = 3.63 min, M¯ = 5.96 and F = 89.31%.Table 4.2: Case example from Example 4.4.2. The magnitude regulatoryindex is 3.785.oscillation is Tp = 3.63 min, M¯ = 5.96 and F = 89.31%. The DOH is shown inFigure 4.5.4.5 ConclusionAlgorithms for detection of oscillation for both off-line (see Section 4.2.3) andon-line/real-time (see Section 4.3.1) were discussed. Unlike existing methods, thealgorithm can detect multiple oscillations, ignore specific oscillation frequencies,61and determine the dominant oscillation. The frequency, magnitude and fitness of allmeasured oscillations is also provided. The fitness of the oscillation can be used toreject small oscillations. The limitation of the algorithm discussed in this chapter,however, is that it requires a signal length of 10 times the presumed oscillationperiod. However, this is much less than Thornhill’s method that requires 50 timesthe presumed oscillation period. In the next Chapter, the frequency of the dominantwill be used to auto-tune a PID controller.62Chapter 5Re-Tuning of a PID ControllerThe closed loop feedback mechanism of PID controllers have found use in a varietyof systems, such as process control, motor, and vehicle control, to name a few[6]. The controller contains only three tunable parameters (the proportional k, theintegral ki, and the derivative kd), yet in many situations it can provide a robustsolutions with good performance [6]. Furthermore, this feedback system is wellunderstood and there are numerous implementations and theories to guarantee therobustness and performance.Our research group has worked extensively on PID control, see [10], [61], [62],[52], [54], [43] and [10]. PID controllers are not only of interest to our group, andother researchers have also investigated them, see [19], [31], [2], and [40].With the patient’s safety in mind, our robust controller is tuned to ensure areliable and safe drug administration, see all references above. However, therehave been cases where some oscillations have been observed in the clinical trials.It is therefore beneficial to have a system that would be able to detect oscilla-tion in real-time and automatically alert the anesthesiologist. Oscillation providesvaluable insight into the plant and the control loop [5]. It is possible to use this newinformation to re-tune the controller in real-time and remove oscillation.In this Section, a retuning mechanism that follows the guidelines for a robustPID controller design is discussed. The retuning mechanism is simulated with the44 PKPD models from Appendix A and compared with the original tuning of iCon-trol. In addition to removing the oscillation, the tuned system has met the following63objectives:• The gain margin, phase margin, and peak sensitivity should be the same orbetter than the original design. This translates to a gain margin of more than2 and phase margin of 30◦−60◦.• The system should have an overshoot of less than 10% for set point changeand disturbance rejection.• A rise time of 5-10 minutes is considered appropriate. However, even withthe current implementation of the closed-loop control, rise time performancecriteria is a secondary objective. The rise time of the original parametersshould be comparable to the tuned parameters.This Chapter presents a robust PID tuning method for the currently imple-mented iControl. The feedback controller structure with all the components ofiControl is shown in Section 5.2. The robustness and performance design require-ment is examined in Section 5.3. The tuning rules and optimization are discussedin Section 5.4. Finally, simulation results and comparisons are presented in Section5.7.5.1 ContributionA robust PID auto-tuning algorithm is presented in this Chapter. Using the fre-quency of a measured oscillation in real-time, the patient is identified. Oscillationis generally due to an aggreesive controller. The controller is then tuned to be lessaggressive and the sustained oscillation is removed. The tuning rule follows theguidelines of a robust controller. IE optimization is used as a performance crite-rion.5.2 Overview of Controller StructureConsider the 2-degree-of-freedom PID controller shown in Figure 5.1. r is the ref-erence DOH and y is the outputWAV . l and d are the input and output disturbancesrespectively, and n is the measurement noise. The surgical stimulus is represented64by d. The block diagram Gp is the patient, and Gc and G f f are the two LTI con-troller transfer functions that together describe the PID controller of the form:rGff Gp−Gc+ + +ul dyPID Controller+nFigure 5.1: A 2-degree-of-freedom PID controlleru(t) = k(br(t)− y(t))+ ki∫ t0(r(τ)− y(τ))dτ + kd(− dy(t)dt), (5.1)where k, ki, and kd are the proportional, integral, and derivative control respectively.b is a step response weighting parameter between 0 and 1. The derivative term doesnot act on the set-point since that can cause spikes during step changes [6].Comparing the PID definition 5.1 to the Figure 5.1, it is easy to realize that Gcacts on the signal y and G f f acts on the signal r:G f f (s) = bk+kis,Gc(s) = k+kis+ kds.(5.2)The goal of a PID controller is to track the reference signal r while rejecting anyload disturbance, measurement noise and process uncertainty. The relationshipbetween the four signal r, l, d and n to y are:y(s) =G f fGp1+GcGpr(s)+Gp1+GcGpl(s)+11+GcGpd(s)− Gp1+GcGpn(s) (5.3)The sensitivity function S(s) describes the transfer function from d(s) to y(s) andthe complimentary function T (s) describes the transfer function from r to y(s).65They provide valuable insight: one can design Gc to provide a reasonable distur-bance rejection and robustness to process uncertainty. G f f can then be used so thecontroller meets the performance design criteria. The specification of the robust-ness and performance will be discussed in Section 5.3.5.2.1 iControl Design StructureAn overview of iControl is provided below. More information on iControl can befound in [43] and [54]. The control structure is shown in Figure 5.2.FspFm−GcGp GNS++ + +l d nr rf u WAVcnsyfGffiControl StructureEvFigure 5.2: The iControl Structure5.2.1.1 Measurement and Reference FiltersThe DOH of the patient in the iControl structure, WAV, is measured by the Neu-roSense monitor. To attenuate the high frequency noise, WAV is passed through asecond order low-pass measurement filter with time constant Tm = 15s:Fm =11+ sTm+(sTm)2/2. (5.4)The reference signal is passed through a first-order low-pass set-point filter withtime constant Tsp = 25s to smooth out any step-like changes:Fsp =11+ sTsp. (5.5)The filtered signals r f and y f are given by:r f (s) =Gspr(s),y f (s) =Gmy(s).(5.6)665.2.1.2 Saturation and Integrator Anti-WindupThe infusion pump has a lower bound umin = 0ml/h and an upper bound umax =600ml/h. The controller’s output v is limited to the saturation values umin and umax.The saturation block diagram from Figure 5.2 is a non-linear dynamic. It canhowever be modeled by an ideal describing function N as defined in [7]:N =12(f1(au+δav)+ f1(au−δav)), (5.7)where the function f1 is given byf1(ρ) =1 ρ > 12pi(arcsinρ +ρ√1−ρ2) −1≤ ρ ≤ 1−1 ρ <−1,(5.8)where the constants au, av, and δ are given by:au = 0.5(umax−umin),av = 0.5(vmax− vmin),δ = u0− v0,u0 = 0.5(umax+umin),v0 = 0.5(vmax+ vmin),(5.9)where vmax and vmin denote the maximum/minimum of the controller’s action out-put v during an oscillation. This dynamic is only active when the controller issaturated.To prevent the windup that will result from the saturation, a classical trackinganti-windup scheme is implemented. The time constant Tt of the tracking anti-windup is 60 second.675.2.1.3 NeuroSense MonitorThe EEG signal E(t) can be translated to the DOH index WAV by the NeuroSensemonitor. The signal E(t) runs from 0 to 1 and the WAV spans from 0 to 100, with0 corresponding to an iso-electric EEG, and 100 corresponding to the fully awakestate. The dynamics of this monitor are described in a very simple LTI transferfunction [54]:GNS =1(8s+1)2. (5.10)5.2.1.4 Patient ModelThe patient Gp can be replaced by the PKPD model from Appendix A. The blockdiagram is shown again for reference in Figure A.2. The PK and PD models areLTI functions. The Hill function however, is a non-linear sigmoid function. Forcontrol design purposes, it needs to be linearized around the reference point. Thelinearized PKPD model of the patient is described in Appendix A and is given by:Gp(s) =kd · γ4 ·V1 ·EC50 ·(s+ k12) · (s+ k31)(s+ p1) · (s+ p2) · (s+ p3) · (s+ kd) · e−Tds. (5.11)5.2.1.5 Reference WeightingIn the first version of iControl, suppression of oscillations and rejection of distur-bance were prioritized over the performance of the system. The reference weight-ing b was set to zero. The reference signal only entered the control signal lawthrough the integrator action [54]. To improve performance and reduce inductiontime, the system was redesigned with a unity reference weighting (i.e. b= 1)[61].5.2.1.6 Current iControl ParametersThe iControl structure has 7 design parameters: k, ki, kd , b, Tt , Tsp and Tm. Thecurrent implementation of iControl has constant values for the anti-windup and the68filter constants:Tt = 60,Tsp = 25,Tm = 15.The reference weighting parameter b is set to 1. The PID parameters are based onthe age and weight of the patient. The Lean Body Mass (LBM) is defined in [27]as:LBM(w,h) =0.3281 ·w+0.33929 ·h−29.5336, if Male,0.29569 ·w+0.41813 ·h−43.2933, if Female,(5.12)where w is the weight in kilogram and h is the height in centimeter. The PIDparameters are defined as:c f = LBM ∗0.03,k = 0.081 · c f ,ki = 0.0055 · c f ,kd = 45 · c f .(5.13)5.3 Robustness and Performance DesignFollowing the iControl structure in Section 5.2, the loop transfer function is definedas:L(s) = Gc ·Gp ·GNS ·Fm ·N. (5.14)The controller parameters are matched to the patient model Gp. The actualpatient model may be different to the modeled PKPD Gp. It is important that thecontroller parameters not to be too sensitive to this process variability.The Nyquist plot of a loop function is shown in Figure 5.3 in blue. From theFigure, the amplitude (or gain) margin Am describes how much the gain of the loopfunction L(s) can change before the system become unstable. The phase marginφm describes how much the phase of the L(s) can change before instability is seen.69Figure 5.3: The Nyquist Plot with region of stabilityThe only uncertain function of the loop L(s) is the patient. Therefore, Am and φmquantify the upper bounds of how much the patient model can change before thesystem becomes unstable.The relationship from d to y is called the sensitivity function:S(s) =11+L(s). (5.15)It describes the amplification of the disturbance as a function of frequency. Themaximum, or the peak, of the S(s),Ms = max0≤ω≤∞| S(iω) |, (5.16)quantifies the worst-case amplification of the disturbance. This quantity is relatedto the gain margin Am and phase margin φm and is provided shortly.The Nyquist stability criterion defines the point where the function L(s) crossesthe negative x-axis at −1 as the instability point. To account for the uncertainty ofthe patient model, a circle of radius 1/Ms (the red circle) centered at -1 is intro-70duced. If the loop function is kept outside of this red circle, then the closed-loopsystem is guaranteed to have the specified gain margin Am and phase margin φm.There is a relationship between the gain and phase margin and the peak sensi-tivity:Am >MsMs−1 ,φm > 2arcsin12Ms.(5.17)Typical values of Ms are between 1.3-2. The typical values of Am are then between2-5, and the φm is between 30◦ and 60◦ [5].5.4 PID Auto-Tuning RulesGiven an oscillation of frequency and magnitude ωu andMu, where the subscript ustands for unstable, following Nyquist’s Stability Theorem, the loop function L(s)from Equation 5.14 is crossing the negative real-axis:L(iωu) =−1 (5.18)The idea for retuning is simple: the patient model is identified at this oscillationfrequency. Call this point Pu. The retuning of the PID controller is used to shape theloop function to a stable point at the same frequency. Call this point Ps. In Figure5.4, the dashed line represents the unstable loop function (prior to retuning). Thecontroller tuning then shapes the loop to the blue line, outside of the red circleregime.The tuning mechanism that follows is motivated by A˚stro¨m et al. in [6]. In theirbook, the authors assume the plant is known at a given point, and the controller istuned to shape the loop to the desired stable point. In the implementation here, theplant is identified at a point from the oscillation.71Figure 5.4: The model is identified at the unstable point Pu and the con-troller is tuned to take the loop function to the stable point Ps.5.4.1 Auto-Tuning for RobustnessLet the points Pu and Ps be given by their polar representations:Pu = 1eipi ,Ps = rseiφs .(5.19)The point Ps is outside of the red circle from Figure 5.3 of radius 1/Ms. Themagnitude and the phase of the stable point relate to the gain Am and phase φmmargins by the following relationships [22]:Am =1rs,φm = φs(5.20)where the gain and phase margins are defined in Section 5.3 and are related to Msvia Equation 5.17. The values used will be discussed in Section 5.5.Let the transfer functions of the loop function be described by the polar repre-72sentation of a complex system at the oscillation frequency ωu:Gc(iωu) = rceiφc ,Gp(iωu) = rpeiφp ,GNS(iωu) = rNSeiφNS ,Fm(iωu) = rmeiφm ,N(v) = rv.(5.21)The ideal describing function N has zero phase [7]. It is assumed that the controlaction v will also oscillate and is saturated (i.e. the describing function’s dynamicneeds to be considered).The magnitudes and phases of Gc (from 5.2), GNS (from 5.10), Fm (from 5.4),and N(v) (from 5.7) can be computed. The identified rp and φp is given by:rp =1rc · rNS · rm · rv ,φp =−(φc+φNS+φm).(5.22)A new PID controller Gct is determined such that it moves the loop function to thepoint Ps at the same frequency ωu:Ltuned(iωu) = Ps = rseiφs = Gct (iωu) ·Gp(iωu) ·GNS(iωu) ·Fm(iωu). (5.23)There are two observations:1. The tuned PID controller action will no longer oscillate. The non-linear sat-uration dynamics N does not need to be considered.2. Since the new loop is still computed at the frequency ωu, the functionsGp(iωu),GNS(iωu), and Fm(iωu) have the samemagnitude and phase as 5.21.Substituting the values of Gp(iωu),GNS(iωu), and Fm(iωu) from 5.21 into Equation5.23 allows us to solve for the new controller’s magnitude and phase:rct =rsrp · rNS · rm · rv =rsrc · rv ,φct = φs−(φp+φc+φNS+φm)= φs−φc.(5.24)73The PID controller has three tuning parameters, k, ki and kd . The tuning ruleabove provides two equations; to get a unique solution, a third equation is needed.In A˚stro¨m et al, it is suggested to define a ratio between ki and kd and set the valueof the ratio by trial and error [6]. In the next section, a condition on the performanceof the system is imposed on the tuning rule instead.5.4.2 Auto-Tuning for PerformanceIn the preceding section, a set of robust rules were defined for auto-tuning thecontroller. By defining an appropriate gain and phase margin, the tuned controllerwill be robust and will have good output disturbance rejection.A common performance criterion is the ability to reject the load disturbance.The Integrated Error (IE) can be used as the metric to measure this performance.A small value of IE indicates a fast load disturbance rejection and a small steadystate error. It has also been shown that a step-like reference applied at the process’sinput is directly related to the PID’s integrator gain, ki [6]:IE =∫ ∞0(r(t)− y(t))dt = 1ki. (5.25)The third tuning rule is to maximize the integrator gain ki so as to minimize the IEand obtain a good load disturbance rejection.5.4.3 Bumpless Parameter ChangeUpon system retuning, parameter change will naturally change the controller’s out-put. This would cause a bump as the system’s states prior and after the parameterschange may not coincide. Care must be taken if a bumpless parameter change isrequired.To ensure a bumpless controller action, it is shown in [6] that is it sufficient toensure the controller output due to the proportional and the integral component (la-beled P+ I) is invariant to the parameter change. This can be achieved by requiringthe state of the integrator to change as:Inew = Iold + kold(r− y)− knew(r− y), (5.26)74where kold and knew are the old and new proportional gain respectively [6].5.5 Auto-Tuning ImplementationThe auto-tuning problem is a minimization problem:Minimize Equation 5.25 given the constraints 5.24 for the given gainAm and phase φm margins and/or the peak sensitivity Ms.To allow flexibility for this optimization problem, the gain and phase require-ments are provided as a range: an Ms range of 1.3 to 2 is used to guarantee a gainmargin of 2.11 to 4.33 and a phase margin of 31◦ to 45◦.The optimization also needs to be bounded. Otherwise, the tuned parameterscan become unbounded (ki is maximized to minimize IE). Root causes of oscil-lation were discussed in Chapter 4. It is assumed that the PID is initially properlytuned. Further, it is assumed there is no pump stagnation. The cause of oscillationis assumed to be due to unmatched model uncertainty: the PID controller is simplytoo aggressive for the patient.Let k, ki and kd be the current PID parameters, as defined by 5.13. Let the k′, k′iand k′d be the tuned parameters. The tuned parameters are expected to be smallerthan the original parameters. This is certainly true for the proportional gain k,however the integral and the derivative gain may need to increase slightly to satisfythe constraints 5.24. The upper and lower bounds of the new PID parameters aregiven as:(k′lower ,k′upper) = (0.75k,k),(k′ilower ,k′iupper) = (0.75ki,1.1ki),(k′dlower ,k′dupper) = (0.75kd ,1.1kd).(5.27)The lower bounds are set to prevent a slow response to stimulation and rise time.The exact values for these parameters were determined via simulation for the 44PKPD models of [10].The optimization is solved using MATLAB’s fmincon interior-point algorithmusing MaxFunEvals= 1e10, MaxIter = 1e3,755.6 Summary of Auto-Tuning AlgorithmThe steps below summarizes the auto-tuning algorithm discussed in this chapter.Prior to running the algorithm, two stability points based on Ms of 1.3 and 2 aredefined. These will be the upper and lower acceptable stability points.1. Run the real-time oscillation detection algorithm discussed in Section 4.3.1to determine the dominant oscillation frequency.2. The magnitude and phase of the controller, filters and dynamics of the Neu-roSense monitor can now be computed using Equations 5.2, 5.10, 5.4 and5.7.3. Identify the patient model at the oscillation frequency using Equation 5.22.4. Solve for the new PID parameters by minimizing the function 5.25 subject tothe conditions of 5.24. Set the upper and lower bounds of the PID parametersas defined in 5.27.5.7 Simulation Examples and ResultsTo assess the tuning robustness and performance, the system is simulated. A stepchange to 50 WAV is applied at time zero. Measurement noise modeled by Solteszet al. from [52] is applied at time 75 minutes until time 85 minutes. Disturbance,also modeled by Soltesz et al. is applied from time 95 minutes to time 135 minutes.The disturbance models a surgical stimulus. The system initially starts with thetuned parameters. After induction of anesthesia is complete, the system is tunedto be unstable. The oscillation detection algorithm from Chapter 4 detects theoscillation and the algorithm from this Section is used to re-tune the system.Figures 5.5 to 5.8 shows a simulation examples of one patient model from eachof the four groups in the PKPD models of Bibian [10]. The black line shows theoutput of the controller with the original PID tuning and the blue line shows theoutput of the controller with auto-tuning. The blue line starts with the same PIDparameters as the black line. After induction is complete, the tuning is turned to beunstable, causing instability in the output. The controller is auto-tuned and the blue76Figure 5.5: Group 1, Case 10: The original tuning has Am = 8.72 and φm =60.29. The re-tuned system has Am = 8.62 and φm = 53.81.Figure 5.6: Group 2, Case 17: The original tuning has Am = 6.27 and φm =61.20. The re-tuned system has Am = 6.13 and φm = 55.64.77Figure 5.7: Group 3, Case 33: The original tuning has Am = 7.57 and φm =65.21. The re-tuned system has Am = 7.43 and φm = 58.96.Figure 5.8: Group 4, Case 38: The original tuning has Am = 6.89 and φm =61.14. The re-tuned system has Am = 6.86 and φm = 55.71.78signal approaches the black signal. The behavior of all four examples is similar andit highlights the effectiveness of the auto-tuning algorithm to effectively tune thecontroller and achieve the same disturbance and noise cancellation as well as thesame performance as the original tuning.Tables D.1 to D.3 in the Appendix show the robustness, output disturbance re-jection, and step-change response for all 44 models. For robustness (Table D.1),the amplitude and phase margins are compared. For output disturbance rejection(Table D.2), the maximum overshoot assesses controller’s initial response to a 20%WAV disturbance. The settling time Ts measures the time it takes for the WAV toreach within 10% of the set-point following the overshoot. Finally, the IAE mea-sures how well the system rejects the disturbance. For setpoint-change response(Table D.3), rise time Tr measures the time it takes for the WAV to reach 80% of theset-point change. Settling time Ts is similar to the output disturbance rejection.The response of the auto-tuned cases is slightly slower than the original tuning,but otherwise follows them very closely. The slower response is to counter theaggressive tuning that was imposed on the original tuning. Table D.4 shows thetuning parameters. The proportional gain k in all cases is lower than the originaltuning.In the tuning algorithm used, no prior knowledge of the patient model was used.The only assumption is that the initial PID parameters are properly tuned using theprior knowledge of the patient. The measured oscillation is due to unmatchedpatient uncertainty, though the unmatched parameter is unknown; there is no newinformation about the patient model.The median (min, max) of the amplitude margin of the original tuning is 8.55(4.72, 17.36). The median (min, max) of the phase margin is 61.6◦ (50.3◦, 67.2◦).The median (min, max) of the amplitude margin of the re-tuned system is 8.46(4.75, 12.21). The phase margin is 56.4◦ (48.8◦, 62.1◦). The re-tuned systemis well within the minimum robustness requirement and agrees with the originaltuning [5].The retuning does not come at a great cost to the output disturbance rejection.The median (min, max) of the IAE of the original tuning is 202.7 (137.5, 304.9).The median (min, max) of the re-tuned system is 209.7 (143, 293.0). The distur-bance rejection has slightly increased.79The median (min, max) of the overshoot of the original tuning is 5.26% (3.0%,6.57%). For the re-tuned system, it is 5.99% (3.01%, 7.63%). All overshoots arestill under 10% and is in-line with the original design criteria.5.8 ConclusionIn this Chapter, oscillation is measured in real time. The frequency of the dominantcomponent is used to automatically re-tune the controller to remove the oscillation.The tuning rule is inspired by A˚stro¨m et al. [6] is based on a robust design. Theplant is identified at the oscillating frequency. A new controller is tuned to shapethe loop function to a stable region with a gain margin of more than 2 and a phasemargin of 30◦− 60◦. Disturbance rejection is guaranteed by minimizing the IE.Percent Overshoot (OS) is kept under 10%. Using the PKPD models of Bibian, thetuned controller is shown to be comparable to the original iControl.The tuned system guarantees stability at the measured frequency only. Thesystem may still be unstable and oscillate at another frequency. The tuning algo-rithm should therefore keep the record of all recorded oscillations. On each succes-sive retuning, the optimization constraints should include the list of all previouslyrecorded oscillations.80Chapter 6Conclusions6.1 Summary and ContributionsThe objectives of this thesis were to 1) assess the applicability of the novel L1-ACas applied to closed-loop control in anesthesia using the WAV index as the controlsignal; 2) design an oscillation detection algorithm that can detect multi-periodoscillation in real time; and 3) develop a tuning algorithm that can re-tune the PIDcontroller used in iControl to remove the detected oscillation.Following the requirements for a fast adaptive algorithm with guaranteed ro-bustness, the L1 Adaptive Control (L1-AC) was reviewed. This controller claimsfast adaptation while maintaining robustness. The fast adaptation is achieved byusing a high gain feedback and robustness is achieved by filtering out the high-frequency components of the feedback law using a low-pass filter. It was shownthat L1-AC loses its adaptivity as the gain of the adaptation law increases. Further,the resulting limiting controller can be achieved using an implementable LTI systemwhose dynamics depend only on the reference model, and not on the patient’s un-known parameters. Furthermore, the loss of adaptivity was mathematically shownto be a consequence of the well-known inversion of nonlinearity due to high-gainfeedback.The majority of oscillation detection algorithms currently in practice cannotguarantee the detection of oscillation if multiple oscillation frequency exists. More-over, most of these algorithms are boolean and can only determine whether or not81an oscillation exists. In a complex systems, multiple oscillations can develop si-multaneously. The algorithm of Chapter 4 can 1) determine multiple oscillations;and 2) determine the frequency, magnitude and fitness of each measured oscilla-tion. The fitness of the oscillation can be used to reject insignificant oscillations.The frequency and magnitude of the oscillations can identify the plant. The algo-rithm requires signal for at least 10 times the presumed oscillation period, whichmay limit its feasibility for short surgeries. For instance, a presumed oscillationperiod of 3 minutes can only be accurately measured after 30 minutes. While thisis not clinically relevant, it is still better than what is required by other methods -some require 50 times the length of the presumed oscillation period.One of the biggest challenges of closed-loop control of anesthesia is the inter-patient drug-response variability. The uncertainty in the PKPD model of patientsis a challenge for designing a controller than can be both robust and perform wellby rejecting surgical stimuli as well as following step-responses. This may leadto an aggressive controller that can cause oscillations. The iControl system wasmodified to automatically re-tune itself when oscillation was detected subject tothe following design objectives:• The gain margin, phase margin, and peak sensitivity are the same or betterthan the original design. This would be a gain margin of more than 2 andphase margin of 30◦−60◦.• The system must have an overshoot of less than 10% for set point changeand disturbance rejection.• A rise time of 5-10 minutes is considered appropriate. However, even withthe current implementations of a closed-loop control, rise time performancecriteria is a secondary objective. The rise time of the original parameterswith the tuned parameters should be comparable.The robustness and performance of the re-tuned PID controller was comparedwith the original iControl tuning. The PID re-tuning was applied to 44 PKPD modelsby Bibian. In all cases, the controller was initially properly tuned according to thelatest iControl version. After induction was completed, the system was re-tunedto be unstable and cause oscillation. The dominant oscillation was automatically82measured and the system was then re-tuned. In all 44 cases, it was shown that there-tuned controller had similar robustness and performance behavior as the originaliControl tuning.The current standard of performance evaluation of closed-loop control is theset of Varvel measures. These measures are not adequate for the current struc-tures of closed-loop control since they cannot be used design criteria. A set ofproposed measures by Soltesz et al. were shown to correlate with the Varvel mea-sures. Unlike Varvel, however, these proposed measures are accepted within thecontrol community and are used as performance criteria. The Integrated Error (IE),Percent Overshoot (OS), and Induction Phase Duration (ID) were used as designobjectives for re-tuning the PID controller.6.2 Future WorkThe possible future directions based on the work of the thesis are outlined below.The works are separated as ”Imminent” and ”Distant” directions.6.2.1 Imminent Future DirectionA limitation of the oscillation detection is the required length data of 10 timesthe presumed oscillation period. However, the system can moderately predict andwarn the anesthesiologist with less data. For example, the system could alert theanesthesiologist by monitoring only 3 or 4 times the presumed oscillation period.This could be displayed as a ”probable” oscillation. As more data is gathered, theconfidence on the measured oscillation can increase to ”certain”. Visually, the DOHcan be in green when no oscillation is present. It can turn to yellow when there isa probable oscillation, and finally red when oscillation is detected.The re-tuning algorithm for the PID controller can only guarantee robustnessfor the measured oscillation frequency. Another oscillation may occur at a differentfrequency at a later time. The tuning algorithm should therefore keep a record of allmeasured oscillations and on every successive oscillation, include all the recordedoscillations as constraints on the re-tuning optimization. Such a system can bebeneficial for use in the ICU, where a patient may be sedated for a few days.The oscillation detection algorithm can form the basis of a measurement that83can give a score to how oscillatory a signal is. This can complete the proposedSoltesz alternatives to Varvel metrics. The proposed measures currently does notquantify oscillation and Varvel’s Wobble metric does not have a substitute.The proposed alternative measures to the Varvel metrics requires more studyand verification. This is needed to create a set of measures that is acceptable byboth the clinicians and control engineers. Without these, no two closed-loop con-trollers can be compared. These measures can also facilitate communication be-tween the clinicians and the control engineers and can act as excellent diagnosticmetrics.6.2.2 Distant Future DirectionThe L1-AC in its current form cannot guarantee the fast adaptation it claims, norcan it guarantee adaptivity at a high adaptation speed. However, adaptive con-trollers may be the only feasible solution to a completely autonomous closed-loopcontroller that is truly robust and performs well despite surgical stimuli. Model-Predictive-Controller (MPC) is a promising adaptive controller that has been studiedby other researchers. Our own research team is also working on MPC.Finally, the phrase ”closed-loop control of anesthesia” has been loosely used todescribe systems that only control the Depth of Hypnosis (DOH). A true control ofDepth of Anesthesia (DOA) requires monitoring several physiological signals in-cluding the EEG, blood pressure, heart rate, respiratory rate, heart rate, etc. Thereare also a variety of hypnotic, opioid, and neuromuscular drugs that are adminis-trated to the patient to achieve a full state of anesthesia. From the controller’s pointof view, this corresponds to a MIMO system. The ideal controller should measureDOA (both hypnotic and analgesic state) as well as level of paralysis, and automati-cally control the infusion of all anesthetic drugs. This MIMO controller is currentlynot available, but should be the holy grail of the closed-loop control of anesthesia.84Bibliography[1] A. Absalom and G. Kenny. Closed-loop control of propofol anaesthesiausing bispectral index: performance assessment in patients receivingcomputer-controlled propofol and manually controlled remifentanilinfusions for minor surgery. British Journal of Anaesthesia, 90(6):737–741,2003. → pages 12, 13[2] A. Absalom, N. Sutcliffe, and G. Kenny. Closed-loop control of anesthesiausing bispectral index: Performance assessment in patients undergoingmajor othopedic surgery under combined general and regional anesthesia.Anesthesiology, 96:67–73, 2002. → pages 4, 7, 12, 13, 63[3] G. Agrawal, S. Bibian, and T. Zikov. Recommended clinical range forwavCNS index during general anesthesia. In Annual Meeting of the AmericanSociety Anesthesiologists, 2010. → pages 100[4] N. Ahmed, T. Natarajan, and K. Rao. Discrete consine transform. IEEETransactions on Computers, 1:90–93, 1974. → pages 47[5] K. J. A˚stro¨m and T. Ha¨gglund. Advanced PID Control. ISA - TheInstrumentation, Systems, and Auomation Society, 2006. → pages 14, 63,71, 79[6] K. J. A˚stro¨m and R. M. Murray. Feedback Systems: An Introduction forScientists and Engineers. Princeton University Press, 2008. → pages iii, 63,65, 71, 74, 75, 80[7] D. Atherton. Nonlinear Control Engineering. Van Nostrand ReinholdLondon, 1982. → pages 67, 73[8] M. Avidan, L. Zhang, B. Burnside, K. Finkel, A. Searleman, J. Selvidge,L. Saager, M. Turner, S. Rao, M. Bottros, C. Hantler, E. Jacobsohn, andA. Evers. Anesthesia awareness and the bispectral index. New EnglandJournal of Medicine, 2008. → pages 8, 985[9] S. Babji, U. Nallasivam, and R. Rengaswamy. Root cause analysis of linearclosed-loop oscillatory chemical process systems. Industrial & EngineeringChemistry Research, 51:13712–13731, 2012. → pages 13[10] S. Bibian. Automation in Clinical Anesthesia. PhD thesis, University ofBritish Columbia, 2006. → pages viii, xii, 3, 5, 9, 10, 14, 63, 75, 76, 92, 94,96, 97, 98[11] S. Bibian, T. Zikov, G. A. Dumont, C. Ries, H. Puil, H. Ahamdi,M. Huzmezan, and B. A. MacLeod. Estimation of anesthetic depth usingwavelet analysis of electroencephalography. In EMBC InternationalConference, 2001. → pages 2, 9[12] S. Bibian, T. Zikov, G. A. Dumont, C. Ries, H. Puil, H. Ahamdi,M. Huzmezan, and B. A. MacLeod. Method and apparatus for the estimationof anesthetic depth using wavelet analysis of the electroencephalogram,2004. URL http://www.google.ca/patents/US20040010203. → pages 2[13] J. D. Boskovic and R. K. Mehra. Performance analysis of a simple L1adaptive control. In American Control Conference (ACC), 2013. → pages23, 33[14] C. Cao and Hova. Stability margins of L1 adaptive controller. IEEETransactions on Automatic Control, 55(2):480–487, 2010. → pages 35[15] C. Cao and N. Hovakimyan. L1 adaptive controller for systems withunknown time-varying parameters and disturbances in the presence ofnon-zero trajectory initialization error. International Journal of Control, 81(2):586–591, 2008. → pages 35[16] CareerCast.com. Jobs rated 2013: Ranking 200 jobs from best to worst.http://www.careercast.com/jobs-rated/best-worst-jobs-2013, 9 2013.Accessed: 2014-09-04. → pages 7[17] CDC. National hospital discharge survey: 2010 table, procedures by selectedpatient characteristics - number by procedure category and age. Technicalreport, Centers for Disease Control and Prevention, 2010. → pages 1[18] D. J. Cullen, E. Eger, W. C. Stevens, N. T. Smith, T. H. Cromwell, B. F.Cullen, G. A. Gregory, S. H. Bahlman, W. M. Dolan, R. K. Stoelting, andH. Fourcade. Clinical signs of anesthesia. Anesthesiology, 36:21–37, 1972.→ pages 186[19] K. Ejaz and J. S. Yang. Controlling depth of anesthesia using pid tuning: Acomparative model-based study. In International Conference on ControlApplication, 2004. → pages 63[20] B. D. Franklin, K. OGrady, P. Donyai, A. Jacklin, and N. Barber. The impactof a closed-loop electronic prescribing and administration system onprescribing errors, administration errors and staff time: a before-and-afterstudy. International Journal of Healthcare Improvement, 2007. → pages 7[21] C. W. Frei, A. Gentilini, M. Derighetti, A. H. Glattfelder, M. Morari,T. Schnider, and A. M. Zbinden. Automation in anesthesia. In AmericanControl Conference, 1999. → pages 11[22] M. Friman. Automatic re-tuning of pi controllers in oscillating controlloops. Technical report, Abo Akademi University, 1997. → pages 72[23] G. Goodwin, S. Graebe, and M. Salgado. Control System Design. PrenticeHall, NJ, USA, 2001. → pages 33, 108[24] M. Gruenewald, C. Ilies, J. Herz, T. Schoenherr, A. Fudickar, J. Hcker, andB. Bein. Influence of nociceptive stimulation on analgesia nociception index(ani) during propofolremifentanil anaesthesia. British Journal ofAnaesthesia, 110:1024–1030, 2013. → pages 2[25] T. Ha¨gglund. A control-loop performance monitor. Control EngineeringPractice, 3:1543–1551, 1995. → pages 13, 46[26] N. Hovakimyan and C. Cao. L1 Adaptive Control Theory: GuaranteedRobuRobust with Fast Adaptation. Society for Industrial & AppliedMathematics, 2010. → pages ix, xii, 3, 14, 22, 24, 26, 28, 104, 106[27] R. Hume. Prediction of lean body mass from height and weight. Journal ofClinical Pathology, 19:389–391, 1966. → pages 69[28] P. A. Ioannou and J. Sun. Robust Adaptive Control. Dover Publications, Inc,2012. → pages 14, 22[29] S. Jafari, P. Ioannou, and L. Rudd. What is L1 adaptive control. In AIAAGuidance, Navigation, and Control (GNC) Conference, 2013. → pages 23,33[30] M. Jeanne, C. Clement, J. de Jonckheere, R. Logier, and B. Tavernier.Variations of the analgesia nociception index during general anaesthesia forlaparoscopic abdominal surgery. Journal of Clinical Monitoring andComputing, 26:289–294, 2012. → pages 287[31] G. Kenny and H. Mantzaridis. Closed-loop control of propofol anaesthesia.British Journal of Anaesthesia, 83:223–228, 1999. → pages 12, 63[32] S. Kerra, M. Jelali, M. N. Karim, and A. Horch. Detection and diagnosis ofstiction in control loops, chapter Detection of oscillating control loops, pageChapter. 14. Springer, 2009. → pages 46[33] E. Kharisov and N. Hovakimyan. Generalization of L1 adaptive controlarchitecture for switching estimation laws. In American Control Conference,2012. → pages 37, 38[34] E. Lavretsky and T. E. Gibson. Projection operator in adaptive systems.Technical report, Cornell University Library, 2011. → pages 25[35] E. Lavretsky and K. A. Wise. Robustness and Adaptive Control. Springer,2013. → pages 22[36] T. Ledowski, W. Tiong, C. Lee, B. Wong, T. Fiori, and N. Parker. Analgesianociception index: evaluation as a new parameter for acute postoperativepain. British Journal of Anaesthesia, 111:627 to 629, February 2013. →pages 2[37] T. L. Lemke, D. Williams, R. V. F., and S. W. Zite. Foye’s Principles ofMedicinal Chemistry. Lippincott Williams & Wilkins, 2012. → pages 2[38] X. Li, J. Wang, B. Huang, and S. Lu. The dct-based oscillation detectionmethod for a single time series. Journal of Process Control, 20:609–617,2010. → pages 54[39] N. Liu, T. Chazot, T. Chazot, A. Gentr, A. Landais, A. Restoux, K. McGee,P. Laloe, B. Trillat, L. Barvis, and M. Fischler. Titration of propofol foranesthetic induction and maintenance guided by the bispectral index:closed-loop versus manual control: A prospective, randomized, multicentralstudy. Anesthesiology, 2006. → pages 7, 11, 12, 97, 100, 101, 102[40] N. Liu, T. Chazot, S. Hamada, A. Landais, N. Boichut, C. Dussaussoy,B. Trillat, L. Beydon, E. Samain, D. Sessler, and M. Fischler. Closed-loopcoadministration of propofol and remifentanil guided by bispectral index: arandomized multicenter study. Anesthesia & Analgesia, 112:546–557, 2011.→ pages 63[41] L. Ljung. System identification toolbox 7: Getting started guide.MathWorks, 2010. → pages 5488[42] S. Locher, K. Stadler, T. Boehlen, T. Bouillon, D. Leibundgut,P. Schumacher, R. Wymann, and A. M. Zbinden. A new closed-loop controlsystem for isoflurane using bispectral index outperforms manual control.Anesthesiology, 101(3):591–602, 2004. → pages 12[43] A. Martinez. Robust control: Pid vs fractional control design, a case study.Master’s thesis, University of British Columbia, 2005. → pages 9, 63, 66,92, 97[44] T. Matsuo, H. Sasaoka, and Y. Yamashita. Detection and diagnosis ofoscillations in process plants. In Process of Seventh Int. Conf.Knowledge-Based Intelligent Information and Engineering Systems, pages1258–1264, UK, 2003. → pages 46[45] R. Ortega. On the effect of input filtering and fast adaptation in modelreference adaptive control. In American Control Conference, 2013. → pages23[46] M. Ralph, L. Beck, and M. Bloom. L1-adaptive mmethod for control ofpatient respones to anesthesia. In American Control Conference, 2011. →pages 13[47] T. Schnider, C. Minto, P. Gambus, C. Andersen, D. Goodale, S. Shafer, andE. Youngs. The influence of method of administration and covariates on thepharmacokinetics of propofol in adult volunteers. Anesthesiology, 88:1170–1182, 1998. → pages 9[48] J. Schu¨ttler and H. Imsen. Population pharmacokinetics of propofol: Amulticenter study. Anesthesiology, 92:727–738, 2000. → pages viii, 9, 94,95[49] L. Sheiner, S. D.R, S. Vozeh, R. Miller, and J. Ham. Simultaneouslymodeling of pharmacokinetics and pharmacodynamics: application tod-tubocurarine. Clinical Pharmacology and Therapeutics, 25:358–371,1979. → pages 10, 96[50] J. Sigl and N. G. Chamoun. An introduction to bispectral analysis for theelectroencephalogram. Journal of Clinical Monitoring, 10(6):392–404,1994. → pages 2, 8[51] J. Smith. The best- and worst-paying jobs in america.http://www.forbes.com/sites/jacquelynsmith/2013/05/13/the-best-and-worst-paying-jobs-in-america-2/,5 2013. Accessed: 2014-09-04. → pages 789[52] K. Soltesz, K. van Heusden, G. A. Dumont, T. Ha¨gglund, C. L. Peterson,N. West, and M. Ansermino. Closed-loop anesthesia in children using a pidcontroller: A pilot study. In IFAC, 2012. → pages 2, 3, 7, 14, 63, 76[53] K. Soltesz, G. A. Dumont, and M. Ansermino. Assessing controlperformance in closed-loop anesthesia. InMediterranean Conference onControl & Automation, 2013. → pages 4, 15, 101[54] K. Soltesz, J. Hahn, T. Ha¨gglund, G. A. Dumont, and M. Ansermino.Individualized closed-loop control of propofol anesthesia: A preliminarystudy. Biomedical Signal Processing and Control, 8:500–508, 2013. →pages 3, 63, 66, 68[55] R. Srinivasana, R. Rengaswamy, and R. Miller. A modified empirical modedecomposition(emd) process for oscillation characterization in control loops.Control Engineering Practice, 15:1135–1148, 2007. → pages 46[56] M. Struys, T. Smet, L. F. M. Versichelen, S. Van de Velde, R. Van denBroecke, and E. Mortier. Comparison of closed-loop controlledadministration of propofol using bispectral index as controlled variableversus ”standard practice” controlled administration. Anesthesiology, 95:6–17, 2001. → pages 12[57] K. Talebian, K. Soltesz, G. A. Dumont, and M. Ansermino. Clinicalassessment of control performance in closed-loop anesthesia. In AmericanSociety of Anesthesiologists, 2013. → pages 13, 15[58] N. Thornhill and T. Ha¨gglund. Detection and diagnosis of oscillation incontrol loops. Control Engineering Practice, 5:1343–1354, 1997. → pages45, 57[59] N. Thornhill, B. Huang, and H. Zhang. Detection of multiple oscillations incontrol loops. Journal of Process Control, 13:91–100, 2003. → pages 52[60] K. van Heusden and G. A. Dumont. Analysis of linear L1 output feedbackcontrol: Equivalent lti controllers. In 16th IFAC Symposium on SystemIdentification, pages 1472–1477, Brussels, Belgium, 2012. → pages 33[61] K. van Heusden, G. A. Dumont, K. Soltesz, C. L. Peterson, A. Umedaly,N. West, and M. Ansermino. Design and clinical evaluation of robust pidcontrol of propofol anesthesia in children. IEEE Transactions on ControlSystems Technology, 22:491–499, 2014. → pages 15, 63, 6890[62] K. van Heusden, N. West, A. Umedaly, M. Ansermino, R. N. Merchant, andG. A. Dumont. Safety, constraints and anti-windup in closed-loopanesthesia. In IFAC, 2014. → pages 63[63] K. van Heusden, K. Talebian, and G. A. Dumont. Analysis of l1 adaptivestate feedback control. why does it approximate an implementable lticontroller? European Journal of Control, 23:1– 7, May 2015. → pages 23,33[64] J. Vannes, E. Kharisov, and Hov. L1 adaptive control with proportionaladaptation. In American Control Conference, 2012. → pages 23, 35, 36, 37[65] J. R. Varvel, D. L. Donoho, and S. Shafer. Measuring the predictiveperformance of computer-controlled infusion pumps. Journal ofPharmacokinetics and Biopharmaceutics, 20:63–94, 1998. → pages 3, 99[66] F. C. Vasella, P. Frascarolo, D. R. Spahn, and L. Magnusson. Antagonism ofneuromuscular blockade but not muscle relaxation affects depth ofanaesthesia. British Journal of Anaesthesia, 94:742–747, 2005. → pages 2[67] J. Wang, B. Huang, and S. Ku. Improved dct-based method for onlinedetection of oscillations in univariate time series. Control EngineeringPractice, 21:622–630, 2013. → pages iii, 13, 46, 49, 50, 54, 55, 56, 57[68] N. West, G. A. Dumont, K. van Heusden, C. L. Peterson, S. Khosravi,K. Soltesz, A. Umedaly, E. Reimer, and M. Ansermino. Robust closed-loopcontrol of induction and maintenance of propofol anesthesia in children.Pediatric Anesthesia, 23:712–719, 2013. → pages 3[69] T. Zikov. Monitoring the anesthetic-induced unconsciousness (hypnosis)using wavelet analysis of the electroencephalography. Master’s thesis,University of British Columbia, 2002. → pages 2, 991Appendix APropofol PKPD ModelingThe pharmacokinetic (PK) and pharmacodynamic (PD) models are described here.PK models the distribution of the drug in the body and predicts the blood plasmaconcentration (Cp) of the drug. The PD models the effect of the drug from the drugplasma concentration. Propofol is the fastest anesthetic agent currently available[10]. A closed-loop control system with a fast reacting agent is more ideal. The PKand PD models discussed here are for this hypnotic agent.In this Appendix, a quick summary and the mathematical models are provided.More in depth discussion may be found in [43] and [10].A.1 PharmacokineticsThe drug uptake, distribution and elimination can be expressed mathematicallyby a pharmacokinetic model. There are a few models available; the exponentialand mamillary compartment models are the most common [10]. The mamillarycompartment model is discussed below.The body is divided into three compartments: 1) a central compartment con-sisting of the blood, brain and liver; 2) a larger compartment consisting of mus-cle and viscera; and 3) a third compartment consisting of bones and fat. The 3-compartment model is shown in Figure A.1.The drug is administrated into the central compartment intravenously. It iseliminated from the body according to the rate k10 through hepatic and/or renal92RapidlyequilibratingcompartmentCentralcompartmentSlowlyequilibratingcompartmentDrug AdministrationEliminationFigure A.1: The 3-compartment pharmacokinetic model. The rapidly equi-librating compartment models the muscles and viscera. The centralcompartment models the blood, brain and liver. The slowly equilibrat-ing compartment models the bones and fat.extraction. The concentration in the central compartment comes to an equilibriumwith the muscle-viscera compartment through the rate constant k21 (and the reverserate k12). The concentration in the central compartment also comes to an equilib-rium with the bones-fat compartment through the rate constant k13 (and the reverserate k31). These rates are usually provided in the units of min−1. The concentrationof the central compartment (C1) increases following the bolus but rapidly decreasesas the concentration of the muscle-viscera (C2) and bones-fat (C3) increases to bal-ance the equilibrium.The blood plasma concentration of the drug is the concentration of the centralcompartment; the mass-balance representation of this compartment in the state-93space is given by [10]:C˙1(t)C˙2(t)C˙3(t)=−(k10+ k12+ k13) k21 k31k12 −k21 0k13 0 −k31 ·C1(t)C2(t)C3(t)+1V1(t)00 · I(t),Cp(t) =[1 0 0]·C1(t)C2(t)C3(t)(A.1)where V1 is the volume of the central compartment; I(t) is the drug infusion rate;and by definition Cp(t) =C1(t) is the concentration of the central compartment.The clearance is the rate at which the drug is removed from a compartment andis expressed in [ml ·min−1]. The total body clearance Cl1 is given by:Cl1 =V1 · k10 (A.2)Likewise, the inter-compartmental clearances Cl12,Cl21,Cl13, and Cl31 are givenby Cli j =Vi · ki j. It is easy to realize that Cl12 =Cl21 =Cl2 and Cl13 =Cl31 =Cl3.The parameters ki j of the equation A.1 are computed according to the studypublished by Schu¨ttler et al. [48]. This population-based study relates the 3-compartment clearances and the volumes to the patient’s body weight and age aswell as the administration type (bolus vs infusion) and the sampling site (venousvs arterial). The values of the clearance and the volumes are shown in Table A.1.The estimates of the intermediates parameters of the Table are given in the TableA.2. The relationship of 3-compartment clearances and the volumes to the plasmaconcentration parameters are given in Equation A.3.k10 =Cl1V1k12 =Cl2V1k21 =Cl2V2k13 =Cl3V1k31 =Cl3V3(A.3)94Table A.1: Propofol PK parameters from [48]. BW stands for body weight,ven= 0 is for arterial sampling, ven= 1 is for venous sampling, bol = 0is for infusion administration, and bol = 1 is for bolus administration.Table A.2: Parameter estimates of the the PK model of Table A.1 from [48].The state-space representation of the PK model of A.1 can be given as a Single-Input/Single-Output (SISO) transfer function PK(s):PK(s) =Cp(s)I(s)=1V1· (s+ k12) · (s+ k31)(s+ p1) · (s+ p2) · (s+ p3) (A.4)where pi are the poles of the system.95A.2 PharmacodynamicsThe pharmacological response of a drug as a function of the drug plasma con-centration can be expressed mathematically by a pharmacodynamic model. Anydrug can target multiple organs in the body, resulting in multiple effects; there isnot a single unique pharmacological response for a given plasma concentration ofa drug. Here, the intake is the hypnotic drug propofol and the pharmacologicaleffect modeled is the depth of hypnosis.The full pharmacological model can be represented by a LTI transfer functionPD(s) and a non-linear sigmoid-type function known as Hill-equation [10].The LTI element models the effect site concentration Ce(s) from the dynamicsof the drug-receptor interaction. There will also be a delay for the drug to reachthe effect site from the plasma concentration. This LTI model is given by the firstorder time-delayed transfer function PD(s):PD=Ce(s)Cp(s)=12EC50· kds+ kd· e−Tds, (A.5)where EC50 is the plasma concentration which yields 50% of the maximum ef-fect; kd expresses the rate of the transfer of plasma concentration to the effect siteintroduced in [49]; and Td is the arm-to-brain delay.The non-linear Hill function models the dynamics of observed effect E(s) to theeffect-site concentration Ce(s). The observed effect E(s) runs from 0 (no hypnoticeffect) to 1 (fully awake). This model is given by:H(s) =E(s)Ce(s)= E0+Emax · Cγe (s)ECγ50+Cγe, (A.6)where E0 and Emax are the minimum and maximum effects, and γ is a measure ofthe steepness of the dose-response curve.96Figure A.2: The full PKPD model introduced in [10].PK(s) PD(s) H(s)ICp Ce EPharmacokinetics PharmacodynamicsA.3 The PKPD ModelThe complete model dynamics relating the drug administration I(s) of propofol tothe observed effect E(s) is shown in Figure A.2 can now be introduced:PKPD(s) =E(s)I(s)= PK(s) ·PD(s) ·H(s) (A.7)This structure represents a nonlinear SISO transfer function that relates the infusionof Propofol I(s) to the observed effect E(s). This effect can then be converted to aDOH index through a EEG monitor, such as the NeuroSense Monitor.The nonlinearity H(s) can be linearized around the DOH of interest. In mostcases, a DOH of 50 is considered adequate [39]. This value represents an E(s) =0.5 and is a logical operating point to linearize the Hill function around it. Thisassumption is only valid for the maintenance phase of anesthesia. The linearizedHill function around E = 0.5 is γ/2. The full detail of the linearizion can be foundin [43].The linearized PKPD(s) of A.7 can now be expressed as:PKPD(s) = Kpkpd · (s+ k12) · (s+ k31)(s+ p1) · (s+ p2) · (s+ p3) · (s+ kd) · e−Tds (A.8)where the patient’s model gain Kpkpd is:Kpkpd =kd · γ4 ·V1 ·EC50 (A.9)Bibian [10] analyzed the induction of 44 patients with a single bolus adminis-tration of propofol and measured the WAV. The PK and PD parameters of this studyare provided in Table A.3 and are used in Chapters 3, 4, and 5 as the model.97Table A.3: PK and PD parameters from the Bibian study [10].98Appendix BControl Performance inClosed-Loop AnesthesiaThe definition of Varvel and the proposed measures are discussed. In Figure B.1,an example of the DOH from a closed-loop control system is shown. The black linerepresents the induction phase, the blue represents the maintenance phase, and themagenta represents the emergence phase. Table B.1 displays the numerical valueof all the discussed measures for reference.B.1 Varvel MeasuresThe Varvel performance measures is constituted of 4 metrics: MDPE, MDAPE, Di-vergence, Wobble [65]. They are all based on the Percent Error (PE), defined as:PE = 100Cm−CpCp, (B.1)where Cm is the measured plasma concentration and Cp is the corresponding esti-mate of the Cm. Vectors could be in real-time (function of t) or in discrete (stepunit of h). The four Varvel metric measurements are now defined. In the context ofclosed-room control, Cm is usually replaced with the measured DOH y and the Cpis replaced by the set-point r.1. Median Performance Error (MDPE) measures the bias and is calculated as99the median of all the PE:MDPE = median(PE). (B.2)2. Median Performance Absolute Error (MDAPE) measures the inaccuracy andis calculated as the median of the absolute of PE:MDPE = median(|PE|). (B.3)3. Divergence measures whether the error is getting bigger or smaller as timeprogresses and is calculated as the slope of the linear regression of the abso-lute PE against time:MDPE =tT |PE|−nt¯|PE|‖t‖2−Nt¯2 , (B.4)where N is the size of the signal; t¯ is the mean of the signal; tT is the trans-pose of the vector; and ‖t‖2 is the square of the norm defined as tT t. Apositive Divergence signals an unstable control system.4. Wobble measures the variability of the estimator and is calculated as themedian of the absolute difference between PE and MDPE:MDPE = median(|PE−MDPE|). (B.5)In [39] a fifth parameter was introduced in an attempt to provide a single scalarscore to the overall performance of an EEG-guided DOH control system. The GlobalScore (GS) is then defined as:GS=MDPA+Wobblefraction of time DOH ∈ (40,60) . (B.6)The interval (40,60) for the DOH is clinically recommended [3] for maintenancephase of anesthesia.100B.2 Proposed Control Performance MeasuresThe proposed performance measures by Soltesz et al. [53] provide different met-rics for the three temporal phases of anesthesia: induction, maintenance, and emer-gence.A. Induction Phase Metrics1. Induction Phase Duration (ID) is adopted from [39]. Traditionally, IDis defined as the time it takes from the beginning of administration ofdrug to the time the DOH falls and remains below 60 for a duration of30 seconds. This definition only applies if the set-point (r) 50 is used.Instead, 60 is replaced by requiring the DOH value to be between r±10for 30 seconds.2. Percent Overshoot (OS) is defined as:OS= 100 ·min r− yE0− y , (B.7)where r is the reference, y is the measured DOH and E0 is the awakebaseline DOH. Typical value of E0 are in the range of 90 < E0 < 100.If E0 is not available, then the value 100 can be used.The maximum overshoot usually occurs after the end of the inductionphase. The signals r and y are then taken as the signals from the startof induction to 10 minutes after the end of induction phase.B. Maintenance Phase Metrics1. Integrated Error (IE) is introduced to replace MDPE (or the bias) of thesystem. It is calculated using the trapezoid approximation rule:IE =N∑k=1tk+1− tktN− t1 ·(rk+1− yk+1)− (rk− yk)2, (B.8)where the signals are for the duration of the maintenance phase only.The quantity tN − t1 is the length of the maintenance phase; IE is nor-malized to the length of the maintenance phase.1012. Integrated Absolute Error (IAE) is introduced to replace MDAPE (or theinaccuracy) of the system. It is calculated using the trapezoid approxi-mation rule:IE =N∑k=1tk+1− tktN− t1 ·|rk+1− yk+1|− |rk− yk|2. (B.9)3. Variability Index (VI) is introduced to replace Divergence of the system.It is calculated as the relative difference between IAE and IE:VI =IAE− IEIAE. (B.10)4. Percentage of Time Outside Adequate Range is calculated as the per-centage of the time the signal y is outside of the adequate range. Ade-quate range is defined as r±10 [39]. The sign of an error is of clinicalimportance. It is justified to provide two percentages, E+ for the timewhen r− y is more than +10, and E− when it is less than −10.C. Emergence Phase Metrics1. Emergence Phase Rise Time (ER) is the time it takes for the DOH toexceed r1+(1− e−1)(E0+ r1), where r1 is the set-point when the ad-ministration of the hypnotic drug was terminated. If the awake baselineE0 is not available, then E0 = 100 can be used as the default value.1020 50 100 150 200 2502030405060708090100DOH InductionMaintenanceEmergence0 50 100 150 200 2500200400600800100012001400Time (min)Infusion Rate (ml/h)Figure B.1: Representative example from a closed-loop DOH control sys-tem. The induction phase is shown in solid black, the maintenancephase is shown in solid blue, the emergence phase is shown in solidmagenta, the reference is shown in thick green, and the r± 10 boundsare shown in dashed black. The red dot represents the overshoot.Table B.1: Varvel and proposed measures of the example from Figure B.1.103Appendix CLimiting Behavior of L1 AdaptiveControlIn Chapter 3, it was claimed that the closed-loop response of the system G(s) ap-proaches that of the reference system Gre f (s) as Γ → ∞. Moreover, it was claimedthat the limiting controller is an implantable, non-adaptive LTI system. This im-plantable LTI system is independent of the system’s unknown parameters ω , θ , σ .This Section will show the proof of these claims.C.1 Problem Formulation and The L1 AdaptiveControllerThis section will provide a summary of the L1-AC. A more detailed description ofthe control structure is available in Chapter 3.Consider the following dynamics state-feedback controller G(s) (see Chapter2.2 of [26]):x˙(t) = Amx(t)+b(ω(t)u(t)+θT (t)x(t)+σ(t)),y(t) = cT x(t),(C.1)where x(t) ∈ Rn is the measured state of the system; u(t) ∈ R is the control input;y(t) ∈ R is the output; b,c ∈ Rn are assumed known constant vectors; Am is a104n× n Hurwitz matrix corresponding to the desired closed-loop dynamics; ω ∈ Ris an unknown constant but with known sign; θT (t) ∈ Rn is a vector of unknownparameters; and σ(t) ∈R models input disturbances. The dynamics of the desiredsystem M(s) are given by:x˙m(t) = Amxm(t)− kgbr(t), xm(0) = x0,ym(t) = cT xm(t)(C.2)where kg ,−1/(cTA−1m b) and r(t) is the reference.The state predictor is given by:˙ˆx(t) = Amxˆ(t)+b(ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ (t)),yˆ(t) = cT xˆ(t).(C.3)The adaptation laws are given by:˙ˆθ(t) =−Γ ·Proj(θˆ (t),−x˜TPbx(t)), θˆ (t) = θ0,˙ˆσ(t) =−Γ ·Proj(σˆ(t),−x˜TPb), σˆ(t) = σ0,˙ˆω(t) =−Γ ·Proj(ωˆ(t),−x˜TPbu(t)), ωˆ(t) = ω0,(C.4)where x˜(t) = xˆ(t)− x(t), Γ ∈ R+ is the adaptation gain, and P = PT > 0 is thesolution of the algebraic Lyapunov equation ATmP+PAm = −Q for arbitrary Q =QT > 0.The L1 control signal is defined as:u(s) =−kD(s)(ηˆ(s)− kgr(s)), (C.5)where r(s) and ηˆ(s) are the Laplace transforms of r(t) and ηˆ(t) respectively andηˆ(t), ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ(t). (C.6)k > 0 is a feedback gain and D(s) is a strictly proper transfer function such that105kg kD(s)ηˆ(t) = ωˆ(t)u(t) + θˆT (t)x(t) + σˆ(t)x˙(t) = Amx(t) + b(ωu(t) + θTx(t) + σ(t))y(t) = cTx(t)x˙(t) = Amxˆ(t) + b(ωˆu(t) + θˆTx(t) + σˆ(t))yˆ(t) = cT xˆ(t)˙ˆθ(t) = −ΓProj(θˆ(t),−x˜TPbx(t))˙ˆσ(t) = −ΓProj(σˆ(t),−x˜TPb)˙ˆω(t) = −ΓProj(ωˆ(t),−x˜TPbu(t))r +−+−Figure C.1: L1-AC as formulated in [26].they lead to a strictly proper stable filterC(s):C(s) =ωkD(s)1+ωkD(s). (C.7)C.2 Removal of the Internal Feedback over ηˆ(t)Figure C.1 shows the L1-AC architecture as shown in [26]. The internal feedbackover the signal ηˆ(t) is confusing and unnecessary. To analyze the limiting behaviorof the system, this loop needs to be taken out. The multiplication of the system’sadaptive parameters ωˆ , θˆ , and σˆ with the state feedback x(t) and the controller’soutput u(t) is the nonlinearity that is present in the control architecture. Using theformulation of ηˆ(t) from the Figure C.1, it follows that this signal be defined asthe output of the following dynamic system:˙ˆθ(t) =−ΓPro j(θˆ (t),−x˜TPbx(t)),˙ˆσ(t) =−ΓPro j(σˆ(t),−x˜TPb),˙ˆω(t) =−ΓPro j(ωˆ(t),−x˜TPbu(t)),ηˆ(t) = ωˆ(t)u(t)+ θˆT (t)x(t)+ σˆ (t).(C.8)The control architecture is a continuous system. The feedback signal u(t) overηˆ(t) from Figure C.1 is the same signal u(t) that feeds into the adaptation lawsblock. The same is true for the signal x(t) that feeds into the predictor and theadaptation laws block. It follows that ηˆ(t) is the only input signal to the predictor.106kg kD(s)x˙(t) = Amx(t) + b(ωu(t) + θTx(t) + σ(t))y(t) = cTx(t)x˙(t) = Amxˆ(t) + bηˆ(t)(t))yˆ(t) = cT xˆ(t)˙ˆθ(t) = −ΓProj(θˆ(t),−x˜TPbx(t))˙ˆσ(t) = −ΓProj(σˆ(t),−x˜TPb)˙ˆω(t) = −ΓProj(ωˆ(t),−x˜TPbu(t))ηˆ(t) = ωˆ(t)u(t) + θˆT (t)x(t) + σˆr +−+−Figure C.2: Equivalent architecture of the L1-AC with removed internalfeedback over ηˆ(t).This leads to a more compact block diagram shown in Figure C.2. The state-spacerepresentation of the predictor can now be replaced by its single input transfer func-tion representation H(s) since ηˆ(t) is the only input, andH(s) = (sI−Am)−1b. Thesame is true for the plant whose only input is u(t), and the state-space representa-tion of it can be replaced by its transfer function G(s). The block diagram is furthersimplified to Figure C.3.kg kD(s) G(s)H(s)˙ˆθ(t) = −ΓProj(θˆ(t),−x˜TPbx(t))˙ˆσ(t) = −ΓProj(σˆ(t),−x˜TPb)˙ˆω(t) = −ΓProj(ωˆ(t),−x˜TPbu(t))ηˆ(t) = ωˆ(t)u(t) + θˆT (t)x(t) + σˆr +−+−Figure C.3: Simplified architecture of the L1-AC to a more coherent struc-ture.107C.3 Linearizing the L1 Controller with GenericAdaptation LawsThe L1-AC control structure creates an internal feedback loop in the controller.The forward loop consists of nonlinear functions with the adaptive gain Γ, and thefeedback loop consists of the LTI function H(s). High-gain feedback over a plantresults in the approximate inversion of LTI function (see [23], Chapter 2.6). Thisinversion can be approximated with a linear model. The dynamic system in (C.8)is rewritten for generic adaptation laws and is shown in Figure C.4:˙ˆθ(t) = Γ f1(x, xˆ),˙ˆσ(t) = Γ f2(x, xˆ),˙ˆω(t) = Γ f3(x, xˆ,u),ηˆ(t) = g(θˆ , σˆ , ωˆ ,x, xˆ,u).(C.9)Taking the gain Γ out of the function fi clarifies the effect of increasing it. Let{xQ, xˆQ,uQ, θˆQ, σˆQ, ωˆQ, ηˆQ; t ∈ R} correspond to any set of equilibrium points ofthe closed-loop system, i.e. ˙ˆθ(t) = f1(xQ, xˆQ) = 0. Define:∆x(t) = x(t)− xQ,∆xˆ(t) = xˆ(t)− xˆQ,∆u(t) = u(t)−uQ,∆θˆ (t) = θˆ (t)− θˆQ,∆σˆ(t) = σˆ(t)− σˆQ,∆ωˆ(t) = ωˆ(t)− ωˆQ,∆ηˆ(t) = ηˆ(t)− ηˆQ,(C.10)where ηˆQ , g(θˆQ, σˆQ, ωˆQ,xQ, xˆQ,uQ).The linearization of nonlinear adaptation laws in equation (C.9) in close vicin-108kg kD(s) G(s)H(s)˙ˆθ(t) = Γf1(x, xˆ)˙ˆσ(t) = Γf2(x, xˆ)˙ˆω(t) = Γf3(x, xˆ, u)ηˆ(t) = g(x, xˆ, u)r +−+−Figure C.4: Simplified architecture of the L1-AC with generic adaptivelaws.kg kD(s) G(s)H(s)r +−+Fxˆ(s) Fx(s)Fu(s)+Figure C.5: Linearized L1-AC with generic adaptation laws.ity of the equilibrium states are given by a first order Taylor series:∆ ˙ˆθ(t) = Γ[∂ f1∂x∣∣∣xQ∆x(t)+∂ f1∂ xˆ∣∣∣xQ∆xˆ(t)],∆ ˙ˆσ(t) = Γ[∂ f2∂x∣∣∣xQ∆x(t)+∂ f2∂ xˆ∣∣∣xQ∆xˆ(t)],∆ ˙ˆω(t) = Γ[∂ f3∂x∣∣∣xQuQ∆x(t)+∂ f3∂ xˆ∣∣∣xQuQ∆xˆ(t)+∂ f3∂u∣∣∣xQuQ∆u(t)],∆ηˆ(t) =∂g∂ θˆ∣∣∣ηˆQ∆θˆ (t)+∂g∂ σˆ∣∣∣ηˆQ∆σˆ(t)+∂g∂ωˆ∣∣∣ηˆQ∆ωˆ(t)+∂g∂x∣∣∣ηˆQ∆x(t)+∂g∂u∣∣∣ηˆQ∆u(t).(C.11)109Let the notation Fix represent the Laplace transform of∂ fix∂x. Then:∆θˆ (s) = ΓF1x∆x(s)+ΓF1xˆ∆xˆ(s),∆σˆ(s) = ΓF2x∆x(s)+ΓF2xˆ∆xˆ(s),∆ωˆ(s) = ΓF3x∆x(s)+ΓF3xˆ∆xˆ(s)+ΓF3u∆u(s),∆ηˆ(s) = Gθˆ ∆θˆ(s)+Gσˆ ∆σˆ(s)+Gωˆ∆ωˆ(s)+Gx∆x(s)+Gu∆u(s),(C.12)whereGθˆ ,Gσˆ ,Gωˆ ,Gx, andGu are the Laplace transform of∂g∂ θˆ∣∣∣ηˆQ, ∂g∂ σˆ∣∣∣ηˆQ, ∂g∂ωˆ∣∣∣ηˆQ, ∂g∂x∣∣∣ηˆQ,and ∂g∂u∣∣∣ηˆQrespectively. Substituting the intermediate signals ∆θˆ (s), ∆σˆ(s), and∆ωˆ(s), the LTI system ∆ηˆ(s) can now be written as:∆ηˆ = Fx(s)∆x(s)+Fxˆ(s)∆xˆ(s)+Fu(s)∆u(s), (C.13)where Fx(s), Fxˆ(s), and Fu(s) are the linearized transfer functions between ∆x, ∆xˆ,∆u around the equilibrium points xQ, xˆQ,uQ, θˆQ, σˆQ, ωˆQ, ηˆQ and correspond to:Fx(s) = Γ(GθˆF1x+GσˆF2x+GωˆF3x),Fxˆ(s) = Γ(GθˆF1xˆ+GσˆF2xˆ+GωˆF3xˆ),Fu(s) = ΓGωˆF3u+Gu.(C.14)kg kD(s) G(s)H(s)r +−+Fxˆ(s)Fx(s)Fu(s)kgkD(s)−kD(s)++u xuxˆηˆFigure C.6: Equivalent form of Figure C.5 of the linearized L1-AC withgeneric adaptation laws.110kg kD(s) G(s)r +−Fx(s)1−Fxˆ(s)H(s)+kD(s)Fu(s)u xFigure C.7: Final form of the linearized adaptation laws for L1-AC withgeneric adaptation laws.The block diagram for this linearized controller is shown in Figure C.5. Thisblock diagram is equivalently shown in Figure C.6. This structure was realized bytaking out the signal u that feeds to Fu out of the internal loop. It now follows thatthe feedback controller, i.e. the relation between ∆x and ∆u, is derived by realizingthat ∆xˆ= H(s)∆ηˆ and ∆u=−kD(s)∆ηˆ and substituting ∆ηˆ from C.13:∆u(s) = K∆x(s)∆x(s) =−Fx(s)kD(s)1−Fxˆ(s)H(s)+Fu(s)kD(s)∆x(s). (C.15)This linearized LTI controller corresponds to a two-degree of freedom LTI con-troller. The transfer function between ∆r(t) and ∆u(t) has been omitted for sim-plicity, but it can also be derived easily. Figure C.7 shows the block diagram forthis LTI controller. The transfer functions Fx, Fxˆ, Fu can be replaced by their defi-nition C.14; the gain Γ can then be taken to infinity to yield the limiting behaviorfor the L1 controller with generic adaptation laws.In the next Section, the projection operator in the adaptation laws used in theHovakimyan’s implementation of the L1-AC is linearized and the limiting behavioris computedC.4 Linearization of the Projection Operator in the L1Adaptive ControlThe previous Section, the L1-AC with generic adaptation laws, was linearized. Inthis Section, the case of projection operator for the adaptation laws is linearized.111Following the linearization definition C.11, the adaptation laws C.4 are linearized:∆ ˙ˆθ(t) = Γ[− x˜Q(Pb)T∆x(t)+ xQ(Pb)T∆x(t)− xQ(Pb)T∆xˆ(t)],∆ ˙ˆσ(t) = Γ[(Pb)T∆x(t)− (Pb)T∆xˆ(t)],∆ ˙ˆω(t) = Γ[uQ(Pb)T∆x(t)−uQ(Pb)T∆xˆ(t)− x˜Q(Pb)T∆u(t)],∆ηˆ(t) = uQ∆ωˆ(t)+ωQ∆u(t)+ θˆTQ∆x(t)+ xTQ∆θˆ (t)+∆σˆ(t).(C.16)At the equilibrium points, it follows that f1(xQ, xˆQ) = 0, f2(xQ, xˆQ) = 0, andf3(xQ, xˆQ,uQ) = 0. This leads to xQ = xˆQ and x˜Q = xˆQ−xQ = 0. The system abovesimplifies to:∆ ˙ˆθ(t) = Γ[xQ(Pb)T∆x(t)− xQ(Pb)T∆xˆ(t)],∆ ˙ˆσ(t) = Γ[(Pb)T∆x(t)− (Pb)T∆xˆ(t)],∆ ˙ˆω(t) = Γ[uQ(Pb)T∆x(t)−uQ(Pb)T∆xˆ(t)],∆ηˆ(t) = uQ∆ωˆ(t)+ωQ∆u(t)+ θˆTQ∆x(t)+ xTQ∆θˆ (t)+∆σˆ(t).(C.17)Assume the initial conditions of these differentials are all zero, i.e. ∆θˆ0 = 0,∆σˆ0 = 0, and ∆ωˆ0 = 0. The system above can be written in the form of (C.13) as:∆ηˆ(s) =[Γs(xTQxQ+u2Q+1)(Pb)T + θˆTQ]∆x(s)−Γs(xTQxQ+u2Q+1)(Pb)T∆xˆ(s)+wQ∆u(s).(C.18)The transfer function between ∆x and ∆u is then written as:∆uˆ(s) =−kD(s)[Γs(xTQxQ+u2Q+1)(Pb)T + θˆTQ]1+ Γs(xTQxQ+u2Q+1)(Pb)TH(s)+ kD(s)ωQ∆x(s). (C.19)This controller is stable if and only if the L1-norm condition is satisfied and thelimit of the controller as Γ → ∞ exists. In this case, the limit of this controller112yields:ulim(s) = limΓ→∞−kD(s)[Γs(xTQxQ+u2Q+1)(Pb)T + θˆTQ]1+ Γs(xTQxQ+u2Q+1)(Pb)TH(s)+ kD(s)ωQx(s)=−kD(s)(Pb)T(Pb)TH(s)x(s).(C.20)113Appendix DRobustness and Performance ofiControlThe following data are from the 44 simulation case studies of Chapter 5. The ro-bustness, output disturbance rejection, and set-point response of the auto-tunedPID controller in response to a detected oscillation is compared to the current im-plementation of iControl. The PID parameters for both the auto-tuned controllerand the iControl are also provided.114Table D.1: Robustness comparison of the iControl vs the auto-tuned algo-rithm of Chapter 5 for the 44 PKPD models.115Table D.2: Output disturbance rejection comparison of the iControl tuningvs the auto-tuned algorithm of Chapter 5 for the 44 PKPD models.116Table D.3: Set-point response comparison of the iControl tuning vs theauto-tuned algorithm of Chapter 5 for the 44 PKPD models.117Table D.4: PID Parameters of the iControl tuning and the auto-tuned algo-rithm of Chapter 5 for the 44 PKPD models.118Appendix EOscillation Detection MATLABOscillationDetectionAlgorithm1 % Input ’ s order of s i gna l s :2 % 1. Output WAV3 % 2. Time ( second ) . W i l l conver t to min ; o s c i l l a t i o n measured i n minutes4 f u n c t i o n p m f = O s c i l l a t i o n D e t e c t i on A lg o r i t h m ( i npu t )56 % We requ i re to save a l l the output y7 % So we can go over i t and de tec t the o s c i l l a t i o n when enough data po in t s8 % are found9 p e r s i s t e n t s i g n a l ;10 g loba l toRunParam ;1112 % Reset on t ime = 0;13 i f i npu t ( 2 ) == 014 s i g n a l . Data = [ ] ;15 s i g n a l . Time = [ ] ;16 end1718 % Store the t ime and output y19 y = inpu t ( 1 ) ;20 s i g n a l . Data = [ s i g n a l . Data ; y ] ;21 s i g n a l . Time = [ s i g n a l . Time ; i npu t ( 2 ) / 6 0 ] ;2223 % This w i l l run the p a t i e n t case i n ” r e a l t ime ” to determine a l l24 % o s c i l l a t o r y components of the s i g n a l .25 % Please set a l l parameters i n the next sec t i on .26 % Time i s i n minutes !27 %% Test Parameters %%28 t e s t . per iod = 2 . 5 ; % Period Test11929 t e s t . magnitude = t e s t . per iod / 1 . 1 ; % Magntinude Test30 t e s t . f i t n e s s = 25; % Fi tness Test (%)31 %% Define the Parameters32 toRunParam . tes tParameters = t e s t ;33 toRunParam . windowPeriod = [ 1 : 3 0 ] ; % Period l l / u l bound ( i n minute )34 toRunParam . pe r i od Inc = 2; % The increment i n per iod from LL to UL35 toRunParam . windowSize = 10; % Number of per iods to be analyzed3637 toRunParam . mergePeriod = 0 . 5 ; % The d i f f e r e n c e between two consecut iveper iods to be considered same38 toRunParam . mergeTime = 10; % The d i f f e r e n c e between two consecut ive t imeto be considered same39 %% Loop through a l l requested o s c i l l a t i o n bounds to de tec t o s c i l l a t i o n40 % We need at l e a s t 2 data po in t s to c a l c u l a t e sampling t ime ( T 2 − T 1 )41 p m f = [0 0 0 ] ;42 i f ( l eng th ( s i g n a l . Time ) > 1)43 f o r windowPeriod=toRunParam . windowPeriod( 1 ) : toRunParam . pe r i od Inc : toRunParam .windowPeriod ( 2 )44 r = segmen t Osc i l l a t i on ( windowPeriod , s i g n a l ) ;45 i f sum( abs ( r ) ) ˜= 046 p m f = r ;47 end48 end49 end5051 end5253 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%54 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5758 %% Segment O s c i l l a t i o n59 f u n c t i o n p m f = segmen t Osc i l l a t i on ( windowPeriod , s i g n a l )60 % This w i l l f i n d ALL o s c i l l a t i o n s upto cu r ren t t ime61 % O s c i l l a t i o n s need to be merged a f t e r to remove r e p e t i t i o n62 % NOTE: We w i l l on ly go up u n t i l cu r ren t t ime . I f an o s c i l l a t i o n i s missed ,63 % i t w i l l no longer be detec ted .64 % WindowPeriod i s lower bound of the o s c i l a t i o n t h a t we need to de tec t . The65 % upper bound i s +toRunParam . pe r i od Inc6667 g loba l toRunParam sampleSet pat ientModel ;68 p e r s i s t e n t l a s t S e t ;6970 i f isempty ( l a s t S e t )71 l a s t S e t = −1;72 end1207374 Ts = s i g n a l . Time ( 2 ) − s i g n a l . Time ( 1 ) ;75 o s c i l l a t i o n S i g n a l s = [ ] ;76 per iodIn IndexLL = round ( windowPeriod / Ts ) ;77 periodInIndexUL = round ( ( windowPeriod+toRunParam . pe r i od Inc ) / Ts ) ;7879 iEnd = leng th ( s i g n a l . Time ) ;80 i S t a r t = iEnd − per iodIn IndexLL ∗ toRunParam . windowSize ;8182 p m f = [0 0 0 ] ;8384 % Enough data po in t to analze85 i f i S t a r t > 08687 s P a r t i a l = [ ] ;88 s P a r t i a l . Time = s i g n a l . Time ( i S t a r t : iEnd ) ;89 s P a r t i a l . Data = s i g n a l . Data ( i S t a r t : iEnd ) ;9091 % This runs the ODA and f i n d s the o s c i l l a t o r y p a i r92 O s c i l l a t i o n P a i r = ODA( s P a r t i a l , windowPeriod) ;9394 % An o s c i l l a t i o n was obsorved95 i f ( ˜ isempty ( O s c i l l a t i o n P a i r ) )96 t e s t = O s c i l l a t i o n P a i r . F i tness > toRunParam . tes tParameters . f i t n e s s ;97 rangeLL = windowPeriod <= O s c i l l a t i o n P a i r . Period ;98 rangeUL = O s c i l l a t i o n P a i r . Period < windowPeriod + toRunParam . pe r i od Inc ;99100 % The r e g u l a r t o r y t e s t s are passed101 % The per iod i s w i t h i n the LL & UL range102 i f ( t e s t && rangeLL && rangeUL )103 % Save the s t a r t / end t ime ind i ces104 % And save the data to the master s i g n a l105 O s c i l l a t i o n P a i r . i S t a r t = i S t a r t ;106 O s c i l l a t i o n P a i r . iEnd = iEnd ;107108 % Return the r e s u l t109 p m f = [ O s c i l l a t i o n P a i r . Period O s c i l l a t i o n P a i r . MagnitudeO s c i l l a t i o n P a i r . F i tness ] ;110111 %disp ( [ O s c i l l a t i o n P a i r . Period O s c i l l a t i o n P a i r . MagnitudeO s c i l l a t i o n P a i r . F i tness ] ) ;112 % New pa t i en t , need to create a blank s t r u c t113 i f ( l a s t S e t ˜= sampleSet )114 l a s t S e t = sampleSet ;115116121117 disp ( [ O s c i l l a t i o n P a i r . Period O s c i l l a t i o n P a i r . MagnitudeO s c i l l a t i o n P a i r . F i tness ] ) ;118119120 % Save t h i s ins tance of the model f o r f u t u r e use121 load ( ’ Osc i l l a to ryMode l s . mat ’ ) ;122 wr i teVar = s t r c a t ( ’ Mod i f i edPa t i en t ’ , num2str ( pat ientModel ) ) ;123 pa t i en tVa r = s t r c a t ( ’ PKPDPatient ’ , num2str ( pat ientModel ) ) ;124 eva l ( [ ’ g l oba l ’ pa t i en tVa r ’ ; ’ ] ) ;125 readVar = s t r c a t ( pat ientVar , ’ ( ’ , num2str ( sampleSet ) , ’ ) ’ ) ;126127128129 i f ˜ e x i s t ( wr i t eVar )130 eva l ( [ w r i t eVar ’ = [ ] ; ’ ] ) ;131 end132 eva l ( [ w r i t eVar ’ ( end+1) . t f = ’ readVar ’ . t f ; ’ ] ) ;133 eva l ( [ w r i t eVar ’ ( end ) .gamma = ’ readVar ’ .gamma; ’ ] ) ;134 eva l ( [ w r i t eVar ’ ( end ) . E0 = ’ readVar ’ . E0 ; ’ ] ) ;135 eva l ( [ w r i t eVar ’ ( end ) . bwt = ’ readVar ’ . bwt ; ’ ] ) ;136 eva l ( [ w r i t eVar ’ ( end ) . age = ’ readVar ’ . age ; ’ ] ) ;137 eva l ( [ w r i t eVar ’ ( end ) . bht = ’ readVar ’ . bht ; ’ ] ) ;138 eva l ( [ w r i t eVar ’ ( end ) . gdr = ’ readVar ’ . gdr ; ’ ] ) ;139 eva l ( [ w r i t eVar ’ ( end ) . study = ’ readVar ’ . s tudy ; ’ ] ) ;140 eva l ( [ w r i t eVar ’ ( end ) . PKtype = ’ readVar ’ . PKtype ; ’ ] ) ;141 eva l ( [ w r i t eVar ’ ( end ) . Td = ’ readVar ’ . Td ; ’ ] ) ;142 eva l ( [ w r i t eVar ’ ( end ) . EC50 = ’ readVar ’ . EC50 ; ’ ] ) ;143 eva l ( [ w r i t eVar ’ ( end ) . Kd = ’ readVar ’ . Kd ; ’ ] ) ;144145 save ( ’ Osc i l l a to ryMode l s . mat ’ , wr i teVar , ’−append ’ ) ;146 end147 end148 end149 end150 end151152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%156157 %% ODA %%158 f u n c t i o n o s c i l l a t o r y P a i r = ODA( s i gna l , windowPeriod)159 g loba l toRunParam ;160 % Computes the o s c i l l a t o r y p a i r161 % F i r s t i t detec ts a high / low p a i r162 % Then i t w i l l perform the t e s t to determine i f i t i s o s c i l l a t o r y .122163 % This w i l l r e t u r n the h ighes t f i t n e s s as the main o s c i l l a t i o n164165 Ts = s i g n a l . Time ( 2 ) − s i g n a l . Time ( 1 ) ;166 N = leng th ( s i g n a l . Time ) ;167168 x = s i g n a l . Data − mean( s i g n a l . Data ) ;169 y = dc t ( x ) ;170171 %% Get the SL Components172 Sy = std ( y ) ;173 yh = seaLevel (3∗Sy , y ) ;174 y l = seaLevel ( Sy , y ) ;175 Yi = ithDCT ( yh ) ;176 Yj = ithDCT ( y l ) ;177 Xi = i d c t ( Yi ) ;178 Xj = i d c t ( Yj ) ;179 [ dump, I ] = s i ze ( Yi ) ;180 [ dump, J ] = s i ze ( Yj ) ;181182 %% Find the pa i rs of x i and x j t h a t match up183 maxFitness = − I n f ;184 maxPair = n u l l ( 1 ) ;185 f o r i =1: I186 x i = Xi ( : , i ) ;187 y i = Yi ( : , i ) ;188 mi = max( abs ( y i ) ) ;189190 high = n u l l ( 1 ) ;191 low = n u l l ( 1 ) ;192 f o r j =1:J193 x j = Xj ( : , j ) ;194 y j = Yj ( : , j ) ;195 mj = max( abs ( y j ) ) ;196197 i f ( mi == mj )198199 high = generateSigna l ( s i g n a l . Time , x i , x ) ;200 low = generateSigna l ( s i g n a l . Time , x j , x ) ;201 break ;202 end203 end204205 % Now perform the t e s t s to see i f t h i s i s o s c i l l a t o r y206 i f ˜ isempty ( low )207208 % Only t e s t f o r the o s c i l l a t i o n s t h a t are w i t h i n the l i m i t209 rangeHighLL = windowPeriod <= high . per iod .mean ;123210 rangeHighUL = high . per iod . mean < windowPeriod + toRunParam . pe r i od Inc ;211 rangeLowLL = windowPeriod <= low . per iod . mean ;212 rangeLowUL = low . per iod . mean < windowPeriod + toRunParam . pe r i od Inc ;213214 i f ( rangeHighLL && rangeHighUL && rangeLowLL && rangeLowUL )215 hTest = high . per iod . t e s t >= toRunParam . tes tParameters . per iod ;216 l Tes t = low . per iod . t e s t >= toRunParam . tes tParameters . per iod ;217 i f ( hTest && lTes t )218219 % We only want the maximum f i t n e s s value .220 % So only s e l e c t t h i s p a i r i f the f i t n e s s i s the h ighes t value .221 i f ( low . f i t n e s s > maxFitness )222 p a i r = [ ] ;223 p a i r . high = high ;224 p a i r . low = low ;225226 maxFitness = low . f i t n e s s ;227 maxPair = p a i r ;228 end229 end230 end231 end232 end233234 % I f a maximum p a i r was found , then perform the magnitude t e s t .235 % We need to determine which component ( high or low ) to use f o r the236 % per iod / magnitude .237 % Whichever p a i r has the h igher p e r i o d i c r e g u l a t o r y value , w i l l then be238 % selec ted as the candidate .239 o s c i l l a t o r y P a i r = n u l l ( 1 ) ;240 i f ( ˜ isempty ( maxPair ) )241 high = maxPair . high ;242 low = maxPair . low ;243244 % Use the low component245 i f ( low . per iod . t e s t > high . per iod . t e s t )246 prd = low . per iod ;247 mag = low . magnitude ;248 % Use the high component249 else250 prd = high . per iod ;251 mag = high . magnitude ;252 end253254 % Fina l t e s t : magnitude r e g u l a t o r t e s t must a lso be s a t i s f i e d .255 mTest = mag. t e s t >= toRunParam . tes tParameters . magnitude ;256 i f ( mTest )124257 o s c i l l a t o r y P a i r = maxPair ;258 o s c i l l a t o r y P a i r . F i tness = low . f i t n e s s ;259 o s c i l l a t o r y P a i r . Magnitude = mag. mean;260 o s c i l l a t o r y P a i r . Period = prd . mean;261 end262 end263 end264265266267268269270271272273274275276277 %% Sea Level %%278 f u n c t i o n suppressed = seaLevel ( SL , func )279 % This f u n c t i o n supresses the values below Sy and re tu rns a vec to r o f280 % same diment ion as y , but w i th supressed values .281 N = leng th ( func ) ;282 tmp = zeros (1 ,N) ;283 index = f i n d ( abs ( func ) >= SL) ;284 tmp ( index ) = func ( index ) ;285 suppressed = tmp ;286 end287288 %% i t h D i sc re te Cosine Transform289 f u n c t i o n output = ithDCT ( y f )290 % Generates the i t h DCT component o f the vec to r sub jec t to the f o l l o w i n g291 % c r i t e r i a :292 % y i ( k ) = yf , i ( k ) f o r ks , i <= k <= ke , i ; o therwise 0293 % where294 % y f ( ks , i ) != 0 && y f ( ks , i−r ) = 0 f o r r = 1295 % y f ( ke , i ) ˜= 0 && y f ( ke , i + r ) = 0 f o r r = 1 ,2 ,3 ,4296 % ks , i <= ke , i297 % I t re tu rns a mat r i x o f k by N where k i s the number of segments t h a t298 % match the c r i t e r i a .299 N = leng th ( y f ) ;300 Yi = [ ] ;301 s = 2;302 whi le ( s ˜= N )303 i f ( y f ( s ) ˜= 0 && y f ( s−1) == 0 )125304 f o r e=s :N−4305 % match found . Find iDCT f o r the i t h component306 i f ( y f ( e ) ˜= 0 && leng th ( f i n d ( y f ( e+1:e+4) == 0) ) == 4 )307 y i = zeros (N, 1 ) ;308 y i ( s : e ) = y f ( s : e ) ;309 Yi = [ Yi y i ] ;310 s = e + 1;311 break ;312 end313 end314 else315 s = s +1;316 end317 end318 output = Yi ;319 end320321 %% Period Sequence322 f u n c t i o n T = periodSequence( time , func )323 % Calcu la tes the per iod sequence from the o r i g i n a l s i g n a l324 z = zeroCrossingSequence( time , func ) ;325 L = leng th ( z ) ;326 per iod = [ ] ;327 f o r l =1:L−1328 per iod ( l ) = 2∗( z ( l +1)−z ( l ) ) ;329 end330 T = per iod ;331 end332333 %% Period Test334 f u n c t i o n R = per iodTes t ( periodSequence)335 % Calcu la tes the per iod of a p a r t i a l ( iDCT) .336 % Based on the work Wang 2013337 alpha = 0.0027;338 N = leng th ( periodSequence) ;339 N = 8;340 CV = std ( periodSequence) / mean( periodSequence) ;341 x = ch i2 i nv (1−alpha /2 ,N−1) ; % We have df = L−1, and L =N+1342 f = s q r t ( x / ( N−1) ) ;343 R = f /CV;344 end345346 %% Magnitude Sequence347 f u n c t i o n magnitude = magnitudeSequence ( Ts , func , per iod )348 % Calcu la tes the Fi tness of a p a r t i a l ( iDCT) .349 % Based on the work Wang 2013126350 N = leng th ( func ) ;351 magnitude = [ ] ;352 I n t e r v a l = round ( per iod / Ts ) ;353 f o r l =1: I n t e r v a l :N−I n t e r v a l354 m = max( func ( l : l + I n t e r v a l ) ) − min ( func ( l : l + I n t e r v a l ) ) ;355 magnitude = [ magnitude m] ;356 end357 end358359 %% Magnitude Test360 f u n c t i o n R = magnitudeTest ( magnitudeSequence )361 % Calcu la tes the modi f ied r e g u l a t o r index .362 % R value > 2.73 denotes an o s c i l l a t i o n .363 alpha = 0.0027;364 N = leng th ( magnitudeSequence ) ;365 CV = std ( magnitudeSequence ) / mean( magnitudeSequence ) ;366 x = ch i2 i nv (1−alpha /2 ,N−1) ; % We have df = L−1, and L =N367 R = s q r t ( x ) / ( s q r t (N−1)∗CV) ;368 end369370 %% Zero Crossing Sequence371 f u n c t i o n z = zeroCrossingSequence( time , func )372 % Calcu la tes the zero−cross ing of a f u n c t i o n .373 N = leng th ( func ) ;374 z = t ime ( f i n d ( func ( 1 :N−1) .∗ func ( 2 :N) < 0 ) ) ;375 end376377 %% Fi tness Test378 f u n c t i o n F = f i t n e s s T e s t ( p a r t i a l , x )379 % Calcu la tes the Fi tness of a p a r t i a l ( iDCT) .380 % Based on the work Wang 2013381 F = 100∗(1−norm ( p a r t i a l − x ) / norm ( x ) ) ;382 end383384385 %% Generates a s i g n a l w i th a l l the needed components386 f u n c t i o n s i g n a l = generateSigna l ( t ime , p a r t i a l , x )387 % Charac ter izes the s i g n a l by d e f i n i n g388 % s i g n a l . x % time−domain s i g n a l389 % s i g n a l . y % DCT s i g n a l390 % s i g n a l . t ime391 % s i g n a l . maxDct392 % s i g n a l . magnitude % time−domain s i g n a l393 % s i g n a l . magnitude . s i g n a l394 % s i g n a l . magnitude . mean395 % s i g n a l . magnitude . s td127396 % s i g n a l . per iod % time−domain s i g n a l397 % s i g n a l . per iod . s i g n a l % time−domain s i g n a l398 % s i g n a l . per iod . mean % time−domain s i g n a l399 % s i g n a l . per iod . s td % time−domain s i g n a l400 % s i g n a l . zeroCrossing % time−domain s i g n a l401 % s i g n a l . index402 % s i g n a l . index . p403 % s i g n a l . index .m404 % s i g n a l . index . f405406407 per iod = [ ] ;408 per iod . s i g n a l = periodSequence( time , p a r t i a l ) ;409 per iod . mean = mean( per iod . s i g n a l ) ;410 per iod . s td = s td ( per iod . s i g n a l ) ;411 per iod . t e s t = per iodTes t ( per iod . s i g n a l ) ;412413414415 magnitude = [ ] ;416 magnitude . s i g n a l = [ ] ;417 magnitude . mean = [ ] ;418 magnitude . s td = [ ] ;419 magnitude . t e s t = [ ] ;420 i f ( ˜ isnan ( per iod . mean) )421 magnitude . s i g n a l = magnitudeSequence ( ( t ime ( 2 )−t ime ( 1 ) ) , p a r t i a l , per iod .mean) / 2 ;422 magnitude .mean = mean( magnitude . s i g n a l ) ;423 magnitude . s td = s td ( magnitude . s i g n a l ) ;424 magnitude . t e s t = magnitudeTest ( magnitude . s i g n a l ) ;425 end426427 s i g n a l . per iod = per iod ;428 s i g n a l . magnitude = magnitude ;429 s i g n a l . f i t n e s s = f i t n e s s T e s t ( p a r t i a l , x ) ;430 end128Appendix FPID Tuning AlgorithmPID Controller1 f u n c t i o n u r uP u I uD v = PIDforSwi tch ing ( s i gna l s )2 % Two DOF PID − same implementat ion as i n i C o n t r o l3 % Strange a n t i windup implementat ion . . .4 % Outputs the i n f u s i o n ra te i n ml / hr567 p e r s i s t e n t I y1 y2 r1 ysp18 p e r s i s t e n t K a dKick dKickmag ra1 rb1 rb2 K0 Ki0 Kd0 ;9 p e r s i s t e n t unstableParams switchedToUnstable ;1011 per iod = s igna l s ( 1 ) ;12 magnitude = s igna l s ( 2 ) ;13 f i t n e s s = s igna l s ( 3 ) ;1415 ysp = s igna l s ( 4 ) ;16 y = s igna l s ( 5 ) ;17 ub = s igna l s ( 6 ) ;18 l b = s i gna l s ( 7 ) ;1920 induct ionComplete2min = s igna l s ( 8 ) ;2122 r = ysp ;23 h = 5;2425 % I n i t i a l i z a t i o n26 i f isempty (K)27 load PIDparams . mat K a dKick ra1 rb1 rb2 uns tab le2812929 % Keep the ins tance of the o r i g i n a l s tab l e values f o r l a t e r comparison30 K0 = K;31 Ki0 = a ( 3 ) ;32 Kd0 = a ( 4 ) ;3334 I = 0;35 y1 = y ;36 y2 = 0;37 r1 = y ;38 ysp1 = y ;39 r = y ;40 ysp = y ;4142 % Unstable parameter swi tch43 swi tchedToUnstable = f a l s e ;44 unstableParams = uns tab le ;45 end46 i f dKick > 147 dKickmag = dKick ;48 dKick = 0 . 5 ;49 end50 i f dKick == 0 . 5 ;51 i f abs ( ysp − r1 ) > 0.0152 dKick = 0;53 y2 = dKickmag ;54 end55 end5657 % reference f i l t e r58 r = −ra1∗r1+rb1∗ysp+rb2∗ysp1 ;59 r1 = r ;60 ysp1 = ysp ;6162 % Measurement f i l t e r63 y2 = a ( 1 )∗y2 + a ( 2 ) ∗( y−y1 ) ;64 y1 = y1 + y2 ;6566 % Wait u n t i l i n d u c t i o n i s complete f o r 2min ,67 % Then swi tch to uns tab le68 i f ( induct ionComplete2min && ˜ swi tchedToUnstable )69 swi tchedToUnstable = t rue ;70 K = unstableParams ( 1 ) ;71 a ( 3 ) = unstableParams ( 2 ) ;72 a ( 4 ) = unstableParams ( 3 ) ;73 end7475 % CLP Dec 2 2009: Increased Cp l i m i t from 7 to 813076 % i f Cp>8, I =0;end % Ju l y 25 th 2012 , Klaske : This needs to be taken out ! ! !77 uP = K∗( r−y1 ) ;78 uI = I ;79 uD = −a ( 4 )∗y2 ;8081 % Ju ly 25 th 2012 , Klaske : Use u n f i l t e r e d y f o r p r o p o r t i o n a l e r r o r ?82 v = uP + uI + uD ;83 u = v ;8485 % Upper / lower l i m i t86 i f v < l b87 u = l b ;88 d i sp l ay ( ’ Lower Bound ’ ) ;899091 end92 i f v > ub93 u = ub ;94 d i sp l ay ( ’ Upper Bound ’ ) ;9596 end9798 I = I + a ( 3 ) ∗( r−y1 ) + a ( 5 ) ∗(u−v ) ;99100 % O s c i l l a t i o n Detected101 % Retune c o n t r o l l e r102 i f per iod > 0103 w = 2∗ p i / ( per iod ∗ 60) ;104105 % Descr ib ing Func t ion magnitude106 de l ta = ( ub − l b ) / 2 ;107 vv = v − de l ta ;108 i f abs ( vv ) <= de l ta109 N = 1;110 else111 alpha = as in ( de l ta / abs ( vv ) ) ;112 N = 1/ p i ∗(2∗ alpha + s in (2∗ alpha ) ) ;113 end114115 % New Stab le po in t116 Ms = 1 . 3 ;117 r s = (Ms−1) /Ms;118 ph i s = 2∗as in (1 / (2∗Ms) ) ;119120 % Current PID params121 P0 = −K;122 I0 = −a ( 3 ) / h ;131123 D0 = −a ( 4 )∗h ;124 i n i t s t a t e = [ P0 I0 D0 ] ;125126 % This c a l c u l a t e s what the PID i s at t h i s s ta te !127 [ dump, out ] = PIDTuningNLConstraints ( i n i t s t a t e , w, 5 , [0 0 ] ) ;128 r p = 1 / ( out ( 1 )∗N) ;129 ph i p = − out ( 2 ) ;130131 r c = r s / r p ;132 ph i c = ( ph i s − ph i p ) ;133134 opt ions = opt imset ( ’ MaxFunEvals ’ , 1e10 , ’ Max I ter ’ , 1e3 , ’ A lgor i thm ’ , ’i n t e r i o r−po in t ’ ) ;135 NLC = @ ( arg ) PIDTuningNLConstraints ( arg , w, 5 , [ r c ph i c ] ) ;136 r e s u l t = fmincon ( @PIDTuningObjective , i n i t s t a t e , [ ] , [ ] , [ ] , [ ] , [ 0 . 5 0.01 h] , [ 1 0 0.05 200] , NLC, opt ions )137138 disp ( r e s u l t ) ;139 pause140 end141142 disp ( ’−−−−−− ’ ) ;143144 % assemble output vec to r145 u r uP u I uD v = [ u r uP I uD v y1 y2 ] ;146147 endPIDTuningNLConstraints1 f u n c t i o n [ c , ceq ] = PIDTuningNLConstraints ( arg , w, N, c o n d i t i o n )2 %f u n c t i o n [ c , ceq ] = PIDTuningNLConstraints ( arg )3 %w = 2∗ p i / ( 3 8 ) ;4 %N = 5;5 %c o n d i t i o n = [2 .8982 4 .0464 ] ;6 % This w i l l r e t u r n the i n e q u a l i t y ( c ) and e q u a l i t y ( ceq ) c o n s t r a i n t s7 %8 % The inpu t args are Kp , Ki , Kd parameters9 % We w i l l rede f i ne K, Ti , Td to work w i th10 % The PID so l v i ng i s of the form :11 % U( s ) = K(1 + 1 / ( T i ∗ s ) + H( s ) ∗ Td ∗ s )12 % where13 % H( s ) = 1 / (1 + Kd / ( Kp ∗ N) ∗ s ) = Kp ∗ N/ ( Kp ∗ N + Kd ∗ s )1415 Kp = arg ( 1 ) ;16 Ki = arg ( 2 ) ;17 Kd = arg ( 3 ) ;1321819 K = Kp;20 Ti = Kp / Ki ;21 Td = Kd / Kp ;2223 alpha = N∗wˆ2∗Tdˆ2 / (Nˆ2+wˆ2∗Td ˆ 2 ) ;24 beta = w∗Td∗Nˆ2 / (Nˆ2+wˆ2∗Td ˆ 2 ) ;2526 rPa r t = 1+alpha ;27 i P a r t = beta − 1 / (w∗Ti ) ;2829 cmp = rPar t + i P a r t∗ s q r t (−1) ;3031 gain = K∗abs (cmp) ;32 phase = angle (cmp) ;3334 c = [ ] ;35 ceq ( 1 ) = gain − c o n d i t i o n ( 1 ) ;36 ceq ( 2 ) = phase − c o n d i t i o n ( 2 ) ;37 end133