ROBUST C O N T R O L : P1D vs. F R A C T I O N A L C O N T R O L DESIGN, A C A S E S T U D Y by Arturo Martinez B .Sc , Institute Tecnologico y de Estudios Superiores de Monterrey, Mexico, 2000 A THESIS SUBMITTED IN P A R T I A L F U L F I L L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF M A S T E R OF APPLIED SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES ( E L E C T R I C A L A N D C O M P U T E R ENGINEERING) T H E UNIVERSITY OF BRITISH C O L U M B I A November 2005 © Arturo Martinez, 2005 Abstract This thesis presents a systematic procedure to design PID and CRONE controllers to regulate the hypnotic state of anesthesia with the intravenous administration of propofol. The hypnotic state of the patient is assessed by means of the WA V index, a processed hypnotic index based on a deterministic analysis of the E E C The objective of the controllers is to provide an adequate drug administration regime of propofol, depending on the hypnotic set point, to avoid under or over dosage of the patients. Propofol is assumed to be administered by a commercial Graseby 3400 infusion pump. The two model-based controllers are designed to compensate for the patients inherent drug-response variability (uncertainty), to achieve good output disturbance rejection, and to attain good set point response. The drug-response model consists of two parts: a pharmacokinetic model that characterizes the drug movement in the body, and pharmacodynamic model that relates the drug concentration in the brain to the WAV hypnotic index. An anti-windup scheme is also implemented to protect the system against performance degradation in the event of actuator saturation. The performance of the controllers is assessed by calculating typical time domain measures (overshoot, settling and rise times), and using the median performance error, median absolute performance error, divergence, and wobble of the system's output. Overall, the CRONE controller was shown to have better output disturbance rejection, and to be less sensitive to the drug-response variability of the patients. ii Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii i Acknowledgments x 1 Introduction 1 1.1 Motivation 1 1.2 Objective and Scopes 3 1.3 Thesis Organization 5 2 Background 7 2.1 Drug Administration and Effect Modeling 8 2.1.1 Pharmacokinetics of Propofol 9 2.1.2 Pharmacodynamics of Propofol 10 2.2 Monitoring Depth of Anesthesia 11 2.3 Automatic Control of Depth of Anesthesia 12 3 Anesthetic Drug Modeling 15 3.1 Pharmacokinetic Models 15 3.2 Pharmacodynamic Models 19 3.3 Propofol PKPD Models 22 3.3.1 PKPD Linearization 23 3.3.1.1 Linearization for the Maintenance Phase of Anesthesia 25 3.3.1.2 Linearization for the Induction Phase of Anesthesia 25 iii 3.3.2 Optimal Nominal PKPD Models 26 4 PID Design 30 4.1 Controller Structure 31 4.2 Design Method and Optimization of Controller Parameters 34 4.3 Pre-Filter Design 37 4.4 Anti Wind-up implementation 39 4.5 Controller Parameters and Step Simulation Results 41 5 CRONE Design 48 5.1 Controller Structure 49 5.2 Design Method and Optimization of Controller Parameters 51 5.2.1 Optimization of the Generalized Template 52 5.3 Controller Frequency Domain Identification 54 5.4 Pre-filter Design 55 5.5 Anti-Windup Implementation 56 5.6 Controller Parameters and Step Simulation Results 56 6 Controller Performance Evaluation 63 6.1 Performance Evaluation Measures 63 6.2 Simulation Experiment 65 6.3 Controller Performance and Results 67 7 Conclusions 73 7.1 Summary and Contributions of the Thesis 73 7.2 Future Work 75 Bibliography 77 iv Appendix A : Laryngeal Mask Airway Study Population 82 Appendix B: Numerical Solution Procedure 84 B . l :PID Controller 84 B . 2: CRONE Controller 88 Appendix C: Set Point and Disturbance Response Analysis 94 C. 1: PID Response Data 94 C. 2: CRONE Response Data 98 Appendix D: Performance Simulation Plots 102 DA: PID Controller 102 D. 2: CRONE Controller 105 Appendix E: Acronyms 108 v List of Tables Table 3-1: Propofol PK parameters from [25], were BW=Body Weight, ven=\ for venous samples, ven=0 for arterial sample, bol=\ for bolus data, and bol=0 for infusion data.... 18 Table 3-2: Optimum identified nominal PKPD(s) parameters 29 Table 4-1: PID controller parameters 41 Table 4-2: PID step set point response (see Appendix C for complete data) 45 Table 4-3: PID output disturbance response (see Appendix C for complete data) 45 Table 5-1: CRONE design constraints 53 Table 5-2: Optimum CRONE design parameters (frequencies expressed in [rad/sec]).... 57 Table 5-3: Rational CRONE controllers 57 Table 5-4: CRONE step set point response (see Appendix C for complete data) 61 Table 5-5: CRONE step output disturbance response (see Appendix C for complete data). 61 Table 6-1: PID's predictive performance summary 68 Table 6-2: CRONE's predictive performance summary 68 Table A - l : L M A Population and PD parameters (Group 1 & 2). As presented in [8] 82 Table A-2: L M A Population and PD parameters (Group 3 & 4). As presented in [8] 83 Table C - l : Complete PID set point & output disturbance response data (Group 1 & 2). 94 Table C-2: Complete PID set point & output disturbance response data (Group 3 & 4). 95 Table C-3: Complete PID set point & output disturbance response data (Entire populations, patients 1-27) 96 Table C-4: Complete PID set point & output disturbance response data (Entire populations, patients 28-44) 97 vi Table C-5: Complete CRONE set point & output disturbance response data (Group 1 & 2). 98 Table C-6: Complete CRONE set point & output disturbance response data (Group 3 & 4). 99 Table C-7: Complete CRONE set point & output disturbance response data (Entire population, patients 1-27) 100 Table C-8: Complete CRONE set point & output disturbance response data (Entire population, patients 28-44) 101 vii List of Figures Figure 3-1: Three Compartment Pharmacokinetic Model 16 Figure 3-2: Pharmacodynamic Block Diagram (adapted from [8]) 20 Figure 3-3: PKPD block diagram. Adapted from [8] 23 Figure 3-4: Hi l l saturation linearization (Patient 1, y = 4.7) 24 Figure 3-5: Unstructured multiplicative uncertainty representation 27 Figure 4-1: PID closed-loop structure 32 Figure 4-2: Nyquist curve and robustness constraint (Group 1) 35 Figure 4-3: Simulink implementation of PID with Anti Wind-up 40 Figure 4-4: PID step set point and output disturbance response (Group 1) 42 Figure 4-5: PID step set point and output disturbance response (Group 2) 42 Figure 4-6: PID step set point and output disturbance response (Group 3) 43 Figure 4-7: PID step set point and output disturbance response (Group 4) 43 Figure 4-8: PID step set point and output disturbance response (Entire population) 44 Figure 5-1: CRONE closed-loop structure 49 Figure 5-2: Generalized template in the Nichols plane 50 Figure 5-3: Frequency response envelopes (Group 1) 54 Figure 5-4: CRONE controller Anti Wind-up scheme 56 Figure 5-5: CRONE step set point and output disturbance response (Group 1) 58 Figure 5-6: CRONE step set point and output disturbance response (Group 2) 59 Figure 5-7: CRONE step set point and output disturbance response (Group 3) 59 Figure 5-8: CRONE step set point and output disturbance response (Group 4) 60 Figure 5-9: CRONE step set point and output disturbance response (Entire Population). 60 viii Figure 6-1: Noxious Stimuli profile 67 Figure 6-2: PID performance experiment simulation (Entire population) 71 Figure 6-3: CRONE performance experiment simulation (Entire population) 71 Figure B-1: Graphical illustration of robustness constraint in Equation (B.3) 85 Figure D - l : PID performance experiment simulation (Group 1) 102 Figure D-2: PID performance experiment simulation (Group 2) 103 Figure D-3: PID performance experiment simulation (Group 3) 103 Figure D-4: PID performance experiment simulation (Group 4) 104 Figure D-5: PID performance experiment simulation (Entire population) 104 Figure D-6: CRONE performance experiment simulation (Group 1) 105 Figure D-7: CRONE performance experiment simulation (Group 2) 106 Figure D-8: CRONE performance experiment simulation (Group 3) 106 Figure D-9: CRONE performance experiment simulation (Group 4) 107 Figure D-10: CRONE performance experiment simulation (Entire population) 107 ix Acknowledgments This thesis has been completed with the support of the following people. First, I would like to thank my research supervisor Dr. Guy Dumont for his guidance, advice, and financial support through the duration of this thesis. I would also like to give special thanks to Dr. Mark Ansermino for always having time to answer my questions and giving great input to this work. Most importantly, I would like to thank my parents, Raquel and Maurilio, for their support and motivation during my early education years, my professional carrier, and the duration of my graduate studies at U B C . 1 owe everything to their love and support. In addition, I also thank the support of my two brothers Alfonso Martinez and Alberto Martinez. The encouragement of my family and other immediate relatives has made my graduate school experience enjoyable. Finally, I would like to thank Silvia Guzman, Jens Heymer, Chris Brouse, Sean Delfel, Tim Patterson, and the assistance of all the staff and students at the Pulp and Paper Centre. Thank you all! x Chapter 1 1 Introduction 1.1 Motivation These days undergoing a surgical procedure has become more and more common. These procedures range from life saving operations to physical appearance reconstruction. In the United Sates of America alone, 28 million patients are reported to undergo surgery annually [33]. Furthermore, it is estimated that the number of surgical procedures wil l increase from 28 million to 40 million per year over the next three decades [33]. In any surgical intervention, anesthesiologists play a very important role before {pre-operative), during (intra-operative), and after {post-operative) surgery [9], [10]. These phases include the assessment of the patient's medical condition and history before the operation, consultation with the surgical team on the requirements of the procedure, providing pain control and life support functions during surgery, supervising care and recuperation after surgery, and medically discharging the patient from the recovery room. The principal tasks of the anesthesiologist during surgery are to administer different kind of drugs to prevent the awareness of pain, attenuate the body's natural response to injury, and to provide optimal operating conditions for surgeons. In order to induce and maintain this adequate state of clinical anesthesia, anesthesiologists use a combination of hypnotic anesthetics, opioids and neuromuscular agents to induce a controlled state of hypnosis, analgesia, and muscular relaxation. These states are defined as the three functional components of anesthesia [29]. 1 Chapter 1: Introduction 2 Depth of anesthesia (DOA) during surgery can be assessed by evaluating the patient's clinical signs such as blood pressure, heart rate, respiratory depth or rate, lacrimation, movement, facial grimacing, and diaphoresis. Since these signs depend mainly on skeletal muscle activity, their reliability was lost after the introduction of neuromuscular agents. Therefore, monitoring of D O A relies on other physiological measures. Given that anesthetics, opioids and neuromuscular agents act via different mechanisms and influence different components of anesthesia, there is still no reliable single variable that can qualitatively assess the patient's DOA. A lot of research effort has therefore focused on assessing the hypnotic state of anesthesia. It is well known that certain intravenous anesthetics {propofol, thiopental, thiamylal, methahoxital, etomidate, etc...) produce a sedative-hypnotic effect through an interaction with the gamma-aminobutyric acid ( G A B A ) in the central nervous system (CNS) [29]. Therefore, it is possible to estimate the state of hypnosis by monitoring the patient's cortical activity. The ability of the brain to process information can be evaluated by the patterns of the electroencephalogram (EEG). Interpretation of the raw E E G pattern is very involved and requires a trained neurophysiologist for evaluation. These drawbacks have motivated the research community to use different signal processing techniques to derive a more concise measure of hypnosis. The E E G patterns have been analyzed using: time domain methods (zero-crossing frequency and mean frequency), frequency domain analysis (spectral edge frequency, median frequency, and relative delta power), Bispectral analysis (BIS™ monitor), and evoked potentials (Midlatency auditory evoked potentials, auditory steady state response, coherent frequency of auditory evoked response, and P300 component). Chapter 1: Introduction 3 See [4], [8], and [20]. Most recently, our research group at The University of British Columbia has developed a new technique that uses wavelet analysis of the E E C This novel Wavelet-based Anesthetic Value (WAV) compares very well with the Bispectral Index (BIS) to measure the state of hypnosis ([5], [6], and [8]). During the operation procedure, there are profound changes in the patient's condition due to surgical incisions or any other stimulating surgical event. Therefore, the anesthesiologist must carefully adjust the anesthesia regime of patients according to their medial condition, response to anesthesia, and requirements of the surgery. Their assessment of the clinical state of the patient relies on their experience, continuous monitoring and control of vital life functions (heart rate and rhythm, breathing, blood pressure, body temperature and body fluid balance), and interpretation of the E E G or its index for depth of hypnosis. Efficient and rapid evaluation of these signs is necessary to avoid over or under dose of the patient. Given the nature of the activities of the anesthesiologist during maintenance of anesthesia, it can be thought that he or she adopts the role of a feedback controller. Based on the patient's specific target value and monitor readings, adequate drug regime can be achieved by adjusting the settings of drug infusion devices and breathing systems. Therefore, it seems rational to tackle the problem from a control engineering perspective. 1.2 Objective and Scopes The objectives of this thesis are to design and evaluate the performance of two robust closed-loop controllers for the hypnosis state of anesthesia using the Wavelet-based Anesthetic Value (WAV) as a control variable. The aim is to provide an adequate drug Chapter 1: Introduction 4 regime for the administration of propofol, depending on the targeted depth of hypnosis, to avoid under or over dosing. There are four main design criteria for the controller: 1) It must be designed using well-known published pharmacokinetic models for intravenous administration of propofol [25], and clinically identified pharmacodynamic models using WAV as a measurement of hypnosis [8]. 2) It must be robust with respect to model uncertainty to compensate for the intra-and inter-patient variability during the maintenance of anesthesia. 3) It must minimize the response of the closed-loop system to output disturbances. To compensate for surgical incisions and noxious stimuli (i.e., stimuli associated with the transmission of nerve pain signals). 4) Finally, yet importantly, it must have good performance for set-point response. Even though the induction of anesthesia is often performed by the anesthesiologists, it is desired that the closed-loop system have little overshoot. It was found in literature that the majority of the controllers used in anesthesia are: fuzzy-controllers, Proportional-Integral-Derivative {PID) controllers, and Model Predictive Controllers (MPC). Most of these published articles do not explain the details of the controller's design criteria. In particular, it was found that a PID controller designed for anesthesia produced severe oscillations around the set point and its parameters were poorly tuned [1]. Therefore, we decided to investigate two robust loop-shaping techniques to design a PID, and a third generation CRONE controller. These controllers were chosen to compare the performance of a well-tuned low and high order controller. Chapter 1: Introduction 5 1.3 Thesis Organization The thesis describes in full detail the design and evaluation of the two proposed robust controllers. The remaining Chapters are organized as follows: Chapter 2. Background This Chapter gives a brief overview of the practice of anesthesia. The main concepts, terminology, and issues are covered. This is followed by a review of the pharmacokinetic (PK) and pharmacodynamic (PD) models found in literature, the tools for the measurement of hypnosis, and previous attempts of closed-loop control of anesthesia. Chapter 3. Anesthetic Drug Modeling This Chapter presents the background in how the PK and PD models are constructed. The main objective here is to adapt these models so they can be used in closed-loop control framework. This is the very first step necessary for a successful design of the proposed controllers. Chapter 4. PID Design This Chapter is devoted to the design of a robust PID controller. The main design criterion and structure is presented, followed by an explanation of the optimization procedure, pre-filter design, and the anti-windup implementation to compensate for saturation. Chapter 5. CRONE Design Similar to Chapter 4, the design criterion for the CRONE controller is presented. The optimization procedure, controller frequency domain identification, pre-filter design, and the anti-windup structure are also described. Chapter 1: Introduction 6 Chapter 6. Controller Performance Evaluation This Chapter elaborates on the results of the two proposed controllers. Furthermore, the predictive performance of the controllers is evaluated according to literature's methodology to test automatic control of anesthesia [32]. Chapter 7. Conclusions In this Chapter, the main contributions of the thesis are summarized and some related problems and future research prospective are presented. Chapter 2 2 Background During the past decades, engineering knowledge has been applied to solve many problems in medicine. For example, signal processing tools and techniques have been used to calculate the hypnotic state of patients from the patterns of the EEG, computer aided diagnosis, digital image processing, and robotics for tele-operated surgery are just a few. More recently, automation of anesthesia has caught the attention of the research community. The use of closed-loop controllers can improve the quality of drug administration. Closed-loop control can offer several potential benefits to the anesthesia practice. Because of more frequent sampling of the control variable (i.e., depth of hypnosis) and the more frequent adjustments in the infusion rate, automatic control can provide a drug administration regime that is adequate for each patient and avoid over or under dosage. Furthermore, controllers can take advantage of the drugs fast metabolism, elimination, and synergy to provide a smoother infusion type titration. Further, they can minimize the inter- and intra-patient variability to deliver an enhanced anesthesia profile according to the surgical stimulations of each procedure. By taking control of these routine tasks during surgery, the anesthesiologist can then concentrate on other signs related to the patient's safety. At the same time, the costs related to drug usage, bed occupancy, and recovery time can be significantly reduced. Due to these benefits, the clinical community has become more interested in the automation of anesthesia [16]. Even though there is still much research to be done, 7 Chapter 2: Background 8 literature shows there have been several works on modeling the effect of drugs, and monitor and control of DOA. To have a successful design of a closed-loop controller for DOA, it is necessary to have an accurate model of the process to be controlled (in this case the patient), and a reliable control variable that measures the relevant drug effect (i.e., hypnotic state). In the following sections, a brief literature review on the necessary components for automation of anesthesia is presented. First, modeling of the drug effect is described, followed by the different techniques for monitoring D O A . Finally, we wil l conclude this Chapter by reviewing prior attempts at closed-loop control of anesthesia. 2.1 Drug Administration and Effect Modeling Anesthesiologists use a combination of drugs to control the DOA of the patient. Each one of these drugs acts in different ways and has an effect on one of the three functional components of anesthesia (hypnosis, analgesia, and muscle relaxation). For instance, anesthetics are most often used to produce a state of hypnosis, while opioids and neuromuscular blocking agents are used to change the state of analgesia, and muscle relaxation, respectively. These drugs can be used in an intravenous or an inhaled titration. In this study, we will only focus on the intravenous administration by means of an infusion pump. It is necessary to understand the relationship between dose and pharmacological effect of the intravenously administered drugs. Each drug has different distribution, metabolism, and effect in the human body. Therefore, these models are drug dependant. The pharmacological study of the drugs is divided in two [28]. First, a pharmacokinetic (PK) model that relates the drug plasma concentration to the Chapter 2: Background 9 administered dose. Second, a pharmacodynamic (PD) model that relates the drug concentration to the observed effect. In the sections below, a review of published PK and PD models for Propofol are presented. 2.1.1 Pharmacokinetics of Propofol There have been a large number of PK models for Propofol published. A thorough review of these studies for the past 20 years can be found in [8]. However, the two most recent and accurate models are the ones identified by Schnider et al. [24], and Schuttler and Ihmsen [25]. A brief description of their work is discussed next. The study performed by Schnider et al. [24], focused on investigating the influence of method of administration (bolus vs. infusion), infusion rate, patient covariates, and disodium edetate (EDTA) on the pharmacokinetics of Propofol. At the time of this study, E D T A was a new addition to the formulation of Propofol emulsion to reduce bacterial growth. They analyzed the blood samples from 25 volunteers ranging from 26-81 years old. Each of the volunteers was administered a bolus dose of Propofol, followed by a 60-minute long infusion profile with different infusion rates. The entire population was studied twice, receiving either propofol with or without EDTA. Their results showed that the PK parameters are age, weight, height, and lean body mass dependant. Furthermore, the method of administration also influences the pharmacokinetics of propofol. On the other hand, E D T A did not show any changes, and the kinetics were linear with respect to the infusion rate. A set of equations is given to calculate the PK parameters of a three-compartment mamillary model. Chapter 2: Background 10 A more extensive study was performed by Schiittler and Ihmsen [25]. They analyzed 4,112 samples of 270 individuals collected from five research groups. The objective of the study was to estimate the pharmacokinetics of propofol with respect to the covariates age, body weight, and gender and to evaluate the inter- and intra-patient variability. The modes of administration (bolus vs. infusion) and sampling sites (venous vs. arterial) were also considered. Similar to [24], they published a set of equations to calculate the pharmacokinetic parameters of a three-compartment model. These results were shown to be age, body weight and administration type dependant, whereas the sampling site only had small influence on the parameters. A more detailed explanation of these models is given in section 3.1. 2.1.2 Pharmacodynamics of Propofol Pharmacodynamic models relate the drug plasma concentration to the pharmacological end effect. Since each drug can have different outcomes, PD models are identified to quantify the plasma concentration to a specific effect. For instance, given that propofol exerts its sedative-hypnotic effect by interacting with specific components of the G A B A receptors in the CNS, it is possible to quantify its hypnotic effect by evaluating any univariate index extracted from the E E G . Such as the Bispectral Index (BIS), Auditory Evoked Potential (AEP), and the Wavelet-based Anesthetic Value (WAV). A detailed presentation of the most recent PD studies can be found in [8]. The majority of the authors derived their PD models for the Bispectral Index (BIS), while only one study was done using auditory Evoked Potentials (AEP), and another with a proprietary E E G index based on semilinear canonical correlation. The major drawback in all of them is the limited number of patients analyzed. Chapter 2: Background 11 More recently, our research group has performed two studies to model the dynamics of thiopental [1], and propofol [8] TO. WAV, by analyzing the induction data for 5 and 44 patients, respectively. In both cases, the PD parameters of a first order plus time delay function and the hill saturation index were estimated using a classical least squares identification method. In contrast to other published work, a time delay is added to the model to account for the arm-to-brain travel time. The PD model identified in [8] wil l be explained further in section 3.2. 2.2 Monitoring Depth of Anesthesia Although as of today there is still no sensor that can effectively quantify D O A , recent work has demonstrated that derivatives from the E E G correlate well with the depth of hypnosis. For instance, time domain, frequency domain, Bispectral, wavelet, and evoked potential analysis are some techniques that have been used to derive a single index to measure and predict the depth of hypnosis (see [4] and [20]). In 1996, Aspect Medical Systems Inc. introduced to the market the BIS™ Monitor. This device records the E E G signal and displays a single value index (BIS index) representing the depth of hypnosis of the patient. BIS values are scaled between 100 and 0, where 100 is associated with the E G G patterns of an awake subject, and 0 represent an isoelectric E E G signal corresponding to the deepest level of hypnosis. So far, this is the only clinical monitor for hypnosis approved by the US F D A [8]. However, the BIS™ monitor has a significant time delay between the changes of the anesthetic state of the patient and the computed index. Although the bispectral index has been successfully used in the operating room and in the design of anesthesia controllers, a monitor that can react faster can further improve drug titration. Chapter 2: Background 12 To alleviate this problem, new signal processing techniques using wavelet analysis have been studied. Bibian et al. ([5], [6]) have proposed and patented a new hypnotic index based on the wavelet decomposition of the raw E E G signal. This technique, referred to as the Wavelet-based Anesthetic Value (WAV), has been shown to correlate very well with the BIS index. However, due to the relative lower computational complexity of wavelet analysis, the WA V index has proved to react faster to changes in the patient's state, reducing the time delay by up to 30 seconds [5]. This benefit can significantly improve the performance of any closed-loop controller, and on time identifications of over or under dosing. 2.3 Automatic Control of Depth of Anesthesia The application of automatic control in anesthesia is not new. Several feedback controllers have been presented in literature [26]. The design criteria and performance of such controllers strongly depend on the accuracy of the mathematical models of the patients, and the sensors used to assess the state of anesthesia. Ideally, the perfect feedback control system would be able to measure the three functional components of anesthesia (hypnosis, analgesia, and muscle relaxation), and steer the administration of multiple drugs simultaneously (i.e., anesthetics, opioids, and N M B ) . However, there are no reliable monitors to measure all these states, and no accurate mathematical models to relate the interaction between drugs. Therefore, most of the controllers designed so far are single input single output (SISO). In terms of automatic control of hypnosis, a large number of controllers have been published [4]. Each one of them mainly differs in their feedback quantity, and their control technique. A brief review of the most recent attempts is given next. ) Chapter 2: Background 13 Two feedback controllers based on midlatency auditory evoked potentials ( M L A E P ) have been published. These papers present two fuzzy rule-based closed-loop systems to control the inhalation concentration of isoflurane [21], and the infusion profile of propofol [17]. Although their performance showed good results when tested on dogs, no clinical studies with patients were made to further validate the proposed design. Furthermore, one major disadvantage of using M L A E P is its very low signal to noise ratio [4]. Frei et al. [11] published the design of a Model Predictive Controller (MPC). The main objective of this study was to control the Mean Arterial Pressure by regulating the inhaled concentration of isoflurane. The controller was tested on 20 volunteers and proved to perform better than manual control because of its more frequent sampling. Later on, following this same line of research, Gentilini et al. ([12], [13]) presented a cascaded internal model control technique to control hypnosis by means of the BIS index with isoflurane. The master controller provides end-tidal concentration references to the slave controller. These references depend on the difference of the measured and set point BIS values evaluated by the master controller. The proposed controller was shown to have good performance over all, and future work will consider the measurement of mean arterial pressure simultaneously. However, one of the drawbacks from this approach is the non-linearity's present in respiratory systems. Mortier et al. have published two articles ([19], [20]) explaining the benefits and results of a closed-loop controller for the administration of propofol. The system that was developed uses the BIS index as a control variable, and an adaptive model-based controller. The proposed controller was tested on ten patients undergoing surgery, and Chapter 2: Background 14 demonstrated that it is able to induce sedation with a minimum amount of drug. Unfortunately, due to patent protection the full details of the closed-loop algorithm are not disclosed. In 2001, Struys et al. [31] performed a comparative study between manual and automatic administration of propofol. The proposed closed-loop system consisted of an adaptive model-based controller guided by the patient's BIS index. The results obtained from the data of 20 patients indicated that the automatic administration of propofol reduced the recovery time compared to standard anesthesia practice. However, no details are given on the controller parameters and design criteria. Finally, Absalom et al. [1] have published the most recent attempt to control anesthesia. The authors of this paper present a PID controller for the intravenous administration of propofol using BIS as the control variable. The closed-loop system was tested on ten patients. Three of them showed severe oscillations around the targeted set point. This instability can easily be attributed to the poor tuning of the controller. As mentioned in the article, the controller parameters were taken from a different design using auditory evoked potential. Therefore, it is expected that the controller's performance will degrade i f the dynamics of the sensor change. Although there have been a large number of attempts to control anesthesia, none of them have convinced the clinical community of their safety and applicability. The controllers' performance must be extensively validated by means of simulations before they can be introduced in the operating room [30]. We believe the design procedure should be reevaluated from a control-engineering point of view. Chapter 3 3 Anesthetic Drug Modeling In this Chapter, we present the dose-response modeling of anesthetic drugs. The analysis of this relationship is divided into two models: a pharmacokinetic (PK) model that represents how the drug is distributed in the body, and a pharmacodynamic (PD) model that relates the drug concentration to the pharmacological effect. In the following sections, we will present the dose-response models to be implemented in a closed-loop control framework. First, in sections 3.1 and 3.2, the Laplace representation of the pharmacokinetic and pharmacodynamic models is covered, followed by the identification of the nominal PKPD(s) models to be used in the proposed controller designs. 3.1 Pharmacokinetic Models A pharmacokinetic model is a mathematical representation of the processes and rates of drug movement within the body. It relates the drug dose administration into the blood, distribution into the body tissues, and elimination by metabolism or excretion. For intravenous administration of anesthetics, this process can be analyzed using a three-compartment mamillary model as shown in Figure 3-1. Body tissues are grouped together into different compartments depending on their blood perfusion and drug solubility characteristics. In general, the body can be divided into three main compartments [28]. First, one small central compartment, termed the vessel-rich or visceral group, consisting of arterial blood and highly perfused tissues such as the brain, heart, kidney, liver, and endocrine glands. Second, a large lean group 15 Chapter 3: Anesthetic Drug Modeling 16 compartment containing muscles, viscera, and skin. Third, a vessel-poor group compartment containing mainly fat and bones. Drug (I) Lean Group * Compartment _ Visceral — Group — Central .4. Compartment V, k,0 Elimination C, k,3 C3 Vessel-Poor Group Compartment Figure 3-1: Three Compartment Pharmacokinetic Model. Immediately after a bolus intravenous injection, drug concentration in the central compartment, C\, rises rapidly, and three processes start simultaneously. First, the drug is distributed from the central to the two peripheral compartments at a rate of ki2 and ki3, respectively. This rate of movement is usually expressed in [min 1], and depends on the rate of blood flow trough the tissues, the mass of the tissues, and the characteristics of drug partitioning between blood and tissues. Second, after the drug concentration in all three compartments (C/, C2, and C3) reaches a state of equilibrium, the distribution process reverses and the drug returns back to the central compartment at a rate of k2i and kn. Finally, at the same time these two processes are taking place, the amount of drug injected into the central compartment is eliminated at a rate kw- This elimination process is usually achieved through metabolism and excretion. Chapter 3: Anesthetic Drug Modeling 17 Given that arterial blood is grouped into the central compartment, the drug plasma concentration time course (Cp(t)) is directly proportional to the concentration of the central compartment (Ci(t)). Thus, as presented in [8], the mass balance representation of the model can be written as: c,(0 c 2 (0 = C3{t) c,(0= [i • (kw + ki2 + £13) k2[ M2 1^3 "-31 2, 0 0 -k 31 C2(t) c3{t) + i(t) (3.1) •c,(/y C2(t) C3{t\ where Vi is the volume of the central compartment expressed in liters [1]. The parameters of equation (3.1) can be calculated according to the population pharmacokinetic models of propofol published by Schuttler and Ihmsen [25]. In this study, they present a set of equations to estimate the three-compartment clearances (Ch, Ch, and C/j) and volumes (Vs, V2, and V3), depending on the patient's body weight and age, and on the administration type (bolus vs. infusion) and sampling site (venous vs. arterial). The parameter set from this study is presented in Table 3-1. Pharmacokinetic volumes describe the extent of dilution of a drug in the body compartments. It is defined as the ratio of the total amount of drug in each compartment, divided by its initial concentration. Whereas clearance represents the ability of the body to eliminate the drug {Cli), and to exchange it between the central and two peripheral compartments (Ch and Ch). Chapter 3: Anesthetic Drug Modeling 18 Table 3-1: Propofol PK parameters from [25], were W = B o d y Weight, ven=\ for venous samples, ven=0 for arterial sample, bol=\ for bolus data, and bol=Q for infusion data. Model Parameters Value Units Cl, Ox ( B W / 7 0 )0 1 ifage<60 [1-min1] Ox { B W / 7 0 ) e i - {age - 6O)-0]O i f age > 60 [1-min1] Ch 03 • { B W / 7 0 f • (l + ven • 0 1 4)• (l + b o 1 ' °\6) [1-mm1] a3 Os • ( B W p O f ' -{l + bol-dj [1-min1] v, 02 • {BW/70f2 • {age/30)0" • (l + bol • 0l5) [1] v2 O, •'(BW/lOf' -(\ + bol-0„) [1] V3 [1] Parameter Estimates Value SE Units e, 1.44 0.09 [1-min1] 9.3 0.9 [1] 03 2.25 0.31 [1-min1] 64 4.2 6.1 [1] e5 0.92 0.15 [1-min1] e6 266 43 [1] e7 0.75 0.06 08 0.62 0.09 09 0.61 0.11 0/0 0.045 0.012 Sn 0.55 0.13 0,2 0.71 0.26 0,3 -0.39 0.15 014 -0.40 0.10 0,5 1.61 0.36 0,6 2.02 0.41 0,7 0.73 0.23 0,8 -0.48 0.12 Chapter 3: Anesthetic Drug Modeling 19 From these two independent parameters (Cl and V), the central and inter-compartment rate constants are calculated as: Cl Cl. k - c h k - C h k - c h JI y (3.2) where the units are expressed in [min 1]. For simplicity, the state space representation of the model expressed in equation (3.1) can be rewritten as a zero-pole-gain transfer function: PK(s) = KPK (s + z,)-(s + z 2) {s + p1)-(s + p2)-(s + p3) (3.3) where KPK represents the gain of the system, and p and z the poles and zeros, respectively. 3.2 Pharmacodynamic Models Pharmacodynamic models of intravenous anesthetics are mathematical expressions that quantify the pharmacological effect of a drug as a function of its plasma concentration. These models are not unique because all drugs have multiple effects and act on different target organs. Thus, PD models are constructed to predict only one endpoint of a specific drug. Chapter 3: Anesthetic Drug Modeling 20 For instance, anesthetic drugs such as benzodiazepines, barbiturates, and propofol are designed to act directly on the brain. These drugs exert their pharmacological effect by interacting with specific G A B A receptors in the CNS [29]. Therefore, their sedative-hypnotic effect is directly related to the drug concentration time course at this effect site (the brain). Although measuring absolute drug concentration at this effect site is not realistic, it is possible to quantify its pharmacological effect by interpreting the patterns of the EEG. For example, the E E G waveforms tend to slow down and increase in amplitude as the drug concentration in the brain increase. Furthermore, these changes can by analyzed using advance signal processing techniques to represent a single variable for hypnosis, (i.e., the WAV index, [8]). Non-Linear Function Scaling Function Pharmacodynamic Model Sensor Figure 3-2: Pharmacodynamic Block Diagram (adapted from [8]). The drug plasma concentration (Cp(t)) and the quantified single WA V variable (Eq(t)) relationship can be analyzed according to the block diagram shown in Figure 3-2. The SISO representation of this model can be divided into two functional groups: one pertaining to the pharmacodynamics of the drug itself, and another one to the dynamics of the sensing equipment. Chapter 3: Anesthetic Drug Modeling 21_ The pharmacodynamics of the drug can by represented by a linear time invariant (LTI) transfer function and a non-linear function. The LTI element accounts for the dynamics between the drug-receptor interaction, and the time the plasma concentration takes to reach its effect site. This PD(s) structure may be represented by any w-order plus dead time transfer function as: fD{s) = £ M = * . • * - + * - • f + • • • + V + V _ _ L _ , m < „ (3.4) where Td represents the arm-to-brain travel time expressed in [sec], and EC50 the steady-state plasma concentration necessary to obtain 50% of the measurable effect expressed in WlmL\-The non-linear function models the saturation between the observed effect (E0bs) and the plasma concentration at the effect site (the brain). The values of Eoos are bounded between 0 and 1, where 0 represent no hypnotic effect (state of a fully awake patient), and 1 is associated to the maximum effect of hypnosis that can be identified from the E E G patters. This saturation block is expressed using a Hi l l saturation equation: Er Enh = (3.5) obs 0.5r+E' where y represent the steepness coefficient of the saturation curve. Finally, the WAV index quantified effect Eq is obtained by considering the dynamics of the sensing equipment, H(s), and the scaling function. In this study the sensor transfer function H(s) is considered to by ideal, and the scaling function will only invert the observed effect value (E0bs) and scale it between 100-0 as: Chapter 3: Anesthetic Drug Modeling 22 Eq=l00-{\-Eobs) (3.6) where now 100 represents an awake state, and 0 the maximum level of hypnosis. 3.3 Propofol PKPD Models The drug-response relationship of propofol can be expressed by combining the pharmacokinetic and pharmacodynamic models presented in the previous sections. To construct these PKPD(s) models, we use the population demographics and PD(s) parameters identified in a Laryngeal Mask Airway ( L M A ) study presented in [8]. This study analyzed the induction data of 44 patients after a single bolus of propofol, and modeled the propofol infusion rate vs. the WA V index relationship. The whole population is divided into four age groups (see Table A - l and Table A-2) to assess the inter-patient variability of automatic infusion of propofol during the maintenance phase of anesthesia. In general, the PKPD(s) differences can be attributed to differences in age, lean body weight, previous diseases, and drug tolerance and sensitivity. The PK(s) part of the models is constructed as explained in section 3.1 by using each patient's age and body weight, and with the administration type and sampling site indicator variables (see Table 3-1) set for bolus data and arterial sampling, respectively. However, a high order PD(s) structure as presented in 3.2 could not be identified in the experiment due to limitations on the L M A data. Therefore, the PD(s) models are expressed as a simple first order plus dead time transfer function: PD(s) Er{s) (3.7) Cp(s) s + kd 2 EC. 50 Chapter 3: Anesthetic Drug Modeling 2 3 The complete PKPD(s) block diagram can be represented as in Figure 3-3 . This structure represents the nonlinear SISO dose-response characteristics of the patient, where Eq(s) is the observed WAV index hypnotic effect, and I(s) is the propofol infusion rate expressed in units of [mg/sec]. Non-Linear Function Scaling Function I ^ PK(s) C„ ^ PD(s) Er [mg/secJ [fxg/mLJ F [-] • [-J • W Sensor Pharmacokinetic Model Pharmacodynamic Model Figure 3-3: PKPD block diagram. Adapted from [8]. 3.3.1 PKPD Linearization The PKPD(s) models discussed in the previous section cannot be directly applied to a linear controller design technique because of the presence of the nonlinear function. Thus, it is necessary to linearize this element. In order to introduce the linearization procedure, let us first represent the Hi l l function as: F(x) = 0.5r +xy (3.8) where x is the input (Er) and F(x) the output (E0bs)- Next, consider the nonlinear gain characteristics of the function depicted in Figure 3-4, and assume that the steady state input to the nonlinearity is at the operating point JC0, as shown. Suppose that an infinitesimally small perturbation, AJC, occurs and now the updated input is at x = xa + A x . Then, the ratio of change in F(x) due to the change in x is approximately: Chapter 3: Anesthetic Drug Modeling 24 AF(JC)« ^ dF{x) dx • Ax (3.9) where the first derivative ofF(x) with respect to x evaluated atx 0 is: dF(x)\ dx r-i 0.5> (0 .5^+^) 2 (3.10) Following this approach, it is clear that the pharmacodynamic gain depends directly on observed effect. Thus, the hill saturation nonlinearity can be substituted by a linear gain around a specific operating point. The linearization procedure must be analyzed separately for the maintenance and induction phases of anesthesia. / ^y / Non-Linear Function F(x) /i fi ft /t /i /* ft Linear Approximation / V / / / / / / / / / / / / / / / / / / / X r -sr Figure 3-4: Hi l l saturation linearization (Patient 1, y = 4.7). Chapter 3: Anesthetic Drug Modeling 25_ 3.3.1.1 Linearization for the Maintenance Phase of Anesthesia. Optimum hypnosis levels during maintenance of general anesthesia are between 40-60% of the WAV index (note that before scaling, this range corresponds to an Er value between 0.6 and 0.4). Therefore, it is reasonable to choose a nominal operating point of 50%. By definition, this level of anesthesia is reached with a steady state plasma concentration of EC50, and as a result, the steady-state pharmacodynamic gain equals 0.5. Thus, evaluating equation (3.10) at this operating point, the linearized gain only depends on the hill steepness coefficient: (3-11) The linear PKPD(s) models can now be constructed by combining equations (3.3), (3.7), and (3.11), and represented as a zero-pole-gain transfer function: PKPlts) = ^ = PK(S).PD(s).t = e-<<is . K (* + *,)•(* + *,) ( 3 1 2 ) I{s) 2 \S + Py)-\S + p2)-[s + p3y\S + p^ where the overall PKPD(s) gain, before scaling, is expressed in [l/mg • sec"1 ] as: K m K " - 2 l ^ - 2 ( 3 ' l 3 ) 3.3.1.2 Linearization for the Induction Phase of Anesthesia The induction phase of anesthesia corresponds to the period when the patient is driven from a fully awake state to a desired D O A . This range can be quantified from values between 0 and 50 percent of the WAV index. Thus, for this case, the linear gain can be approximated as a straight line passing through the origin and the 50 percent WA lva lue . Chapter 3: Anesthetic Drug Modeling 26_ As explained in the previous, section, this point results on a pharmacodynamic gain of 0.5, and gives a constant slope of unity. Following this discussion, the PKPD(s) gain for the induction phase, or set point response, only depends on the pharmacokinetic and pharmacodynamic gains K * - K ~ - T k ; <3- l4) 3.3.2 Optimal Nominal PKPD Models The major difficulty in the design of automatic controllers of anesthesia is the inherent patient variability due to differences in demographics and drug tolerance. These discrepancies are translated into the PKPD(s) dose-response model uncertainty that characterizes the stability of the closed-loop system. A controller is considered robust i f it is insensitive to the differences between the actual patient's drug-response and the nominal model for which the controller was designed. Therefore, identification of a nominal model is mandatory before the two proposed controllers can be designed. From the patient data presented in Appendix A , we have decided to divide the PKPD(s) models into five design groups: four depending on age range, and one for the entire population. This division allows us to minimize the uncertainty bounds [3] and to study the performance of a controller designed for each group and for the entire population. For each group, we can represent the patient's drug-response relationship by a set of LTI models Gp (s) that correspond to each patient in that group, expressed as in equation (3.12), and a nominal model Gn(s) that characterizes that group's behavior without Chapter 3: Anesthetic Drug Modeling 27 uncertainty. In the frequency domain, the PKPD(s) uncertainty can be expressed using a lumped multiplicative structure of the form: Gp{s) = G0{s)\\ + w{s)A(s)) (3.15) which is represented by the block diagram in Figure 3-5. Here, A(s) is any stable transfer function that satisfies the condition | A ( / - © ) | < \ ---\fco and expresses the uncertainty as a unity disc-shaped region in the Nyquist plot, while w(s) is a weight function that quantifies the magnitude of the unstructured uncertainty. A more detailed explanation of uncertainty representation can be found in [27]. Figure 3-5: Unstructured multiplicative uncertainty representation. In this structure, the uncertainty at each frequency is represented as a disc-shaped region, centered at Ga(j • co) and of radius l(a>), that encircles all the parametric states of the plant. From equation (3.15), a rational representation of the weighting function is defined by: \w(j -o}\ = l{(o), Vo) (3.16) where Chapter 3: Anesthetic Drug Modeling 28_ A simple way of reducing the uncertainty expressed by equation (3.15) is to minimize \Gp(j ®)~G0(j-co^. This can be done by fine-tuning the parameters of G0(s), expressed as in equation (3.12), that minimize the difference between the patients and the nominal frequency response over a specific frequency range. This approach yields a realization of the same order as the original patients PKPD(s) models. Each group has as many models,Gp(s), as the number of patients. Thus, the minimizing cost function can be expressed as: J = i\Gn{j-co)-G0{j.a>l -±-<co< T T (3.18) where n is the number of patients in the group, and Tdmax and Tdmin are the maximum and minimum time delays from that group, respectively. The optimum parameters of the nominal model are bounded between: f < T°P< < T c/,min — 1 d — (/.max K < Kop' < K 7 < Z°P' < Z z l , m m - z \ — -M.max < zor" < Z z 2 ,min — ^2 - ^ 2,max / \ m i n ^ PT' ^ A .max P2 ,mi„ ^ PT ^ A . m a x />3 ,min ^ PT' ^ / \ m a x P4 ,min — P\ — P4,max (3.19) This optimization procedure is performed on each of the five groups mentioned earlier, and the identified models will be used later in Chapters 4 and 5 to design the two proposed controllers. A summary of the nominal parameters is presented in Table 3-2. Chapter 3: Anesthetic Drug Modeling Table 3-2: Optimum identified nominal PKPD(s) parameters. 29 Optimum nominal PKPD(s) parameters Group Td K z, Pi P2 Ps P4 [s] [10"4/mg-s_1] [io- 3 i [10"51 po- 2 i [lO"3] [lO"4] [10"51 1 18.6 1.698 1.477 2.572 3.239 6.961 2.803 2.703 2 16.5 1.928 1.478 2.703 3.843 7.735 2.912 2.787 3 8.3 1.438 1.486 3.627 2.870 7.748 2.843 2.121 4 17.8 1.823 1.478 2.651 3.656 9.100 2.962 2.710 A l l 16.4 1.580 1.479 2.577 3.467 8.230 2.962 2.785 Chapter 4 4 PID Design The Proportional-Integral-Derivative (PID) controller is used as a form of feedback for a wide variety of problems, such as process control, motor drives, magnetic and optic memories, automotive, flight control, instrumentation, etc. Its fixed complexity, only depending on three parameters, has proven to have satisfactory performance in many applications i f it is properly designed and well tuned. Furthermore, its implementation does not require very advanced equipment since it has been in the market for several years, and its operation principle is easy to understand. Therefore, the PID or PI is still the most widely used controller in the process industry. More recently, the PID controller has found its way into the anesthesia practice. Absalom et al. [1] published a study on a PID controller for the intravenous administration of propofol using BIS as the control variable (see section 2.3). The proposed controller was tested on ten patients, of which three showed continuous oscillations around the target set point. This instability can be attributed to the poor robustness of the controller. Furthermore, the article does not explain the design criteria used and how the gain and time constants were calculated. For this reason, we believe the PID controller design and structure should be revised using a control engineering approach. The systems under control are the dose-response models (PKPD(s)), explained in Chapter 3, for the L M A population data shown in Appendix A . A PID controller is designed for each of the nominal models identified in section 3.3.2 with the following performance specifications: 30 Chapter 4: PID Design 31 1) The closed-loop system must be robust with respect to model uncertainties to compensate for the variability of the dose-response characteristics of the patients. 2) The system must have a percent overshoot for set point and output disturbance response of less than 5%. This specification is to avoid over dosage. 3) In addition, it must have an initial rise time of around 5 minutes. This specification may be too strict since in common practice the operation procedure does not start until the patient reaches an adequate DOA, usually taking up to 15 minutes. Thus, a rise time between these values is considered appropriate. This Chapter presents a systematic way of designing a PID controller for the intravenous administration of propofol using the WAV index as a control variable. The closed-loop structure and the input-output relationships of the feedback system are described in section 4.1. An explanation of the design method and optimization of the controller parameters, as published in [23], is covered in 4.2. Section 4.3 describes the pre-filter design to improve the set point response of the closed-loop system, followed by the anti-windup implementation, as described in [2], in section 4.4. Finally, the controller parameters and step simulation results are explained in 4.5. 4.1 Controller Structure The closed-loop control structure for the PID design is depicted in Figure 4-1. Note that the plant transfer function, G0(s), represents the nominal PKPD(s) models identified for the five design groups (see section 3.3.2). The PID controller is implemented with two degrees of freedom characterized by two LTI transfer functions: Gc(s) describes the Chapter 4: PID Design 32 feedback from the process output y to the control signal u, and Gff(s) describes the feed-forward from the set point ysp to u. These two LTI transfer functions are expressed as: Gc(s) = kp+^ + kd-s s Gff(s)=kp+^ (4.1) Combining these equations describes the control law in the time domain as: «(')=*„ -MO-*('))+*, • \ ( y j T ) - y { T ) ) - d T + kd t d y ( ^ v dt j (4.2) Note from the equation above that the derivative term is only applied to the output y, this is usually done to avoid large transients in the control signal due to sudden changes in the set point. ysp- G0(s) Gc(s) 9f -> y Figure 4-1: PID closed-loop structure. A number of sensitivity functions can be calculated to describe the properties of the control system. These sensitivity functions describe how the input signals affect the process output. There are four external inputs acting on the closed loop system: set point ysp, load disturbances /, output disturbances d, and measurement noise n. The set point Chapter 4: PID Design 33 (ysp) or reference signal is the hypnotic level set by the anesthesiologist. Load disturbances (/) represent any variables that affect the control signal (w); these may be variations on the infusion rate due to obstruction of the IV tube. Output disturbances (d) are any external perturbations that affect the system output directly. During surgery, these are seen as sudden changes on the DOA of the patient caused by incisions, organ manipulations, body invasion, etc. Finally, measurement noise (n) accounts for the distortion of the feedback variable. The relationship between the four external input signals and the system's output is given by: y = G " ' G f f - y v + — ^ / + - d ^ n (4.3) \+G0Gc y s p \ + Ga-Gc l + G„-Gc \ + G0-Gc Note that in the previous equation the argument s was omitted to simplify the writing. A l l the transfer functions in equation (4.3) have the same denominator \ + G„(s)Ge(s)=\ + L(S) (4.4) where L(s) is the loop transfer function and can be analyzed using the Bode plot, Nyquist or Nichols diagrams to design a stable controller. The closed-loop transfer function from the output disturbance (d) to the output (y) is defined by the sensitivity function: Furthermore, note that by definition the sensitivity of the closed loop system to modeling errors is expressed as the maximum value of this transfer function: M, = max\s{j • co] (4.6) Chapter 4: PID Design 34 Thus, designing a controller with constraints on Ms will allow us to adjust the tradeoff between the performance and robustness of the controller, and to attain good output disturbance rejection. Typical values of Ms are in the range 1.2-2.0; where low values correspond to a robust controller with small percent overshoot, while high values would result in a more aggressive controller with an increased percent overshoot and less robustness. The design procedure is explained in the following section. 4.2 Design Method and Optimization of Controller Parameters The goal of the PID controller design procedure is to achieve a well-damped robust system with good rejection of load disturbances. The performance criterion is to minimize the Integrated Error (IE) at load disturbance. This calculation measures the system's ability to reject disturbances, where a small IE value results in no steady state error and fast rejection of disturbances. Furthermore, it has been shown that after a step disturbance is applied at the process input, the IE value is directly related to the PID's integral gain as [2]: iE = ){yspif)-y^))dt = T (4-7) 0 Ki Thus, the performance criterion of minimizing the integral error, IE, is equivalent to maximizing the integral gain ki. This criterion by itself does not guarantee the stability of the closed-loop system. The stability of the system can be analyzed trough the Nyquist plot of the open loop transfer function L(ja>), as shown in Figure 4-2. For a closed-loop system to be stable, it is desirable that L stays as far away from the -1 point as possible. The shortest distance Chapter 4: PID Design 35 from the Nyquist curve of the open loop transfer function to the critical point -1 is inversely proportional to the maximum of the sensitivity function, M~x. This may be interpreted as a robustness constraint that will maintain L\j • co) outside a circle centered at -1 and of radius \/Ms , that is: (4.8) Real Axis Figure 4-2: Nyquist curve and robustness constraint (Group 1). This constraint allows us to trade off performance and robustness by only varying one design parameter. Small values of Ms would result in a robust closed-loop system with Chapter 4: PID Design 36 low sensitivity to modeling errors, and as by the input-output relationship of equation (4.3), it wi l l correspond to a small output disturbance response. However, as explained in [23], for large values of the derivative gain, kd, the robustness constraint in equation (4.8) may give controller parameters that wil l make the Nyquist curve tangent to the sensitivity circle at two frequencies. This is mainly due to an excessive phase lead of the controller, and results in an oscillatory time response. To avoid these drawbacks two additional constraints are introduced into the design procedure: dco (4.9) < 0 (•2 , -2W with L{jco) = v{co) + jw{co) (4.10) where the dots denote the derivatives with respect to co. The first constraint in equation (4.9) prevents the controller from having excessive phase lead and the second constraint specifies that the Nyquist curve of L(j • co) should have negative curvature. The former becomes too strict for systems with integral or close to integral action; thus, this constraint wil l be omitted because the nominal dose-response models have poles relatively close to zero, as given in Table 3-2. For our implementation, the overall design procedure can be formulated as the parameter optimization that maximizes the integral gain kt, subject to the robustness and phase lead constraints in (4.8) and (4.9), respectively. Although the design procedure Chapter 4: PID Design 37 focuses on load disturbance rejection, it can also be applied for output disturbance rejection since it only depends on the maximum sensitivity value Ms. Hence, designing a controller with small Ms values will result on a closed-loop system with good output disturbance rejection and robustness to modeling errors. The details of the numerical solution procedure are briefly explained in section 1 of Appendix B. As mentioned at the beginning of this Chapter, one of the design objectives is to have an output disturbance response with an overshoot of less than 5%. This value corresponds to a damping ratio, £ , of 0.69. Furthermore, assuming the time domain specification can be investigated from a second order system, the maximum sensitivity value is expressed This value is used to design a PID controller for each of the nominal dose-response models presented in Table 3-2. A summary of the controller parameters and step simulation results follows in section 4.5. 4.3 P r e - F i l t e r Design In order to have good response to set point changes, a two degree of freedom controller is implemented to improve the tracking performance of the closed-loop system. This is done by adding a pre-filter compensator to the reference signal (ysp). The two degree of freedom transfer function relating the set point to process output is expressed by: as [2]: (4.11) y(s) _ G?{s)-Gff(s) •Fsp{s) = Gsp(s).Fxp(s) (4.12) ysp{s) I + G;(4GC(*) Chapter 4: PID Design 38_ where Gspp{s) are the PKPD models for all the patients in each of the design groups. These transfer functions are expressed as in equation (3.12) with the linearized set point gain calculated as (3.14), as explained in section 3.3.1.2. The pre-compensator is designed as a first order low pass filter, expressed as: K(s) =—-— (4-13) i + 5 - r s / , For each design group, there are as many Gp(s) models as patients in that particular group. To compensate for the worst-case set point response, the pre-filter time constant is tuned to the model that has the maximum resonance peak of all the Gsp(s) transfer functions. That is, the G™3" (s) that corresponds to the maximum M value in the design group. Msp =max|G™(/H V G » (4.14) According to the design specifications given at the beginning of the Chapter, the closed-loop set point response is desired to behave as a second order system with damping ratio of 0.69 and natural frequency of 1.07e"2 [rad/sec]. Hence, the pre-filter time constant can be calculated as: Tv=— (4-15) where a>cf is the cutoff frequency at which the magnitude ratio between the reference model, Tref, and the one degree of freedom set point response transfer function, G™"", crosses the -3dB point from above. Chapter 4: PID Design 39 Trefii-Vcf) = -3dB (4.16) 4.4 Anti Wind-up implementation A l l actuators have limitations on their output variables. These limitations may be due to safety restrictions imposed on the actuators to prevent damages to the system, or to inherent physical restrictions of their operating principle. When an actuator reaches a limit, its output level is saturated and no attempt to further increase the control signal wil l give any variations on the plant's input. The main effects of actuator saturation on a feedback system are large percent overshoots, damped oscillations, and even instability to disturbances. Most controllers are implemented with integral action to eliminate steady state errors. However, when an actuator reaches a hard limit and its output saturates, the error continues to be integrated and the state of the integral term becomes very large. This effect is called an integrator wind-up, and results on large transients of the output signal. This problem can be solved by implementing an anti wind-up scheme (see Figure 4-3). This is done by turning off, or resetting, the integrator whenever the output of the actuator saturates. The integral action is reset with a time constant Tt, depending on the difference between the output of the controller (w) and the actuator output (w). When no saturation effect is present, the difference between these signals is zero and thus the controller operates normally. However, i f the actuator saturates the difference is not zero and the integrator's input becomes zero. The time constant T, is calculated as [2]: (4.17) Chapter 4: PID Design 40 C D -Ysp Set point Y Output Gain Integral Gain *9" Integral Reset Time Constant / l / T t \ C D -w Actuator Traking Proportional Gain Kds Kd/(Kp*N).s+l Derivative Filter Integrator - • C D u Control Signal Figure 4-3: Simulink implementation of PID with Anti Wind-up. Note in Figure 4-3 that the derivative term is implemented as a filter. This is usually done to limit the effect of high frequency noise. Typical values of N are in the range 2-20. For our implementation, this value is set to ten. For the automatic infusion of propofol, there are two strict saturation levels. First, once a drug is injected into the patient it is impossible to subtract it, thus an obvious safety limit is to restrict negative infusion rates. This sets a lower saturation limit of zero. The second saturation limit is imposed by the maximum infusion rate that the infusion pump can deliver. We assume that a Graseby 3400 infusion pump is used for this purpose. According to the manufacturer's specifications, this device can pump up to 1200 [mL/Hr] [15]. Hence, with an emulsion containing 10 [mg/mL] of propofol [18], the Chapter 4: PID Design 41 upper saturation limit of the actuator corresponds to a maximum propofol infusion rate of 3.33 [mg/sec]. 4.5 Controller Parameters and Step Simulation Results Following the procedure explained in the previous sections, a PID controller is designed for each of the five nominal models shown in Table 3-2. These models represent the linearized nominal dose-response characteristics of the age groups and the entire study population. The controller parameters are summarized in Table 4-1. Table 4-1: PID controller parameters. Group PID Parameters K */ kd T * sp T t 1 2.610 0.026 65.09 156.81 49.819 2 3.947 0.046 85.29 129.90 43.024 3 10.207 0.107 202.38 111.95 43.397 4 4.455 0.058 104.83 124.96 42.416 A l l 4.455 0.058 104.83 124.96 42.416 To assess the performance of the controllers, the closed-loop system response is simulated for a 50% step set point change at time zero, and a 20% step output disturbance applied 2000 seconds (33.33 minutes) after the beginning of the simulation (see Figures 4-4 through 4-8). Note that in these plots the infusion rate is expressed in [mg/Kg/min], which is the most commonly used unit by anesthesiologists. For the simulations, the nonlinear function in equation (3.5) was implemented; that is, the controllers designed for the linearized models were tested with the nonlinear hill equation. Chapter 4: PID Design 42 100 80 I60 > 40 20 0 -^ -Patients . Y sp Disturbance 10 21) 30 40 Time |minutes| 50 60 30 40 Time [minutes) 50 60 Figure 4-4: PID step set point and output disturbance response (Group 1). 10 20 30 40 Time [minutes] — Patients Y sp Disturbance 50 60 30 40 Time |minutes| Figure 4-5: PID step set point and output disturbance response (Group 2). Chapter 4: PID Design 43 100 p 80 - 6 0 : > 401 Patients 20 h - — Y sp Disturbance 1 L 1 1 - I I 10 20 30 40 50 60 Time jminutes] Time jminutes] Figure 4-6: PID step set point and output disturbance response (Group 3). i i — — T R— • —| 1 i Patients Y sp Disturbance 1 1 1. . i 10 20 30 40 50 60 Time [minutes] — 4 0 10 20 30 40 50 " 60 Time |minutes| Figure 4-7: PID step set point and output disturbance response (Group 4). Chapter 4: PID Design 44 0 10 20 30 40 50 60 Time |minutes | Figure 4-8: PID step set point and output disturbance response (Entire population). To evaluate the closed-loop system response of the controller, we calculate the following four time domain measures. The system's set point response and disturbance rejection characteristics are evaluated by the maximum overshoot (OS) and a 2% settling time (TS). Although the controllers were specifically designed for a 5% overshoot, we consider here any value below 10% to be acceptable. To analyze the design specifications on the set point rise time (7V), we define this as the time required for the response to rise from zero to the desired reference level for the first time. Finally, the Integrated Absolute Error (IAE) is calculated to assess how well the controller corrects for an output disturbance. A summary of these measures is given in Table 4-2 and Table 4-3 for set point and output disturbance response, respectively. Chapter 4: PID Design 45 Table 4-2: PID step set point response (see Appendix C for complete data). Set Point Response Group T [min] T * s [min] OverShoot [%] min max mean min max mean min max mean i 3.99 10.71 6.92 3.88 13.49 7.22 0.65 13.86 3.08 2 2.79 9.05 6.85 4.06 27.19 7.40 0.65 36.05 4.05 3 6.23 7.88 7.22 5.20 5.60 5.42 0.43 0.90 0.62 4 2.84 8.94 5.99 3.26 8.39 5.63 0.64 3.55 1.58 A l l 2.71 10.66 7.87 4.31 33.18 8.18 0.33 37.07 3.30 Table 4-3: PID output disturbance response (see Appendix C for complete data). Output Disturbance Response Group IAE [-] T [min] OverShoot min max mean min max mean min max mean 1 13.97 46.07 24.24 5.97 16.85 8.89 3.88 11.95 5.97 2 15.67 30.65 22.09 5.48 12.29 7.51 3.44 16.93 6.50 3 10.44 22.34 15.75 1.60 8.14 4.12 0.78 3.23 2.12 4 13.21 34.87 24.71 4.87 12.61 8.24 3.62 8.99 6.11 A l l 10.93 42.62 22.09 3.97 14.60 7.58 3.15 19.90 6.20 Even though the dose-response characteristics between patients vary extensively and that the controllers were tested with the presence of the nonlinear Hi l l equation, it is interesting to note that a simple and well-tuned PID controller gives good output disturbance rejection. The main objectives during the anesthesia procedure are to maintain the patient at a desired D O A and to avoid over dosage. From the simulation results shown previously and the data in Table 4-3, it is noticed that all the controllers in each design group maintain an average output disturbance overshoot, or overdose, of less Chapter 4: PID Design 46 than 6.5% and a settling time under 8.89 minutes. Furthermore, the output disturbance responses are stable and well damped with no severe oscillations. The controllers designed for Groups 3 and 4 have the best performance, as shown in Figure 4-6 and Figure 4-7. These controllers behave very well for each patient in their group because the pharmacodynamic parameters are closer to each other (see Table A-2). Thus, there is less PKPD uncertainty and the nominal models capture more adequately the dose-response characteristics. Note as well that for Group 3 the rise time is higher than the settling time. This is mainly due to how we defined Tr, and only means that the system is within 2% of the desired set point before it crosses it for the first time (see Table 4-3). The responses for Groups 1 and 2 show a higher set point and disturbance overshoot for patients #3 and #16, respectively. As shown in Figure 4-4, Figure 4-5, Table 4-2, and Table 4-3. The initial response for patient #16 has a decaying oscillation around the set point and an unacceptable overshoot of 36% (see Table C - l ) . This behavior is due to the patient's long pharmacodynamic time delay shown in Table A - l . This value is almost three times the nominal time delay presented in Table 3-2, and therefore explains the increased settling time and overshoot. However, besides the responses of these two patients, all the set point and output disturbance overshoots of the two Groups are below the maximum allowed value of 10% (see Table C- l ) . The controller designed for the entire population has only four patients with unacceptable overshoots (see Figure 4-8, Table C-3 and Table C-4). Besides the responses for these four patients, the controller behaved well in spite of the very high uncertainty due to the variability in the patients' dose-response characteristics. Chapter 4: PID Design 47 Furthermore, i f we ignore patients #1, #3, #7, and #16, the data in Table C-3 and Table C-4 reveal that the maximum set point overshoot is only 5.32% with an average of 1.18%. Applying the same outlier discrimination analysis for the output disturbance, the data reveals a maximum overshoot of 9.12%) with an average of 5.50%. Thus, the controller satisfies the design specification for 40 out of 44 patients on the design for the entire population. Chapter 5 5 CRONE Design In this Chapter, we present the CRONE controller design strategy. CRONE is the French acronym for 'Commande Robuste d'Ordre Non Entier'. This control system design, based on fractional differentiation, is a frequency domain approach for the robust control of perturbed plants using a unity feedback configuration. The design procedure is divided into two parts. First, for the nominal state of the plant the open loop transfer function is correctly positioned in the Nichols diagram to guarantee the design specifications and to maintain robust stability with respect to plant uncertainties. Second, the optimum controller is then obtained from the ratio of the optimal open loop to the nominal plant transfer functions. Three CRONE control generations have been developed and applied to a variety of problems. In this thesis, we design and implement a third generation controller as explained in [22]. The design specifications are the same ones used for the PID controllers (see Chapter 4). The closed-loop structure and input-output relationships of the feedback system are described in section 1.1. The design method and optimization of the controller parameters are explained in section 5.2, followed by the controller's frequency domain identification in section 5.3. The pre-filter design is covered in section 5.4, and the anti-windup implementation, as described in [14], in section 5.5. Finally, the controller parameters and step simulation results are explained in section 5.6. 48 Chapter 5: CRONE Design 49 5.1 Controller Structure The two degree of freedom CRONE closed-loop structure is shown in Figure 5-1. Four external inputs are acting on the system. The relationship between these signals and the system's output is given: G0C 1 + G -C A sp Ssp ^ 1 -d — 1 + G • C 1 + G • C 1 + G • C (5.1) where G0(s) represents the nominal PKPD(s) models (see section 3.3.2), C(s) is the CRONE controller, and F (s) is the pre-fdter compensator. FsP(s) C(s) o -P(s) G0(s) / ! + Figure 5-1: CRONE closed-loop structure. The open loop transfer function here is defined as: B{s) = G0{s)-C{s) (5.2) The third generation CRONE method is based on a generalized template in the Nichols plot. This template is defined by a straight line of arbitrary inclination around the open loop ultimate frequency cou. The shape of the generalized template is characterized by a Chapter 5: CRONE Design 50 transfer function of complex non-integer order of integration, n = a + i-b, where the real part (a) determines the phase location at frequency cou, that is at a-nil, and the imaginary part (b) determines the inclination with respect to the vertical, as shown in Figure 5-2. Negative values of b correspond to a positive slope of the straight-line segment, while positive values result in a negative slope. In needs to be specified that the imaginary unit /, which is in the integration order n = a + i-b, is independent of the imaginary unit j which is in the operational variable s = a + j -a>. 0 argfTQ-m) Figure 5-2: Generalized template in the Nichols plane. The generalized template is described by a transfer function with complex non-integer order of integration: B{s) = K cosh V ^ J) CO.. V •> J ( f cos V v b-ln\ V s J)) -sign(b) (5.3) The position and orientation of the template can be set within a frequency range [ ^ , < w „ ] by varying the parameters of n. Thus, an optimal template is defined as the generalized template that: Chapter 5: CRONE Design 51 • Is tangent to an overshoot (Qd) or damping (Q) contour, and • Minimizes a quadratic criterion based on the extreme variations of Qd or Q due to the uncertainty from the perturbed states of the plant. 5.2 Design Method and Optimization of Controller Parameters In the version of the third generation CRONE design used in this thesis, the open loop transfer function for the nominal state of the plant, B0(s), takes into account the control specifications at low and high frequencies, the generalized template around unit gain frequency cou, and the non-minimum phase part of the PKPD(s) models. Thus, Ba(s) is defined as [22]: 1 + -COu 1 + ^ -1+-c 1 + -cos 6-In V v 1 + C •-COu 1 + -co * J) \ 2 b J 1 + COu 1 + ^ - 1 + -(5.4) where C = 1 + 2 \ 1 + ^ coD (5.5) Chapter 5: CRONE Design 52 K ensures a unity gain at frequency cou, nb and nh are respectively the asymptotic behavior orders in open loop at low and high frequencies, cor is the resonance frequency, and co'h, cob, coh, and co'h define the transitional frequencies. The time delay of the nominal dose-response models is substituted by a second order Pade approximation. This approximation computes a rational transfer function with two stable poles and two unstable zeros. Thus, the nominal model is divided into its minimum and non-minimum phase parts. G„{s) = Gmp{s)-Gnmp(s) (5.6) Note, however, that the performance limitation of the nominal models, before the Pade approximation, is directly related to the time delay of the model. This value provides an upper bound of a>u <\jTd on the achievable open loop ultimate frequency [27]. The maximum time delay of the entire population would correspond to an ultimate frequency of 2.22e"2 [rad/sec] (see Table A - l ) . For our CRONE controller designs, we have set a>u at 1 .Oe'2 [rad/sec] for the five design groups. 5.2.1 Optimization of the Generalized Template The open loop Nichols locus optimization procedure determines the optimal values of the following seven parameters: the optimal gain K, the optimal real and imaginary orders of integration, a and b , and the optimal transitional frequencies co'b, cob, coh, and co'h. However, due to the conditions of unit gain frequency at CDU and tangency to an overshoot contour, the computation is reduced to only five independent parameters [22]. Chapter 5: CRONE Design 53 The optimal template is achieved by a constrained optimization procedure. The optimal template minimizes the quadratic criterion: J CRONE ~ ("rmax ~ 4d} + (^min _ £ d ) (5-7) where, according to the design specifications of Chapter 4, the desired damping ratio, % d, is equal to 0.69. The optimization procedure is subject to the constraints shown in Table 5-1, where these values were fine-tuned by trial and error. The control effort is limited to avoid saturation of the infusion pump for a 20% output disturbance (see section 4.4). nb and nh are fixed to 1 and 2, respectively. Table 5-1: CRONE design constraints. Complementary Sensitivity Function T(s) G (s) max|r0--^<l 0<co <1.96e~3 Limits sluggishness ^ m a x W 1.57e"2 < co < oo Limits cutoff frequency Gmin(s) mm\T(j-colB>-\ 0<ft><1.57<T3 Limits measurement noise effect Sensitivity Function S(s Gmax(s) max\s(j-co]dij<-5 0 < co < le'3 Improves output disturbance rejection. G.{s) max\S(j-colB<2.\6 V<y Limits Sensitivity to modeling error. Input Sensitivity Function CS(s) G m a x V * ) m a x | C S ( 7 - ^ < 24.43 Vcy Limits control effort In these constraints, Gmax(s) and Gmin(s) are the frequency response envelopes for each design group, as shown in Figure 5-3. The details of the numerical solution procedure are briefly explained in section 2 of Appendix B. Chapter 5: CRONE Design 54 10 o -10 -20 -30 h |G I 1 m a i l to*.*****, Frequency [rad/sec] 0 -50 „ -100 DC -150 X -200 a £ -250 -300 -350 10 arg(G ) ^ o v max' -arg(G . ) o x mm' arg(G o ) *X^V^I \ \ --i 10 Frequency | rad/sec | Figure 5-3: Frequency response envelopes (Group 1). 5.3 Controller Frequency Domain Identification Once the optimal open loop template has been determined, the frequency response of the fractional controller Cf(s) is calculated as: CAs) (5.8) where the optimal frequency response of Bo{s) is obtained from equation (5.4), and G0(s) describes the nominal dose-response characteristics as in equation (5.6). The synthesis of the rational transfer function Cr(s) consists of identifying the parameters of a rational integer model that has the same frequency response as C f (j • a>). The rational model on which the parameter estimation is performed is given by: Chapter 5: CRONE Design 55 Cr (*) = - H = - v — ^ ! 2 , "c< dc (5.9) A[s) s ' + ad^ • s ' +--- + a]-s + a0 The only restrictions on the controller are that it must be realized as a minimum phase and biproper transfer function; that is, it should be stable and have a constant high frequency gain different from zero. This restriction is only imposed to enable implementation of the controller in an anti-windup structure (see section 5.5). The controller parameters were identified using the Frequency Domain System Identification Toolbox in Matlab. A summary of the controller parameters and the step simulation results is given in section 5.6. 5.4 Pre-filter Design The pre-filter is designed in the same way as explained in section 4.3. The only difference is that the closed loop transfer function relating the set point to the process output is now expressed as: jfel_= gyfrKfr) ,F w up w ( 5 I 0 ) yr(s) \ + G?(,)-C,{s) F " [ S > °»WA,W (5-10) Note that this equation, in comparison to (4.12), only depends on the rational representation of the CRONE controller Cr(s), the pre-filter transfer function F (s), and the PKPD models with the linearized gain for the induction phase of anesthesia represented by Gpp(s). Chapter 5: CRONE Design 56 5.5 Anti-Windup Implementation The anti-windup scheme used for the CRONE controller is shown in Figure 5-4. Given that the controllers identified in section 5.3 are biproper and of minimum phase, it is possible to split Cr{s) into a direct feed-through term c„ and a strictly proper transfer function C(s) [14]: C{s) = cK+C{s) (5.11) c* Saturation i „ -i Figure 5-4: CRONE controller Anti Wind-up scheme. Here, the saturation block represents the lower and upper infusion rate limits of the infusion pump (see section 4.4). This structure resets the controller's integrator whenever a saturation effect occurs on the output of the actuator. 5.6 Controller Parameters and Step Simulation Results Following the procedure explained in the previous sections, a rational CRONE controller is designed for each of the design groups. The optimal open loop Nichols locus is obtained by using the linearized nominal dose-response model in Table 3-2, with the constraints and frequency response envelopes explained in section 5.2.1. The optimized CRONE parameters are shown in Table 5-2. Chapter 5: CRONE Design 57 Table 5-2: Optimum CRONE design parameters (frequencies expressed in [rad/sec]). CRONE Design Parameters Group a b K <o't (Ob (Or (»h (O'h T •* sp [ lO 4 ] [lO"4] [lO"3] [ i f / 1 ] [lO"1] 1 1.118 -0.238 34.641 3.987 4.496 3.987 6.646 23.441 88.80 2 1.124 -0.258 31.846 4.301 5.329 4.295 3.194 3.194 67.65 3 1.040 -0.336 32.324 3.511 4.102 3.511 0.351 0.351 49.16 4 1.125 -0.226 31.960 4.272 6.181 4.272 2.863 3.213 72.37 A l l 1.078 -0.201 37.645 3.418 3.418 3.092 1.228 14.697 144.85 The rational transfer function representation of the CRONE controller is then identified as in section 5.3, and given in Table 5-3. The degree of the rational controller is not the same for each group. During the identification procedure, we tried to approximate the ideal frequency response of Cf (j • a>) by the lowest order rational transfer function. Table 5-3: Rational CRONE controllers. Group CRONE Controller Cr(s) 1 6.31-/+1 L54-J7 + 9 . 2 1 / +8.37-/ +1.43 / +5.62-e"2 -s3 +A2\-e*-iL +\2-e1 s+Z.02en ss +1.93- s7+1.76- /+1.5- s5 +0.45- s4 +3.21- e2 • s3 +4.64- es • s2 +1.19- e9 • s-27& e22 2 4 . 6 4 s 6 +3 .66 -S 5 + 5 . 7 7 - i 4 + 031-s3 +3.49-e" 3 -s 2 + 1.04-<T6 -s + 2J6-e~u s6 +0.86 -5 5 +1.28 / +0.17 -5 3 +2.64-e" 4 -s 2 + 735-e~9 - 5 - 2 . 2 4 - e " 2 1 3 1.9-s6 + 0 . 3 1 - / + 9.81-e~3 -s 4 + 7.08-<T5 s3 + 4.99<T8 V +9.22-e" 1 2 -5 + 1.73-e-'6 s6 +0.15-5 5 + 3.28e" 3 -s4 + 6.65e~6-s 3 +3.06e" 9 -s 2 +1.02-<T13 s-2.49-e"27 4 6.76 -s6 + 4 . 5 8 - / + 279.8 • s4 + 18.69 -s3 + 0.19 • s2 +5.94 e - 5 -s + l.5\-e^ s6 + 0.78 • s5 +41.44 -s4 + 7.42 • s3 +1.05 • e'2 • s2 + 2.82 • e~7 • s - 4.46 • e~ 2 0 A l l 4 . 5 2 - / +0.64-s4 +2.28-e" 2 -s 3 +\J9-e~4-s2 + 5.43e~ 8 s + ].38-e-u s5 + 0 . 1 9 - / + 8 .79e" 3 • s3 +1.23 e'5 • s2 +3.13• e" 1 0 • s -3 .59• <T23 Chapter 5: C R O N E Design 58 The same simulation profile and performance measurements used section 4.5 are repeated here to assess the closed loop response of the C R O N E controller. The controllers are simulated with the non-linear hill equation. Figures 5-5 to 5-9 show the step set point and output disturbance responses for each group. The simulations reveal the improved robustness of the C R O N E controller over the P I D controller. Furthermore, the C R O N E controller reduces the infusion rate and thus also the saturation of the actuator. A summary of the set point and output disturbance performance measurements is given in Table 5-4 and Table 5-5, respectively. The complete data are available in Appendix C. — 4 — Patients Y sp Disturbance 10 20 30 4 0 Time |minutes| 30 40 T i m e [minutes| 50 60 Figure 5-5: C R O N E step set point and output disturbance response (Group 1). Chapter 5: CRONE Design 59 i n L Patients . Y sp Disturbance I ! 1 1 1 1 L _ 0 10 20 30 40 50 60 Time (minutes) Time [minutes] Figure 5-6: CRONE step set point and output disturbance response (Group 2). Disturbance ol 1 1 1 1 1 1 0 10 20 30 40 50 60 Time |minutes] Time |minutes| Figure 5-7: CRONE step set point and output disturbance response (Group 3). Chapter 5: CRONE Design 60 Time [minutes] Figure 5-8: CRONE step set point and output disturbance response (Group 4). oi 1 1 1 1 1 1— 0 10 20 30 40 50 60 Time (minutes| Time |ininutrs| Figure 5-9: CRONE step set point and output disturbance response (Entire Population). Chapter 5: CRONE Design 61 Table 5-4: CRONE step set point response (see Appendix C for complete data). Set Point Response Group [min] T [min] OverShoot [%] min max mean min max mean min max mean 1 4.00 7.44 5.64 5.56 12.57 8.33 0.96 13.63 4.04 2 3.16 6.42 4.34 4.49 10.17 7.57 1.45 32.03 6.98 3 2.93 6.14 4.12 4.22 9.86 6.84 0.00 12.35 6.02 4 3.05 6.10 4.64 3.44 10.19 7.20 0.57 12.81 4.70 A l l 3.88 16.39 14.31 6.02 13.01 9.62 0.08 11.67 0.71 Table 5-5: CRONE step output disturbance response (see Appendix C for complete data). Output Disturbance Response Group IAE [-] T [min] OverShoot [%] min max mean min max mean min max mean i 14.84 46.05 24.27 2.52 9.51 5.31 0.33 3.56 1.60 2 17.66 34.69 23.85 3.16 7.88 5.12 0.11 8.80 2.17 3 22.04 57.89 36.55 2.53 18.45 7.84 0.00 1.92 0.45 4 13.93 36.44 25.66 4.38 7.36 6.35 0.42 4.41 1.84 A l l 12.05 44.05 24.36 2.22 10.65 5.45 0.18 9.58 2.12 The output disturbance rejection of these controllers is improved considerably. Each group maintains an average overshoot, or over dosage, of less than 2.17% and a settling time below 7.84 minutes (see Table 5-5). These values clearly show that the design specifications on disturbance rejection and robustness are satisfied for all the design groups. Furthermore, the disturbance responses are stable and well damped with no significant oscillations. The IAE values are slightly higher in comparison to those of the Chapter 5: CRONE Design 62 PID controller; this can be attributed to the reduced overshoot of the CRONE (compare Table 4-3 to Table 5-5). In terms of set point response, the maximum average rise time is 5.64 minutes for the age groups, and 14.31 minutes for the entire population. This reveals that the controllers satisfy the design specification on this parameter. However, the overshoot and settling time are not improved in comparison to those of the PID controller (compare Table 4-2 to Table 5-4). The performance of the controller for set point response was not the main objective of the design because it is assumed the anesthesiologists have full control of this phase. These drawbacks in the set point response could be avoided by designing a higher order pre-fdter. The controller designed for the entire population has an average overshoot for set point and disturbance of 0.71% and 2.12%, respectively (see Table 5-5). Only patients #3 and #16 have a set point overshoot greater than 10% (see Table C-7). Analyzing the complete simulation results shown in Table C-7 and Table C-8, these patients can be considered as outliers since their set point overshoot is more than twice the mean value of the entire population. If we remove the outliers, the data reveals that the remaining patients have a maximum overshoot of less than 3.14% for set point response and 7.02% for disturbance rejection. Thus, the controller satisfies the design specification for 42 out of 44 patients on the design for the entire population. Chapter 6 6 Controller Performance Evaluation This chapter describes the performance evaluation methodology used by the clinical community to test closed-loop controllers in anesthesia. This methodology has been used by Struys et al. [30] to compare the performance of PID and MPC controllers for the automatic administration of propofol using BIS as the control variable. The aim of this analysis was to study the controllers' responses to variations in patient types and interventions. The performance evaluation measures used for this analysis are explained in detail in [32]. In this thesis, we compare the performance of the PID and CRONE controllers explained in Chapter 4 and Chapter 5, respectively. The controllers are simulated for a 50% set point DOA, and a theoretical disturbance profile that emulates the patient's arousal reflexes during the surgical procedure. The content of this chapter is organized as follows. The performance evaluation measures, as published in [32], are explained in section 6.1. The simulation profile used to test the controllers is presented in section 6.2. Finally, the performance evaluation results are shown in section 6.3. 6.1 Performance Evaluation Measures In the anesthesia practice, the performance of a closed-loop system is assessed by calculating the bias, inaccuracy, divergence, and wobble measures. These four measures are based on the percentage performance error (PE), defined as the error vector between 63 Chapter 6: Controller Performance Evaluation 64 the system's output and the set point. This is calculated according to the following formula: PE„ = y's ~ y s p 100 (6.1) where yy is the/* sample of the system's output for patient;', and y is the target D O A . For our simulations, the set point is fixed to 50%. Bias (Median Performance Error) The bias measure describes whether the output value is systematically either above or below the set point. For the i'h patient, this performance measure is given by the median performance error as: MDPEi = median\PE9 J = \,...,N,\ (6.2) where Ni is the number of PE values for the i'h patient. The bias is expressed as a percentage [%]. Inaccuracy (Median Absolute Performance Error) The percent, [%], inaccuracy of the closed-loop system for each patient is calculated by the median absolute performance error: MDAPEj = median\pE.\,j = \,...,N,} (6.3) This performance measure indicates the expected size of the error between the system's output (y) and the set point (ysp). Divergence Divergence reflects the possible time-related trend of the output in relation to the set point. For each patient, it is calculated as the slope of the linear regression equation of Chapter 6: Controller Performance Evaluation 65 \PEij\ against time, and is expressed in units of percent divergence per minute [%/min]. It is calculated according to the following equation: where tu is the time (expressed in minutes) at which the corresponding PEy was determined. A positive value reveals that the output diverges from the set point, while a negative value reveals that the output converges to the set point. Wobble measures the total intra-patient variability of the performance error. This value is directly related to the controller's ability to achieve a stable D O A . For the fh patient, the percentage wobble is calculated as follows: This performance measure quantifies the oscillation of the output around the set point, and is expressed as a percentage [%]. 6.2 Simulation Experiment The predictive performance of the PID and CRONE controllers is assessed using the performance measures explained in the previous section. We have designed a standard experiment to test each of the 10 controllers (five PID's and five CRONE's) under the same conditions. The total experiment time is exactly one hour and ten minutes, and includes a virtual induction, maintenance, and emergence phase. DIVERGENCE, = (6.4) Wobb le WOBBLE\ = median\PEv - MDPEi I, j = 1,..., N,} (6.5) Chapter 6: Controller Performance Evaluation 66 The virtual induction phase is started 30 seconds after the beginning of the simulation. For each patient, a target WAV value of 50% is set as a step set point change and it is held constant until the end of the virtual surgery. Subsequently, one hour into the simulated experiment, the emergence phase begins by changing the set point to 100%; that is, stopping all drug administration to bring the patient back to a fully awake state. We also require a standard theoretical noxious stimuli profile to simulate the maintenance phase of anesthesia. This stimulus profile, or disturbance profile, emulates the patient's arousal reflexes during a surgical procedure. As mentioned in section 4.1, we assume that these disturbances are directly affecting the system's output (y) and can be simulated as an offset of the patient's D O A . The disturbance profile shown in Figure 6-1 was designed to emulate a typical noxious stimuli trajectory of a surgical procedure. Although the time-offset trajectory might seem to be too simple, we believe it emulates three typical stimuli. First, stimulus A simulates the arousal reflex due to the first surgical incision. Second, a decaying exponential shown as B in Figure 6-1 represents how the patient's offset slowly decreases but settles at an onset of 10% during which continuous normal surgical stimulations might occur. Finally, stimulus C simulates the withdrawal of stimulations during the skin-closing period. This stimulus profile is applied to all the patients and is fixed to start 20 minutes after the beginning of the experiment. Chapter 6: Controller Performance Evaluation 67 30 r _ 20 s 55 Vi 3 O 'S o Z 10 0 0 10 20 30 40 50 60 70 T i m e [minutes| Figure 6-1: Noxious Stimuli profile. 6.3 Controller Performance and Results Following the standard experiment explained in the previous section, all 176 virtual surgeries (i.e., 44 patients x 2 controllers x 2 for the age groups and the entire population) were simulated in Matlab Simulink. The simulation plots for each of the 10 controller are shown in Appendix C. The response of each controller is only analyzed for the maintenance phase of anesthesia, as we are assuming the anesthesiologist guides the induction and emergence phases of anesthesia. This assumption is consistent with how previous authors have tested the performance of other control algorithms. For example, Absalom et al. [1] and Struys et al. ([30], and [31]) published three studies where the controllers were run in open loop for the induction phase of anesthesia, and then after the patients have reached the desired Chapter 6: Controller Performance Evaluation 68 D O A the loop was closed and the controllers guided the drug administration regime until the end of the surgery. A summary of the performance results for the PID and CRONE controllers is presented in Table 6-1 and Table 6-2, respectively. Data are presented as mean ± standard deviation (SD) of the total number of patients in each design group. Table 6-1: PID's predictive performance summary. Group MDPE [%] PID's Perfon MDAPE [%1 ance Results Divergence [%/min] Wobble [%] mean SD mean SD mean SD mean SD 1 2 3 4 A l l 0.025 ± 0.357 -0.068 ± 0.128 0.122 ± 0.037 -0.064 ± 0.228 -0.023 ± 0.200 0.774 ± 0.625 0.485 ± 0.204 0.208 ± 0.084 0.629 ± 0.445 0.536 ± 0.452 -0.236 ± 0.067 -0.228 ± 0.040 -0.158 ± 0.034 -0.244 ± 0.045 -0.221 ± 0.059 0.791 ± 0.684 0.538 ± 0.232 0.269 ± 0.107 0.649 ± 0.448 0.571 ± 0.437 Table 6-2: CRONE's predictive performance summary. Group MDPE [%] CRONE's Perfc MDAPE [%] irmance Results Divergence [%/min] Wobble [%] mean SD mean SD mean SD mean SD 1 2 3 4 A l l 0.188 ± 0.403 0.094 ± 0.150 -0.198 ± 0.063 0.204 ± 0.261 0.233 ± 0.327 0.516 ± 0.756 0.295 ± 0.413 0.852 ± 0.248 0.587 ± 0.652 0.518 ± 0.706 -0.254 ± 0.046 -0.266 ± 0.028 -0.276 ± 0.019 -0.270 ± 0.036 -0.253 ± 0.044 0.408 ± 0.495 0.269 ± 0.276 0.681 ± 0.301 0.143 ± 0.233 0.352 ± 0.439 From the tables shown above, it is observed that the PID controllers have a lower mean MDPE value compared to those of the CRONE. As explained in section 6.1, this measure indicates the bias of the controller. The sign of the MDPE represents the Chapter 6: Controller Performance Evaluation 69 direction of the performance error (PE), where a negative value means that the controller tends to overdose the patient, and a positive value would indicate an under dosage of the patient. In our study, it can be seen from Table 6-1 and Table 6-2 that three out of five PID controllers tend to overdose, while only one CRONE controller exhibits the same behavior. Considering only the designs for the entire population and for Groups 1, 2, and 4, the CRONE controllers produced overall better MDAPE results compared to those of the PIDs. However, for Group 3 the PID produced a smaller value compared to the CRONE. The MDAPE value represents the size of the PE and thus a small value describes better performance of the controller in approaching the set point and eliminating control errors. In this regard, the deteriorated performance of the CRONE controller for Group 3 could be due to its higher set point settling time and overshoot as compared to those of the PID (compare Group 3 in Table 4-2 to Table 5-4, and Figure D-3 to Figure D-8). Both controller techniques produced negative divergence mean values, meaning that the system's output converges to the desired set point when using either the PID or CRONE controllers. The absolute value of this measure indicates the speed of convergence of the system. From the data shown in Table 6-1 and Table 6-2, it is identified that the CRONE controllers produce more negative values of divergence than the PID controllers do. These negative values for divergence indicate that the size of the PEs decrease more rapidly when the CRONE controllers were used. Comparing the wobble measures for both techniques, the data reveals that the CRONE controllers have much less overall oscillation when considering only the designs for the entire population, and for Groups 1, 2, and 4. However, as appreciated by the Chapter 6: Controller Performance Evaluation 70 divergence values for Group 3, there is more oscillation for the CRONE controller compared to the PID. This is the same deteriorated performance revealed by the MDAPE for Group 3 and that can be seen by comparing Figure D-3 to Figure D-8. Combining the information on MDPE, MDAPE, divergence, and wobble for the two control strategies, it can be concluded that the CRONE controllers' performance is better overall in four out of five design groups. That is, the CRONE has tighter control of anesthesia for the designs performed on the entire population and for Groups 1, 2, and 4. For Group 3, the performance of the PID controller is better in comparison to the CRONE. As mention earlier, this may only be due to the CRONE's deteriorated set point response. The most interesting results are the ones obtained from the designs on the entire population. The virtual surgery simulation plots for this group of patients using the PID and CRONE controllers are shown in Figure 6-2 and Figure 6-3, respectively. From these figures, it is seen that the CRONE controller provided tighter control of D O A with less oscillation around the set point and responded more rapidly to changes in the patient's D O A due to the noxious stimuli profile. In terms of the set point response, or the virtual induction phase of the experiment, the PID controller produced more oscillation and higher overshoots compared to those of the CRONE. As a result of this deteriorated performance shown in Figure 6-2, the infusion rate time course of the PID controller saturated more at the lower safety limit of zero (see section 4.4). In contrast, it is seen from Figure 6-3 that the infusion rate profile of the CRONE controller does not saturate at all during the virtual induction phase. Furthermore, these results reveal that the PID tends to use more volume of propofol overall. Chapter 6: Controller Performance Evaluation 71 100 80 E 6 0 > 40 20 0 _ 4 10 20 30 40 Time |minutes] 30 40 Time |minutes] 50 50 Patients Y sp Disturbance 60 70 60 70 Figure 6-2: PID performance experiment simulation (Entire population). 100 n 80 60 %l > 40 -20 0 0 10 20 30 40 Time |minutes| — Patients . . Y sp '"• Disturbance 50 60 70 _ 4 30 40 T i m e Intimites] Figure 6-3: CRONE performance experiment simulation (Entire population). Chapter 6: Controller Performance Evaluation 72 In terms of the controllers' performance during the virtual maintenance phase of the experiment for the entire population, the output disturbance rejection of the CRONE is better in comparison to that of the PID (compare Figure 6-2 to Figure 6-3). The response of the CRONE controller was shown to have less oscillation and overshoot. Furthermore, during this phase of the virtual surgery the PID saturated the infusion pump at both the low and high limits, while the CRONE did not saturate the actuator at all. This drawback showed to be one limitation of the PID controller design. From the performance comparison presented in this chapter, and the time domain responses of the PID and CRONE controllers explained in Chapters 4 and 5, it is seen that the CRONE controller performs best overall. This suggests that a higher order controller can compensate more appropriately for the inherent variations in the patients' drug-response characteristics. As mentioned earlier, the most interesting result is that both controllers are able to guide the drug administration regime of propofol for the designs on the entire population. This clearly suggests that good closed-loop performance can be achieved by identifying first a nominal PKPD(s) model that minimizes the uncertainty within a specific operating frequency. Chapter 7 7 Conclusions 7.1 Summary and Contributions of the Thesis The objectives of this thesis were to design and evaluate the performance of two closed-loop controllers for the hypnotic state of anesthesia using WAV as the control variable. The aim of the controller designs was to provide an adequate drug regime for the administration of propofol to avoid over or under dosing the patient. The proposed controller designs were based on the following four main objectives: 1) The controllers were designed based on a linearized approximation of the patient's drug-response characteristics. These models were obtained by combining the three compartment pharmacokinetic models of propofol published by Schiittler and Ihmsen [25], and the clinically identified pharmacodynamic models using WA F a s a measure of hypnosis presented in a doctoral thesis by Bibian [8]. 2) The controllers were designed to be robust with respect to the model's uncertainty to compensate for the inherent inter-patient drug-response characteristics. 3) The controllers' design focused on minimizing the closed-loop output disturbance response due to typical noxious stimulations of a surgical procedure. 4) Finally, a pre-filter compensator was also designed to improve the set point response of the system. The major difficulty in the design of automatic controllers for anesthesia is the inherent inter-patient drug-response variability. This variability results from the extreme 73 Chapter 7: Conclusions 74 uncertainty associated with the pharmacokinetics and pharmacodynamic of the patients. In general, these PKPD(s) differences are attributed to differences in age, lean body mass, previous diseases, drug tolerance and sensitivity, and drug addictions. Nevertheless, anesthesiologists are able to compensate for these inter-individual variability's by adapting the drug administration regime to the individual patient's needs, and therefore guaranteeing an adequate D O A across the population. In order to reduce this uncertainty, this thesis contributes with a systematical procedure to identify a nominal PKPD(s) model that captures more adequately the dose-response characteristics of the five design groups. Although it is impossible to eliminate the patients' variability, these identified nominal models minimized the uncertainty circles and were useful in the design of robust controllers. In this thesis, we have proposed the design of two robust controllers. The first proposed controller is fixed to a typical PID structure. This control structure was chosen for two reasons. First, it was found in literature that a PID controller was designed and evaluated for the administration of propofol using BIS as the control variable [1]. This controller produced continuous oscillations around the set point when it was tested in three out of ten patients undergoing major orthopedic surgery. We believe this deteriorated performance was due to its poor robustness and thus decided to revise its design procedure using a control engineering approach. Second, we wanted to investigate i f a simple PID controller is suitable for the automatic administration of propofol. To compare the performance of a low order to a higher order controller, we decided to investigate the CRONE controller design technique. The CRONE controller is based on a fractional order of integration and allowed us to synthesize a robust high order controller. Chapter 7: Conclusion 75 Based on the two proposed control algorithms, this thesis contributes with a systematic procedure for the design of a PID and CRONE controller following the design specifications given in Chapter 4. The design procedures for each controller are divided in three parts. First, the main design criterion and optimization procedure was presented for each strategy with the objectives of designing a controller with robustness to models uncertainty, and with good output disturbance rejection. Second, to improve the set point response of the closed-loop system a two degree of freedom structure was implemented by explaining the design of a pre-filter compensator. Finally, for each proposed controller we have explained an anti-windup scheme to compensate for the actuator saturation. The systematic design procedures contributed by this thesis can be easily changed to design controllers with different specifications. The design variable in each procedure can be updated to change the performance of the closed-loop system; allowing the engineer, or anesthesiologist, to trade off between performance and robustness. The performance of the PID and CRONE controllers was tested with a standard simulations experiment to emulate the arousal reflexes of the patients during a surgical procedure. Overall, the higher order of the CRONE controller was shown to give better output disturbance rejection than that of PID controller. 7.2 Future Work Some of the possible research directions related to the work presented in this thesis are outlined next. In reality, anesthesiologists assess the patient's D O A by monitoring several physiological signals such as blood pressure, heart rate, breathing rhythm, body temperature, body fluids, and by interpreting the E E G or its single value index for Chapter 7: Conclusion 76 hypnosis. From the evaluation of all these signals and the experience of the anesthesiologist, he or she administers different drugs (e.g., anesthetics, opioids, or neuromuscular agents) to influence one of the functional components of anesthesia (e.g., hypnosis, analgesia, and muscle relaxation). Thus, from this point of view, the automatic control of anesthesia can be treated as a M I M O system where the ideal control system would steer multiple drugs simultaneously to affect all the components of anesthesia. However, this multivariable approach is still not applicable because there are no models that can accurately represent the interaction between drugs and their effects. In terms of SISO control of anesthesia, a master-slave configuration can be implemented. For instance, this configuration can be adopted to have the master controller regulate the hypnotic state of the patient, while the slave controller regulates the drug plasma concentration. Furthermore, a fuzzy logic controller can also be implemented to take into account certain safety restrictions from the standard practice of anesthesia. Finally, yet importantly, an artifact-tolerant structure can be designed to compensate for corrupted measurement signals. In the operating room, the feedback signal can be corrupted by sudden disconnections of the sampling lines, routine calibrations of the sensing devices, or from the high impedance of the electrodes. Although the proposed controllers compensate for the actuator saturation, an artifact protection scheme can be implemented to guarantee the patients' safety when the control variable is compromised. Bibliography [1] A . Absalom, N . Sutcliffe, and G. Kenny, "Closed-loop Control of Anesthesia Using Bispectral Index: Performance Assessment in Patients Undergoing Major Orthopedic Surgery under Combined General and Regional Anesthesia", Anesthesiology, vol. 96, pp. 67-73, Jan. 2002. [2] K . Astrom and T. Hagglung, PID Controllers: Theory, Design, and Tuning. Research Triangle Park, N . C , Instrument Society of America, 1995. [3] S. Bibian, G. Dumont, M . Huzmezan, and C. Ries, "Quantifying Uncertainty Bounds in Anesthesia PKPD Models", in Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biological Science, vol. 1, pp. 786-789, 2004. [4] S. Bibian, C. Ries, M . Huzmezan, and G. Dumont, "Clinical Anesthesia and Control Engineering: Terminology, Concepts and Issues", European Control Conference, Cambridge, U K , Sept. 1-4, 2003. [5] S. Bibian, T. Zikov, G. Dumont, C. Ries, E. Puil, H . Ahmadi, M . Huzmezan, and B. MacLeod, "Estimation of anesthetic depth using wavelet analysis of electroencephalogram", in Proceeding of the 23rd annual EMBS International Conference, Istanbul, Turkey, Oct. 25-28, 2001. [6] S. Bibian, T. Zikov, G. Dumont, C. Ries, E. Puil, H . Ahmadi, M . Huzmezan, and B. MacLeod, "Method and apparatus for the estimation of the anesthetic depth using wavelets analysis of the electroencephalogram", Patent Pending, #60/395,313, January 15 t h 2004. 77 Bibliography 78 [7] S. Bibian, T. Zikov, C. Ries, G. Dumont, and M , Huzmezan, "Intra- and Inter-patient Variability During Induction of Anesthesia", IEEE Conference on Engineering in Medicine and Biology, 2004. [8] S. Bibian, "Automation in Clinical Anesthesia", PhD. Thesis, The University of British Columbia, work in progress. [9] R. Dripps, J. Eckenhoff, and L. Vandam, Introduction to Anesthesia: The Principle of Safe Practice. Philadelphia, Pennsylvania: W. B . Sauders Company, seventh ed., 1988. [10] R. Dripps, J. Eckenhoff, L . Vandam, D. Longnecker, and F. Murphy, Introduction to Anesthesia. Philadelphia, Pennsylvania: W. B . Sauders Company, eighth ed., 1992. [11] C. Frei, A . Gentilini, M . Derighetti, A . Glattfelder, M . Morari, T. Schnider, and A . Zbinden, "Automation in Anesthesia", in Proceedings of the American Control Conference, San Diego Cai., USA, June 1999. [12] A . Gentilini, C. Frei, A . Glattfedler, M . Morari, T. Sieber, R. Wymann, T. Schnider, and A . Zbinden, "Multitasked Closed-Loop Control in Anesthesia: Automatic Controllers Capable of Regulating Multiple Patient Outputs for High-Quality Anesthesia Treatment", IEEE Engineering in medicine and Biology Magazine, vol. 20, no. 1, pp. 39-53, Jan.-Feb. 2001. [13] A . Gentilini, M . Rossoni-Gerosa, C. Frei, R. Wymann, M . Morari, A . Zbinden, and T. Schnider, "Modeling and Closed-Loop Control of Hypnosis by Means of Bispectral Index (BIS) with Isoflurane", IEEE Transactions on Biomedical Engineering, vol. 48, no. 8, pp. 874-889, Aug. 2001. Bibliography 79 [14] G. Goodwin, S. Graebe, and M . Salgado, Control System Design, Upper Saddle River, New Jersey, Prentice Hall, 2001. [15] Graseby Medical. Graseby 3400 Product Information. Cited 2005 http://www.graseby.co.uk/product frame.html. [16] P. Glass, and I. Rampil, "Automated Anesthesia: Fact or Fantasy?", Anesthesiology, vol. 95, pp. 1-2, Jul. 2001. [17] J. Huang, Y . Lu , A . Nayak, and R. Roy, "Depth of Anesthesia Estimation and Control", IEEE Transactions on Biomedical Engineering, vol. 46, no. 1, pp. 71-81, Jan 1999. [18] Information for Health Professionals. Data Sheet: Propofol. Cited 2005. http://www.medsafe.govt.nz/profs/Datasheet/p/Propofol 10m gpermlinj.htm. [19] E . Mortier, M . Struys, T. De Smet, L . Versichelen, and G. Roily, "Closed-loop controlled administration of propofol using bispectral analysis", Anesthesiology, vol. 53, pp. 749-754, Jan. 1998. [20] E. Mortier, and M . Struys, "Monitoring the depth of anesthesia using bispectral analysis and closed-loop controlled administration of propofol", Best Practice & Research Clinical Anesthesiology, vol. 15, no. 1, pp. 83-96, 2001. [21] A . Nayak, and R. Roy, "Anesthesia Control Using Midlatency Auditory Evoked Potentials", IEEE Transactions on Biomedical Engineering, vol. 45, no. 4, pp. 409-421, Apr. 1998. [22] A . Oustaloup and B . Mathieu, La commande CRONE: du escalaire au multivariable, H E R M E S , Paris, 1999. Bibliography 80 [23] H . Panagopoulos, "PID Control: Design, Extension, and Application", PhD Thesis, Lund Institute of Technology, 2000. [24] T. Schnider, C. Minto, P. Gambus, C. Andresen, D. Goodale, S. Shafer, and E . Youngs, " The Influence of Method of Administration and Covariates on the Pharmacokinetics of Propofol in Adult Volunteers", Anesthesiology, vol. 88, pp. 1170-1182, May 1998. [25] J. Schuttler, and H . Ihmsen, "Population Pharmacokinetics of Propofol: A Multicenter Study", Anesthesiology, vol. 92, pp. 727-738, Mar. 2000. [26] H . Schwilden, and H . Stoeckel Eds., Control and Automation in Anaesthesia, Berlin, Springer, 1995. [27] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Analysis and Design. John Wiley & Sons, 1 ed., July 1996. [28] D. Stanki, and W. Watkins, Drug Disposition in Anesthesia, New York, New York USA, Gruner & Stratton, 1982. [29] R. Stoelting, and R. Miller, Basics of Anesthesia. Philadelphia, Pennsylvania: Churchill Livingstone, fourth ed., 2000. [30] M . Struys, T. De Smet, S. Greenwald, A . Absalom, S. Binge, and E. Mortier, "Performance Evaluation of Two Published Closed-loop Control Systems Using Bispectral Index Monitoring: A Simulation Study", Anesthesiology, vol. 100, pp. 640-647, Mar. 2004. [31] M . Struys, T. De Smet, L . Versichelen, S. Van de Valde, R. Van den Broecke, and E. Mortier, "Comparison of Closed-loop Controlled Administration of Propofol Bibliography 81 Using Bispectral Index as the Controlled Variable versus "Standard Practice" Controlled Administration", Anesthesiology, vol. 95, pp. 6-17, Jul. 2001. [32] J. Varvel, D. Donoho, and S. Shafer, "Measuring the Predictive Performance of Computer-Controlled Infusion Pumps, Journal of Pharmacokinetics and Biopharmaceutics, vol. 20, no. 1, pp. 63-94, 1992. [33] R. Wiklund, and S. Rosenbaum, "Anesthesiology", The New England Journal of Medicine, vol. 337, no. 16, pp. 1132-1141, Oct 1997. Appendix A: Laryngeal Mask Airway Study Population Table A - l : L M A Population and PD parameters (Group 1 & 2). As presented in [8]. Patient # Demographics PD Parameters ASA Weight Height Age Gender Td Kd EC so y [Kg] [cm] [yr] [s] rs"'-10"3l [pg/mL] 01 n.c. 100 178 21 Male 22 133.5 3.2 A.l 02 n.c. 59 168 28 Female 4 AAA 3.1 2.5 03 1 90 190 26 Male 44 25.0 2.4 1.9 o 04 1 53 157 21 Female 45 51.5 3.8 1.2 ? 05 1 90 185 19 Male 39 85.7 3.8 2.3 06 1 60 162 28 Female 18 82.5 3.9 2.1 » ft oc 1—1 07 1 78 170 24 Male 32 AAA 2.9 2.8 08 1 68 164 19 Female 12 26.7 1.9 2.3 09 1 70 170 25 Male 7 35.2 3.4 1.9 n. s 10 1 81 180 23 Male 9 32.8 2.8 2.8 o i -O 11 1 83 178 18 Male 17 46.4 2.8 2.3 12 1 67 163 21 Female 4 26.2 2.4 2.5 13 1 88 183 22 Male 9 50.4 2.5 2.6 14 1 59 162 21 Male 18 160.5 3.6 3.9 15 1 72 176 19 Male 20 75.0 2.5 1.9 16 n.c. 73 180 34 Male 44 54.8 3.2 2.7 17 n.c. 71 178 34 Male 29 83.1 4.0 2.3 © 18 1 91 189 35 Male 18 34.4 3.7 2.1 V ce i . 19 1 85 180 38 Male 1 29.6 3.3 1.2 M CU ^ 20 1 86 168 33 Female 1 24.9 3.1 1.5 ft O 21 1 75 176 38 Male 12 35.2 3.9 1.8 22 1 91 178 39 Male 4 24.8 2.7 2.0 f N O. 23 1 78 178 34 Male 12 28.7 2.8 2.2 s o 24 1 65 163 36 Female 5 27.0 2.8 2.1 O 25 1 58 157 34 Female 4 67.2 3.6 2.0 26 1 77 180 33 Male 12 29.3 3.1 1.8 27 1 105 185 39 Male 13 29.1 3.7 1.8 82 Appendix A : Laryngeal Mask Airway Study Population 83 Table A-2: L M A Population and PD parameters (Group 3 & 4). As presented in [8]. Patient # Demographics PD Parameters ASA Weight Height Age Gender Td Kd EC s o 7 [Kg] [cm] [yr] [s] [s _ ,-10 3] [pg/mL] o m 28 2 66 168 46 Female 11 36.6 6.1 1.3 V 29 1 99 182 45 Male 2 32.6 4.7 1.3 « 30 1 77 173 48 Male 12 35.0 4.5 1.4 o 31 98 189 47 Male 10 28.7 3.9 2.0 32 1 79 172 46 Male 12 34.8 3.9 1.8 <*•> a. 33 1 96 193 40 Male 9 36.6 3.2 1.9 s o 34 1 77 176 40 Male 8 35.6 3.4 2.3 O 35 1 97 190 48 Male 13 30.0 3.0 1.9 o 36 n.c. 80 177 53 Male 35.0 38.0 3.3 1.8 VI 37 1 63 165 60 Female 3 31.5 3.5 2.3 I- 38 1 91 185 52 Male 29 42.0 4.4 2.2 v? 39 100 182 52 Male 2 21.8 4.7 1.4 V 1 O ITi 40 1 97 176 52 Male 16 28.8 3.7 1.1 T t 41 1 95 183 50 Male 10 26.4 4.0 1.9 a. s 42 2 56 173 52 Female 6 58.0 5.0 1.5 o k. O 43 2 100 180 54 Male 6 32.2 4.2 1.5 44 3 95 176 60 Male 12 24.3 3.1 1.8 Appendix B: Numerical Solution Procedure B.l: PID Controller In Chapter 4, it was shown that robustness to process variations can be expressed by the maximum sensitivity Ms. This condition specified that the Nyquist curve of the open loop transfer function should avoid the circles enclosing the critical point -1. As explained in section 4.2, this robustness constraint is defined by the following equation: f{kp,ki,kd,co)=\ + \kp-^-(ki-cv2-kd)-G0(j-co) >R2 (B. l ) where G0(joj) is the frequency response of the nominal model, R is the inverse of maximum sensitivity (R = \/Ms ), and kp, kt, and kd are the PID controller parameters. Let the frequency response of the nominal model be expressed by: GQ {j •co) = a{co) + j- p{co) = p{co) • eJ' j<p(a>) (B.2) Thus, by substituting equation (B.2) in (B. l) the robustness constraint can be written as [23]: { R J kP+'(\2 V >i (B.3) For fixed kd and co, the equation above represents the exterior of an ellipse in the kp-k\ plane, as shown in Figure B-1. The ellipse has its center located at kp = - a(a>)/p(w)2 and kt =a>2 -kd -»•/?(&>)/p(tx>)2, and has its axes are parallel to the coordinate axes. 84 Appendix B : Numerical Solution Procedure 85 The horizontal half axis has a length of Rjp(co) and the vertical half axis has a length of Figure B - l : Graphical illustration of robustness constraint in Equation (B.3). When the frequency co ranges from 0 to co, the ellipses generate an envelope which defines one boundary of the sensitivity constraint. Since the function / is quadratic in kt the envelope has two branches, where only one branch corresponds to stable closed-loop systems. To have a stable closed-loop system it is required that the integral gain is positive. The largest value of k, for fixed kd on the envelope occurs at the maximum point on the locus of the lower vertex of the ellipse (see Figure B - l ) . From Equation (B.3), the locus of the lower vertical vertex is given by: co •R/p{co). • kj a(co) p{cof (B.4) ki(co)= co2 -kd co o • B(co) co • R p(cof p{co) Appendix B: Numerical Solution Procedure 86_ The largest value of the integral gain in the envelope can thus be found by maximizing ki on the locus of the lower vertex. Differentiating the equation for ki in (B.4) gives the condition for extremum defined by: h„D{a>) = {R + sm{<p(co)))(^\ " -] " <M- M<M) + 2 -R-kd=0 (B.5) p{co) CO where the dot denotes the derivative with respect to co. Therefore, the solution to equation (B.3) is reduced to a single algebraic equation (B.5) in co. Solving this equation gives the frequency coa for which we can compute the controller parameters kp and k( given by equation (B.4). To ensure that the solution corresponds to a maximum, it is also required that the following condition be satisfied: 7 T W > 0 (B.6) dco Equation (B.5) is solved with the Optimization Toolbox in Matlab 6.5. Note, however, that there may be several solutions for this equation i f the optimization routine is started from different initial conditions. To alleviate this problem, the optimization procedure is divided in two steps. First, a PI controller is designed by fixing kd = 0 in equation (B.5) and finding its roots within the following frequency interval: <y90 < co < co^^w (B.7) where co^ is the frequency where the nominal model has a phase of <p. Second, once the PI controller has been design successfully, the design procedure is performed again with Appendix B: Numerical Solution Procedure 87 a kd value different than zero. For the PID design, equation (B.5) is solved with the constraints in equation (4.9), as explained in section 4.2. To summarize, the PID design problem can be solved using the following procedure: 1) Choose the design parameter Ms and compute R = \/Ms . 2) From the nominal model of the plant, determine the search interval in equation (B.7) to design the PI controller, and determine the number of PID constraints in equation (4.9) that apply to the system under control. As explained in section 4.2, these PID constraints must be treated differently for system with integral or close to integral action. In this thesis, the second constraint which specifies that the Nyquist curve of L(j • co) should have negative curvature was omitted because the drug-response models had poles close to zero. 3) Use the Optimization Toolbox in Matlab 6.5 to solve equation (B.5) with kd = 0 and evaluate condition (B.6). If this Step fails, change the design parameter Ms, recomputed R, and return to Step 2. 4) Once the PI design has been successful, use the found kp and kt parameters to define a new PID search interval given by: co, slarl <co<co stop CO, o CO, slarl 2 {com+CO210) 2 (B.8) CO, stop Appendix B: Numerical Solution Procedure 88_ where dJ0 is the frequency at which the sensitivity function from the PI design has its maximum. a>m and dJ270 are the frequencies where the argument of the open loop transfer function from the PI design is -180° and -270°, respectively. 5) Fix the derivative gain kd to a value different than zero and solve equation (B.5) within the frequency interval given in equation (B.8). 6) Compute the PID constraints in equation (4.9) and test i f they are satisfied. 7) Increase the value of kd until the robustness (B.3) or PID (4.9) constraints are violated. A detailed explanation of this design procedure can be found in [23]. B.2: CRONE Controller In Chapter 5, it was shown that robust stability with respect to plant uncertainty can be achieved by correctly positioning the nominal open loop transfer function in the Nichols diagram. This nominal open loop transfer function takes into account the control specifications at low and high frequencies, the generalized template around unit gain frequency cou, and the non-minimum phase part of the PKPD(s) models (see section 5.2). The optimization of the generalized template consists of determining the seven optimal parameters of the transfer function B(s) (expressed by equations (5.4) and (5.5)) that define the following characteristics: • Its position along the 0 dB axis given by the optimal real order of integration, a, and gain K. Appendix B: Numerical Solution Procedure 89 • Its direction or inclination in relation to the vertical given by the optimal imaginary order of integration b. • Its length around the ultimate frequency CDU given by optimal transitional frequencies coh and coh, and • Its prolongation in low and high frequencies given by the optimal extreme transitional frequencies co'h and co'h. However, the optimization procedure is reduced to only five independent parameters by performing the variable transformation [22]: b - , COu CO, a. cok cor COu (B.9) CO. coL where cor is the resonance frequency of the closed-loop system. From this transformation, the frequency response of the open loop transfer function defined by equation (5.4) and (5.5) can be now written as: B^-cor)^Gnmp{j-co\K-B[{j.cor\Bh{j (B.10) where G (j • cor) is the frequency response of the non-minimum phase part of the nominal model evaluated at cor and Appendix B: Numerical Solution Procedure 90 f , A " ' " 1 + cosl b • ln 1 + j • a! (B .H) i For known values of a and the point M 0 ( X 0 , 7 n ) at which the template tangents to the overshoot contour TQ , the calculations are reduced to only finding the parameters K, a, b and cor. Note here that X0 [deg] and Y0 [dB] are solely defined by the desired magnitude of the closed-loop resonance peak Qd [dB]. For our design, we have set Qd equal to le" 4 [dB]. Following this approach, the first parameter to be determined is the gain K and is calculated from the gain placement Y0 of the template at frequency cor as: K = 10 r„ /20 (B.12) The non-minimum phase part of the plant is ignore here and will be added at a later Step. The second parameter to be determined is the real order of integration, a, and is calculated from the phase placement Xn of the template at frequency a>r as: X0 = argP'b{j-cor) + argJ3b(j cor)+ a r g - c o r ) + arg/?,(j• cor) + arg(3'h(j• cor) (B. 13) where the phase expressions are defined by: Appendix B: Numerical Solution Procedure 91_ ngB'h{j.cor) = nb\-9V + 0'b) ngBb{j-cor) = 0b-0'b avgB{j.cor) = a{eh-0h) (B.14) &rgj3h{j.cor) = 2.(0'h-0h) argB'h{j•cor) = -nh -0h and 0 = arctan(or) (B.15) Thus, substituting the equations in (B.14) into (B.13) the real order of integration is expressed by: a = —^—-{Xo-nb.{-9O° + 0'b)-{0h-0'b)-2.{0'h-0h) + nh-0h) (B.16) The third parameter to be determined is the imaginary order of integration, b, and is calculated from the slope of the generalize template that tangents to the overshoot contour YQ at the point M0. This is done by finding the value of b that satisfies the following equation: n-d2 + 9 1n(10)• At> of, where d,= "* a l ' 1 + a2h \ + al d2—^ —^-(B.18) \ + al„ \ + alb and Appendix B: Numerical Solution Procedure 92 9 1n(l0) cos\—-X 0 +10- K » / 2 0 (B . l 9) Equations (B.l2), (B.l6), and (B.l7) calculate the parameters K, a, and b that values are calculated from the resonance frequency cor and the desired magnitude of the closed-loop resonance peak Qd set by the designer. Finally, the last parameter to be determined is the resonance frequency, cor, and is calculated to guarantee that the open loop transfer function has a gain of 0 dB at the ultimate frequency cou. Given that K, a, and b are already known at this point, the resonance frequency cor can be calculated iteratively to assure the next equation is satisfied: P0{j-co^ = \K\\p[{j.a\\BXj-co^ (B.20) To summarize, the CRONE optimization design problem can be solved using the following procedure: 1) Choose the desired ultimate frequency cou and the magnitude of the closed-loop resonance peak Qd. 2) From the nominal model of the plant, define the asymptotic behavior orders m and rih at low and high frequencies, respectively. 3) Extract the frequency response envelope from all the parametric states of the plant (see section 5.2.1). characterize the open loop frequency response of /?„(/•&>) in equation (5.4). These Appendix B: Numerical Solution Procedure 93 4) Define the optimization vector v = [Y0 a'h ah ah a'h J and choose suitable initial conditions. In out implementation we found the best solution to the optimization problem by setting the initial guess to v0 = [l.5 10 10 0.1 O.l]; which corresponds to a generalized vertical template that is tangent to the overshoot contour at 1.5 dB and has a length of two decades. 5) Solve the optimization problem using the Optimization Toolbox in Matlab 6.5. i . Calculate K, a, and b from equations (B.12), (B.16), and (B.17), respectively. i i . Calculate cor to guarantees that the open loop transfer function has a gain of 0 dB at cou. i i i . Construct the frequency response of the open loop transfer function in equation (B.10) by adding the response of the non-minimum phase part of the plant. iv. Find the parameter vector that minimizes the objective function in equation (5.7) under the constraints given in Table 5-1. 6) Once the optimum open loop template has been found, identify the rational CRONE controller as explained in section 5.3. A detailed explanation of this design procedure can be found in [22]. Appendix C: Set Point and Disturbance Response Analysis C l : PID Response Data Table C - l : Complete PID set point & output disturbance response data (Group 1 & 2). I] H Performance Measurements Patient Set Point Disturbance # Tr T OS T OS IAE [min] [min] [%] [mini [%] [-1 01 4.56 8.20 6.26 6.10 4.19 15.25 02 8.72 6.36 1.15 7.38 3.88 17.53 03 3.99 13.49 13.86 11.47 11.95 38.34 04 6.92 10.08 2.12 16.85 7.20 46.07 7 05 5.17 10.73 9.21 11.16 7.58 33.63 k. 06 6.73 6.04 1.63 10.22 4.65 24.91 « ft oo "H 07 6.25 3.88 1.30 6.72 7.19 22.60 08 10.71 8.15 0.65 6.20 6.42 18.04 i-H 09 6.49 5.92 1.90 11.08 5.39 28.02 C s Q 10 8.85 4.89 1.29 7.31 5.45 19.75 11 5.51 7.72 2.41 9.05 6.14 24.64 12 10.32 6.23 0.90 6.78 5.59 18.28 13 6.57 5.57 1.19 7.71 4.62 18.98 14 5.50 4.95 1.34 5.97 4.02 13.97 15 7.58 6.17 1.04 9.32 5.35 23.66 16 2.79 27.19 36.05 5.68 16.93 24.50 17 7.96 4.51 1.12 6.44 5.36 19.38 o 18 4.38 4.06 1.57 7.47 6.59 23.48 V SB 19 6.64 5.80 1.20 12.29 5.00 30.65 cs 20 6.00 5.13 0.99 9.33 5.85 25.88 21 6.13 4.95 1.16 8.24 5.49 23.59 22 8.89 5.98 0.83 6.44 5.62 18.78 r i a. 23 8.77 6.75 0.74 5.62 6.00 17.41 s o s-24 9.05 6.72 0.65 5.48 4.98 15.67 25 7.97 5.80 0.86 6.63 3.44 16.26 26 8.88 5.17 0.90 7.24 6.26 21.83 27 4.74 6.78 2.59 9.22 6.43 27.69 94 Appendix C: Set Point and Disturbance Response Analysis 95 Table C-2: Complete PID set point & output disturbance response data (Group 3 & 4). Patient # Performance Measurements Set Point Tr Ts OS fminl [min] [%] Disturbance Ts OS IAE [min] [%] [-] o IT) V S B « ft O c s © O 28 29 30 31 32 33 34 35 6.37 6.23 7.03 7.24 7.49 7.60 7.88 7.88 5.24 5.20 5.39 5.37 5.49 5.48 5.58 5.60 0.88 0.90 0.65 0.62 0.53 0.49 0.43 0.44 7.86 8.14 5.96 3.87 1.87 1.89 1.78 1.60 3.23 3.18 2.77 2.27 1.75 1.56 0.78 1.41 22.34 22.31 18.31 14.18 13.61 12.70 10.44 12.07 o VI «x> k. cs ft o m a. a o i-a 36 37 38 39 40 41 42 43 44 2.84 8.94 5.51 5.07 5.95 4.36 6.93 5.50 8.80 5.85 6.48 3.26 8.39 5.25 4.05 5.56 4.95 6.88 3.14 0.64 1.24 3.55 1.04 1.56 0.97 1.33 0.79 5.84 4.87 6.46 11.72 12.61 7.48 8.73 9.98 6.45 8.99 3.62 7.38 6.29 6.04 6.65 4.06 5.37 6.61 24.16 13.21 22.81 33.38 34.87 23.39 22.27 27.61 20.73 Appendix C: Set Point and Disturbance Response Analysis 96 Table C-3: Complete PID set point & output disturbance response data (Entire populations, patients 1-27). I H Performance Measurements Patient Set Point Disturbance # Tr T 1 s OS T OS IAE [mini [min] [%] [min] [%] [-] 01 3.59 19.28 12.95 3.97 7.52 12.43 02 9.62 6.77 0.57 4.97 3.32 12.53 03 2.99 24.87 30.41 13.09 17.04 42.62 04 9.24 7.10 0.82 10.22 9.12 36.60 05 3.44 6.02 4.88 7.11 8.66 26.02 06 8.98 5.68 0.81 6.53 4.45 17.51 07 2.71 7.55 17.49 4.24 10.31 18.80 08 10.66 7.33 0.33 4.29 5.99 13.53 09 9.41 5.29 0.79 7.15 5.53 19.95 10 9.37 7.09 0.63 4.97 5.02 14.18 11 9.52 4.79 0.68 6.03 6.04 17.81 12 9.97 7.31 0.44 4.73 5.11 13.34 13 9.58 6.53 0.60 5.23 4.01 13.49 14 4.10 28.01 5.32 4.04 6.70 10.93 15 9.83 6.52 0.53 6.17 5.09 16.77 16 2.72 33.18 37.07 5.35 19.90 26.15 • 17 8.82 5.70 0.93 6.13 5.14 18.21 18 9.09 4.31 1.08 7.08 6.47 22.22 19 7.24 6.01 0.89 11.40 5.03 29.03 20 9.37 5.45 0.86 8.79 5.84 24.55 21 8.78 5.37 0.98 7.77 5.43 22.35 22 9.38 6.99 0.68 6.16 5.43 17.81 23 9.37 7.11 0.61 5.40 5.65 16.44 24 9.70 7.12 0.54 5.27 4.68 14.84 25 9.05 6.26 0.72 6.32 3.33 15.52 26 9.36 6.64 0.73 6.88 6.13 20.69 27 8.42 4.57 1.25 8.66 6.44 26.26 Appendix C: Set Point and Disturbance Response Analysis Table C-4: Complete PID set point & output disturbance response data (Entire populations, patients 28-44). 97 Performance Measurements Patient Set Point Disturbance # Tr T OS T OS IAE [mini [min] [%] [min] [%] [-1 28 6.37 10.03 2.20 14.37 5.07 35.83 29 6.40 10.81 2.45 14.60 4.97 35.58 30 6.68 5.69 1.19 11.00 5.11 29.43 31 8.76 4.65 1.22 7.72 5.70 23.00 32 8.59 5.62 1.01 7.72 5.09 22.07 a 33 8.82 5.60 0.92 7.26 4.92 20.24 mlatio 34 9.05 6.39 0.80 5.81 4.31 16.13 mlatio 35 9.05 6.59 0.83 6.64 5.45 19.94 mm o 36 8.91 5.96 0.78 6.16 7.72 23.46 M N u JM 37 9.25 6.84 0.70 5.01 3.15 13.61 8 38 8.46 4.47 1.29 6.91 6.38 23.19 39 5.75 9.07 2.62 12.98 5.70 34.90 40 6.91 5.95 1.11 14.05 5.45 36.41 41 5.76 4.67 1.24 8.03 5.89 24.17 42 7.60 6.19 1.09 9.54 3.69 23.36 43 6.38 5.59 1.31 11.00 4.84 28.86 44 9.03 6.90 0.84 6.83 5.81 21.13 Appendix C: Set Point and Disturbance Response Analysis 98_ C.2: CRONE Response Data. Table C-5: Complete CRONE set point & output disturbance response data (Group 1 & 2). H H Performance Measurements Patient Set Point Disturbance # Tr T 1 s OS T 1 s OS IAE [min] [min] [%] [min] [%] [-1 01 5.04 9.74 6.46 2.52 1.71 15.80 02 6.40 5.56 0.96 3.58 0.77 17.81 03 4.38 11.00 13.63 9.51 3.56 32.20 04 7.44 6.61 1.86 9.48 0.33 46.05 05 6.41 12.57 5.45 5.60 1.32 33.36 VI fc. a 06 7.43 6.52 1.16 5.45 0.44 25.67 18<ye; 07 4.32 7.98 6.82 6.64 2.59 21.98 18<ye; 08 4.00 6.26 3.46 6.61 3.26 19.22 09 7.05 6.34 1.71 5.83 0.58 28.60 a s 10 5.02 8.72 3.88 3.24 1.73 20.09 Gro 11 5.61 10.34 4.11 4.24 1.61 25.34 Gro 12 4.59 7.47 3.33 6.32 2.25 19.06 13 5.65 8.63 2.45 3.57 1.26 19.55 14 5.14 8.56 3.21 2.60 1.35 14.84 15 6.14 8.58 2.06 4.49 1.18 24.46 16 3.16 10.17 32.03 5.80 8.80 24.86 17 3.97 7.76 5.21 3.16 1.80 21.15 o •<* 18 4.27 8.62 6.99 3.85 1.76 24.82 V VI fc. 19 6.42 5.63 1.45 7.88 0.11 34.69 CQ W 20 4.92 8.89 3.61 5.16 0.94 27.76 ft o 21 4.76 8.30 3.39 4.53 0.91 25.02 22 3.81 7.14 5.63 6.31 2.15 20.37 a. 23 3.37 6.56 8.42 6.47 3.21 19.70 a © 24 3.42 6.06 5.02 5.90 2.69 17.66 O 25 5.01 4.49 1.70 3.59 0.86 17.86 26 3.98 7.75 5.58 3.63 1.95 23.47 27 4.98 9.46 4.73 5.21 0.85 28.89 Appendix C: Set Point and Disturbance Response Analysis 99 Table C-6: Complete CRONE set point & output disturbance response data (Group 3 & 4). | Patient # Performance Measurements Set Point Tr Ts OS bninl [min] \%\ Disturbance Ts OS IAE fminl [%1 [-1 28 29 30 31 32 33 34 35 I f f 37 38 39 40 41 42 43 44 5.78 6.14 4.65 3.65 3.45 3.34 2.93 3.06 3.05 3.93 3.70 5.76 5.79 4.17 6.10 5.66 3.60 5.78 6.14 4.22 9.86 6.17 6.17 7.61 8.79 7.67 3.44 7.29 10.19 8.12 8.05 5.26 8.36 6.39 0.00 0.00 I. 45 8.81 6.56 7.12 II . 84 12.35 12.81 0.57 9.11 3.25 2.02 5.88 1.10 2.18 5.36 17.75 18.45 9.23 4.02 4.05 3.58 2.53 3.09 6.96 4.38 7.07 6.74 7.36 7.15 4.98 5.64 6.84 0.00 0.00 0.00 0.01 0.00 0.28 1.92 1.36 4.41 2.20 2.84 0.61 0.42 2.10 0.55 0.64 2.76 57.26 57.89 43.40 30.24 29.20 26.17 22.04 26.19 23.79 13.93 23.40 34.92 36.44 24.23 23.84 28.92 21.42 Appendix C: Set Point and Disturbance Response Analysis 100 Table C-7: Complete CRONE set point & output disturbance response data (Entire population, patients 1-27). Performance Measurements a .© _w "s a. o OH 4> J . s w Patient Set Point Disturbance # Tr T OS T OS IAE [mini [min] [%] [min] [%] [-] 01 14.13 7.27 0.18 4.58 2.89 13.39 02 15.37 9.24 0.12 2.22 1.83 14.31 03 4.20 13.01 11.67 8.23 8.85 33.95 04 13.50 9.46 0.17 5.44 1.75 36.24 05 6.81 6.02 0.49 9.33 3.38 28.13 06 14.65 9.12 0.16 3.30 1.49 19.90 07 4.08 8.39 3.14 5.74 7.02 21.64 08 16.01 8.98 0.09 5.50 6.23 17.89 09 14.61 8.96 0.16 3.60 1.87 22.66 10 14.91 9.14 0.13 5.94 3.66 17.34 11 14.43 8.52 0.16 7.53 3.30 21.22 12 15.49 9.18 0.11 5.77 4.37 16.64 13 15.01 8.93 0.13 5.76 2.50 16.02 14 14.54 8.70 0.16 4.14 2.23 12.05 15 14.92 8.96 0.13 7.14 2.54 19.92 16 3.88 11.49 10.22 5.33 9.58 23.97 17 14.77 9.23 0.16 2.97 1.71 20.40 18 14.76 9.16 0.16 3.63 1.77 24.09 19 14.99 10.28 0.14 7.37 0.34 33.37 20 14.95 9.48 0.13 4.89 1.08 27.02 21 14.96 9.55 0.15 4.29 1.01 24.35 22 15.81 9.59 0.10 5.82 2.05 19.60 23 15.66 9.41 0.10 6.09 2.99 18.78 24 15.98 9.46 0.10 5.48 2.47 16.85 25 15.28 9.49 0.13 3.44 0.88 17.31 26 15.24 9.50 0.12 3.42 1.92 22.71 27 14.80 9.25 0.16 4.92 1.02 28.18 Appendix C: Set Point and Disturbance Response Analysis 101 Table C-8: Complete CRONE set point & output disturbance response data (Entire population, patients 28-44). Performance Measurements Patient # Set Point Tr Ts OS [mini [mini \%\ Disturbance Ts OS IAE rmini r%i r-i "3 o, o PH ii "28" 29 30 31 32 33 34 35 15.63 15.35 15.20 15.36 15.32 15.44 15.92 TT36" 11.38 10.52 9.60 9.81 9.54 9.48 9.75 7T2T 0.23 0.17 0.14 0.13 0.12 0.12 0.10 10.65 7.38 4.31 4.40 3.94 2.82 3.42 1H6~ 0.21 0.28 0.91 0.72 1.02 1.47 1.26 " T O T 44.05 33.65 24.49 23.79 22.03 17.68 21.13 3 6 " 37 38 39 40 41 42 43 44 T6U5 -16.37 15.15 15.20 16.03 15.32 15.88 15.66 16.39 T3f j 9.76 9.75 10.77 11.02 9.75 10.62 10.68 9.96 u r n 0.09 0.14 0.19 0.14 0.14 0.15 0.15 0.08 "X9T" 2.49 3.69 9.08 10.10 4.58 6.60 7.57 3.60 W 0.83 1.07 0.22 0.18 0.79 0.24 0.22 1.07 ~23~W 14.26 23.80 40.50 43.01 25.59 27.69 33.40 21.87 Appendix D: Performance Simulation Plots D.l: PID Controller 30 40 Time Intimites] Patients . Y sp Disturbance so 60 70 30 40 Time |minutes| Figure D - l : PID performance experiment simulation (Group 1). 102 Appendix D: Performance Simulation Plots 103 Patients Y sp Disturbance 30 40 Time |minutes| 50 60 70 _ 4 B 1 H r 2-0 10 20 30 40 50 60 70 Time [minutes] Figure D-2: PID performance experiment simulation (Group 2). V Patients Y sp Disturbance 10 20 30 40 Time |minutes| 50 60 70 1 • SI 0 10 20 30 40 50 60 70 Time |minutes| Figure D-3: PID performance experiment simulation (Group 3). Appendix D: Performance Simulation Plots Time [minutes] Figure D-4: PID performance experiment simulation (Group 4). 0 10 20 30 40 50 60 70 Time |minutes] I | 3 -"3D a 2 0 10 20 30 40 50 60 70 Time |minutes] Figure D-5: PID performance experiment simulation (Entire population). Appendix D: Performance Simulation Plots 105 D.2: CRONE Controller 20 .1 10 20 30 40 Time [minutes| 50 Patients . Y sp Disturbance 50 70 30 40 Time (minutes| Figure D-6: CRONE performance experiment simulation (Group 1). Appendix D: Performance Simulation Plots 80 h i 60 10 20 30 40 Time [minutes] 50 — Patients . . Y sp Disturbance 60 70 30 40 Time [minutes] 50 60 70 Figure D-7: CRONE performance experiment simulation (Group 2). 0 10 20 30 40 50 60 70 Time |minutes] Figure D-8: CRONE performance experiment simulation (Group 3). Appendix D: Performance Simulation Plots 1 Patients Y sp Disturbance 50 60 70 50 60 70 Figure D-9: CRONE performance experiment simulation (Group 4). Time jminutesj Time [minutes] Figure D-10: CRONE performance experiment simulation (Entire population). Appendix E : Acronyms A E P Auditory Evoked Potentials BIS Bispectral Index CNS Central Nervous System C R O N E 'Commande Robuste d'Ordre Non Entier' D O A Depth of Anesthesia E D T A Disodium Edetate E E G Electroencephalogram G A B A Gamma-aminobutyric Acid IAE Integrated Absolute Error IE Integrated Error L M A Laryngeal Mask Airway LTI Linear Time Invariant M D A P E Media Absolute Performance Error M D P E Media Performance Error MINO Multiple-Input Multiple-Output M L A E P Midlatency Auditory Evoqued Potentials M P C Model Predictive Controller PD Pharmacodynamic Model PE Performance Error PI Proportional-Integral Controller PID Proportional-Integral-Derivative Controller P K Pharmacokinetic Model 108 Appendix E: Acronyms 109 SISO Single-Input Single-Output W A V Wavelet-based Anesthetic Value
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Robust control : PID vs. fractional control design,...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Robust control : PID vs. fractional control design, a case study Martínez, Arturo 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Robust control : PID vs. fractional control design, a case study |
Creator |
Martínez, Arturo |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | This thesis presents a systematic procedure to design PID and CRONE controllers to regulate the hypnotic state of anesthesia with the intravenous administration of propofol. The hypnotic state of the patient is assessed by means of the WAV index, a processed hypnotic index based on a deterministic analysis of the EEG. The objective of the controllers is to provide an adequate drug administration regime of propofol, depending on the hypnotic set point, to avoid under or over dosage of the patients. Propofol is assumed to be administered by a commercial Graseby 3400 infusion pump. The two model-based controllers are designed to compensate for the patients inherent drug-response variability (uncertainty), to achieve good output disturbance rejection, and to attain good set point response. The drug-response model consists of two parts: a pharmacokinetic model that characterizes the drug movement in the body, and pharmacodynamic model that relates the drug concentration in the brain to the WAV hypnotic index. An anti-windup scheme is also implemented to protect the system against performance degradation in the event of actuator saturation. The performance of the controllers is assessed by calculating typical time domain measures (overshoot, settling and rise times), and using the median performance error, median absolute performance error, divergence, and wobble of the system's output. Overall, the CRONE controller was shown to have better output disturbance rejection, and to be less sensitive to the drug-response variability of the patients. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0064957 |
URI | http://hdl.handle.net/2429/17568 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2006-0075.pdf [ 6.6MB ]
- Metadata
- JSON: 831-1.0064957.json
- JSON-LD: 831-1.0064957-ld.json
- RDF/XML (Pretty): 831-1.0064957-rdf.xml
- RDF/JSON: 831-1.0064957-rdf.json
- Turtle: 831-1.0064957-turtle.txt
- N-Triples: 831-1.0064957-rdf-ntriples.txt
- Original Record: 831-1.0064957-source.json
- Full Text
- 831-1.0064957-fulltext.txt
- Citation
- 831-1.0064957.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0064957/manifest