Received 19 March 2014; revised 19 April 2014; accepted 25 April 2014

ABSTRACT

This paper presents the automatic drug administration for the regulation of bispectral (BIS) index in the anesthesia process during the clinical surgery by controlling the concentration target of two drugs, namely, propofol and remifentanil. To realize the automatic drug administration, real clini- cal data are collected for 42 patients for the construction of patients’ models consisting of pharma- cokinetic and pharmacodynamic models describing the dynamics reacting to the input drugs. A nominal anesthesia model is obtained by taking the average of 42 patients’ models for the design of control scheme. Three PID controllers are employed, namely linear PID controller, type-1 (T1) fuzzy PID controller and interval type-2 (IT2) fuzzy PID controller, to regulate the BIS index using the nominal patient’s model. The PID gains and membership functions are obtained using genetic algorithm (GA) by minimizing a cost function measuring the control performance. The best train- ed PID controllers are tested under different scenarios and compared in terms of control performance. Simulation results show that the IT2 fuzzy PID controller offers the best control strategy regulating the BIS index while the T1 fuzzy PID controller comes the second.

Keywords:

Anesthesia, Propofol, Remifentanil, BIS, PID, Fuzzy, Control, Simulation, Close-Loop

1. Introduction

Anesthesia is a reversible state of people who are temporarily lack of consciousness in the purpose of undergo- ing a surgery without pain. General process of anesthesia consists of induction, maintenance, emergence, and re- covery [1] . The depth of anesthesia (DOA) represents the level of consciousness [2] , which is to be controlled through the process of anesthesia for varieties of control targets such as steady error, settling time, and overshoot. In general, control strategies can be categorized into two classes: open-loop control and closed-loop control. In open-loop control, based on the knowledge and experience, anesthetists manually adjust the drug dosage to maintain the DOA assisted by some clinical indices of patients. In closed-loop control, the drug dosage is auto- matically adjusted according to some indices of DOA, which makes control input continuous and responsive [3] . Closed-loop control is also expected to avoid over- and under-dosage and suppress the adverse effect of inter- individual differences [4] . Despite of its advantages, the stability of closed-loop control needs to be ensured due to the automated process without supervision [5] .

To achieve the stabilization based on closed-loop control, the following components are required [6] : 1) a patient model, with an output as the index of DOA; 2) a controller for stabilization, such as proportional-integral- differential (PID) controllers. Therefore, to begin with, an estimated patient model is required to represent real patients supporting the control design and performance evaluation. A widely employed mathematical model consists of a linear pharmacokinetic (PK) model and a nonlinear pharmacodynamic (PD) model [7] . The PK mo- del uses the drug dosage as the input and the drug concentration as the output. In what follows, the PD model exploits the drug concentration as the input and exports the index of DOA.

As stated above, the control input of the overall patient model is the drug dosage. According to the types of drugs, the anesthesia can be classified to inhalational and intravenous anesthesia. Inhalational anesthesia has been used since the mid-19th century with for instance nitrous oxide, isoflurane, and halothane as commonly used drugs [8] . In contrast, intravenous anesthesia with propofol and remifentanil as commonly used drugs [9] has the following advantages [10] : separate provision of anesthesia from ventilation, reduced atmospheric pollution, rapid and clear-headed recovery, and so on. While inhalational anesthesia is still frequently used for children, in- travenous anesthesia becomes more and more popular due to the rapid and safe transition [11] .

Main difficulties involved in modeling are the determination of PD model parameters and the selection of the output index. In general, the parameters of PK and PD models need to be determined beforehand. For the PK mo- del, parameters can be estimated depending on sex, age, and weight of patients. Nonetheless, for the PD model, it is not possible to estimate the parameters for certain patient. Accordingly, it requires the controller to be robust in a domain of PD model parameters [12] . Although PK model parameters are varied with patients and PD model parameters are even changing for one patient, a general index can be designed to evaluate the DOA for all patients [13] . The bispectral (BIS) index [14] is an extensively accepted index to measure the DOA. It is a uni- variate dimensionless parameter from 0 to 100, which is devised from a large set of electroencephalogram (EEG) data [15] .

With regard to the controllers for stabilization, several control problems are faced for the closed-loop control of anesthesia: stability, which is the basic control objective; robustness, which overcomes the uncertainty of PD model parameters, measurement noise, surgical stimulation, and so forth; adaptiveness, which makes the controller adaptive to different patients rather than only one patient. Commonly used controllers include PID controllers, model-based controllers, and knowledge-based controllers [6] .

To close the loop of control systems, classical PID controller is a doubtless option due to its successful applications in other areas. For stabilization, it is necessary to tune the parameters of PID controllers for a certain patient model. While the tuning process is an empirical manner, some tuning schemes were even developed to gua- rantee the robustness [16] [17] . Combined with patient model identification from the induction phase of anesthe- sia [18] , the adaptiveness can be further achieved. In addition, combined with genetic algorithm (GA), the para- meters of PID controllers can also be online optimized [19] .

Despite of various causes of the changes in PK and PD model parameters, the model-based controller relies on the current model which reflects the patient’s current pharmacological behavior. For this reason, the patient model needs to be updated through the overall process of anesthesia. This online adaptation can be achieved by several approaches such as Kalman filter algorithm [20] , Bayesian-based adaptive control [21] , and adaptive ge- netic fuzzy clustering algorithm [19] . Alternatively, the variability of PK and PD models between individuals as well as surgical stimulation and anesthetic-analgesic interaction can be explicitly considered offline such that the stability is mathematically guaranteed [22] . Additionally, a Lyapunov-based adaptive controller was developed to attain partial asymptotic regulation [23] .

Unlike model-based controllers, knowledge-based controllers do not require a known mathematical model. In fuzzy logic controllers, for example, decisions are made based on fuzzy rules predefined by expert knowledge

and experience. Compared with linear PID controllers and model-based controller, knowledge-based controllers are easier to implement without tuning PID parameters or mathematical derivation. A hybrid control scheme, namely fuzzy-PID control, was developed to combine the merits of both control strategies [24] . Moreover, a multivariable neural-fuzzy controller was proposed to simultaneously administrate both propofol and remifentanil [25] . However, one problem of knowledge-based controllers is that the interaction of each piece of knowledge makes the controller less transparent and makes it difficult to achieve adaptive control [6] . To deal with this problem, a direct adaptive interval type-2 (IT2) fuzzy logic controller was proposed for multivariable anesthesia systems [26] and a genetic fuzzy logic controller was developed to adjust fuzzy rules using GA [19] .

Due to the advantages of closed-loop control, in this paper, we aim to implement the closed-loop control for anesthesia model. The performance of controllers is ensured by simulation and optimization. For the anesthesia model used in this paper, while existing anesthesia model is utilized, the parameters of PD model are specially obtained from clinical data collected from 42 patients. Following the same design procedure, these tailor-made controllers can be redesigned for other patients as long as clinical data have been collected. In the anesthesia model, the co-administration of both propofol and remifentanil is investigated. To realize the drug administration, we propose two groups of strategies: two fuzzy PID controllers and one fuzzy PID controller with scaling factors. For each group of strategies, linear PID controller, type-1 (T1) fuzzy PID controllers and IT2 fuzzy PID controllers are designed. To draw a distinction from existing IT2 fuzzy PID controllers, we combine IT2 fuzzy PID controllers with GA. All PID gains, scaling factors and parameters of membership functions are optimized by GA in an offline manner subject to a performance index (cost function) which quantifies the performance of the controllers. A BIS training profile which sets different local targets for regulation considering the real anesthesia situation is employed for training purposes. The trained PID control strategies will be tested by a testing profile to verify their performance towards unseen working conditions. Comparisons are made among all controllers to demonstrate the characteristics of each control strategy.

Following sections are organized in this sequence: in Section 2, a multivariable anesthesia model used in this paper is described; in Section 3, the control background of linear PID controller, T1 fuzzy PID controllers and IT2 fuzzy PID controllers is introduced; in Section 4, we design all the control strategies and overall procedure and necessary information of the simulation; in Section 5, simulation results are provided and comparison and analysis are carried out; finally in Section 6, a conclusion is drawn.

2. Multivariable Anesthesia Modeling

In this section, the anesthesia model with two drugs is introduced. As shown in Figure 1, the anesthesia model consists of the target controlled infusion (TCI) system, the PK model, and the PD model, where and are plasmatic concentration targets, and are infusion rates, and are plasmatic concentrations, and are effect-site concentrations, and the subscript and are for propofol and remifentanil, respectively. Details of each module are described in the following subsections. The clinical protocol for obtaining the parameters of PD model is also described.

2.1. Pharmacokinetic Modeling

PK model is applied to reckon the relation from the drug infusion rate and the amount of drugs to plasmatic concentration and effect-site concentration. In this paper, we utilize the PK model from Marsh [27] and Minto [28] for propofol and remifentanil, respectively. A three-compartmental structure of these PK models is presented in Figure 2 which can be formulated by the following dynamic system:

(1)

is the continuous time; is the drug infusion rate; is the rate constant for dis- tribution and elimination of the anesthetic drug; is the volume of central compartment; and is a scalar describing the delay on time of the plasmatic-effect-site drug concentration equilibration.

2.2. Pharmacodynamic Modeling

After and are obtained from PK model, PD model maps to the BIS value which reflects the anesthetic condition of patients. In general, PD model is established by a sigmoid model named as Hill curve for single drug. Based on the Hill curve, a response surface model was proposed by Minto et al. [30] to characterize drugs’ interaction for multi-drugs. The PD model for representing propofol and remifentanil can be formulated as follows [31] :

(2)

where is the clinical effect when no drug is infused; is the maximum drug effect; is the steepness of the concentration-response curve; is the ratio of the two drugs defined in Equation (3); and are the normalized and, respectively, defined in Equation (4); is the effect-site drug concentration associated with 50% maximal drug potency; is defined in Equation (5); is a predefined scalar.

(3)

(4)

(5)

In this paper, let us define and as 98 and 0, respectively., , and are patient-dependent parameters which are estimated by clinical data explained in the following subsection.

2.3. Target Controlled Infusion System

TCI system transfers the plasmatic concentration targets to corresponding drug infusion rates. Practically, the infusion can be realized by Alaris PK syringe pumps which incorporates the TCI system. These pumps are able to infuse adequate amount of drugs according to predefined. Furthermore, they can predict the realistic concentration based on PK model.

The TCI system consists of a bolus-elimination-transfer (BET) infusion scheme designed from the PK model so the is achieved. The initial bolus, loading dose (LD), is calculated according to Equation (6) where is the volume of central compartment represented on Figure 2. The LD bolus will ensure the is achieved however it needs to be followed by a continuous infusion, , defined by Equation (7) to ensure is maintained due to elimination and compartmental transfer of the drug [32] .

(6)

(7)

where is the continuous time; is the rate constant for distribution and elimination of the anesthetic drug defined in Section 2.1.

TCI systems are also able to deal with changes. When increases an initial additional loading dose (ADDLD), Equation (6) and Equation (7) will become the following equations to maintain the new.

(8)

(9)

where is the time of change; are the rate constant for distribution and elimination of the anesthetic drug; is the amount of anesthetic drug in the second compartment at the time of target

change; is the amount of anesthetic drug in the third compartment at the time of target change.

The new continuous infusion rate formula takes into account the amount of drug already present in the peripheral compartments and which will be redistributed.

In the case when the decreases, no additional bolus is necessary, however the infusion of anesthetic drug is halted until is achieved on central compartment due to elimination of drug presented on the PK model. Once this is achieved the continuous infusion restarts according to Equation (9) where is the time is achieved on central compartment through the elimination process.

2.4. Clinical Protocol

The parameters of PD model in this paper are determined by clinical data obtained from 42 elderly patients who had major vascular surgery of the lower limb. The information of these patients is summarized as follows (in the format of “mean ± standard deviation”): 35 males and 7 females, age 72 ± 8 years old, height 169 ± 10 cm and weight 73 ± 14 kg.

The process of collecting data from these patients is briefly described here. Before induction, all sensors were settled and a radial arterial line was inserted to receive baseline readings on all clinical monitors. Two Alaris PK syringe pumps were employed to implement the total intravenous anesthesia (TIVA) with propofol [27] and re-

mifentanil [28] TCI. To begin with, was set to be 3 ng/ml. Once reached 2 ng/ml, was altered to 2 ng/ml and was started with 3 mg/ml. From then on, and were adjusted by anesthetists according to readings of different clinical apparatuses to maintain the BIS value within [40,60] for an appropriate DOA. Among these readings, cardiac output (CO) was monitored and stabilized at or near pre-induction level by anesthetists through all periods of the surgery.

3. PID Controllers

In this section, 3 types of PID controllers, namely linear PID controller, T1 fuzzy PID controller and IT2 fuzzy PID controller, are introduced. These 3 types of PID controllers are employed to regulate the output BIS index of the anesthesia model.

3.1. Linear PID controller

A linear PID controller is shown in Figure 3, which consists of 3 elements, namely proportional, integral and derivative blocks. The output of the linear PID controller in discrete time is given as

(10)

where is the sampling time; is the interval of sampling time; is the output; and are constant proportional, integral and derivative gains, respectively, to be determined.

3.2. Type-1 Fuzzy PID Controller

In view of the linear PID controller [33] , as the proportional, integral and derivative gains are constant, it is not able to handle well a highly nonlinear system. It motivates the use of T1 fuzzy PID controller [34] [35] of which the gains are changing according to the operating domains. By applying different sets of gains in different operating domains, more appropriate PID controller is employed to deal with the nonlinear system resulting in an improvement of control performance.

A T1 fuzzy PID control system is shown in Figure 4, which consists of a T1 fuzzy PID controller and a patient’s model (detailed in Figure 1) connected in a closed loop. Unlike the linear PID controller having a set of constant gains, the T1 fuzzy PID controller has a fuzzy inference system providing a set of feedback gains through a reasoning process according to the operating condition.

The behavior of the fuzzy inference system is governed by a set of fuzzy rules of the following format:

(11)

where is a linguistic variable; is the fuzzy term corresponding to the linguistic variable in the rule; is a positive integer; is the output of the fuzzy inference

system; and is the singleton membership function corresponding to the rule. The inferred output is gi- ven as

(12)

where

(13)

(14)

(15)

where denotes the number of rules; is the normalized grade of membership;

is the grade of membership corresponding to the fuzzy term.

The output of the fuzzy inference system Equation (12) is employed to replace the gains of the linear PID controller turning it to become a T1 fuzzy PID controller. More precisely, 3 fuzzy inference systems are required to implement a T1 fuzzy PID controller. The outputs of the 3 fuzzy inference systems will be employed as, and. Consequently, the PID gains are no longer constant but dependent on the operating condition characterized by the membership functions.

3.3. Interval Type-2 Fuzzy PID Controller

Type-2 (T2) fuzzy sets [36] -[39] demonstrate a superior characteristic handing uncertainties compared with the T1 fuzzy sets. Uncertainties are captured by the lower and upper membership functions which form the footprint of uncertainty (FOU). A T2 fuzzy inference system can be considered as a set of infinite number of T1 fuzzy inference systems. Consequently, T2 fuzzy inference system is able to outperform the T1 fuzzy inference system in terms of reasoning and generalization capability. In general, the defuzzification process for general T2 fuzzy sets is computational demanding. By using IT2 fuzzy sets [40] , the computational demand can be significantly reduced.

By employing interval fuzzy sets for the T1 fuzzy inference system Equation (12), it becomes an IT2 fuzzy inference system. The behavior of an IT2 fuzzy inference system is described by a set of rules of the following format:

(16)

(17)

(18)

where

(19)

(20)

(21)

(22)

(23)

in which denote the lower grade of membership, upper grade

of membership, lower membership function and upper membership function, respectively. Using the center-of- set type reducer, the inferred output of the IT2 fuzzy inference system is given as

(24)

where can be obtained using the Karnik-Mendel (KM) algorithms [41] . The final defuzzifized output is given by

(25)

The IT2 fuzzy PID control system [42] can be represented by Figure 4 as well of which the PID gains are given by the outputs of 3 IT2 fuzzy inference systems.

4. Simulation Design

In this section, the simulation environment is presented. First, we design the training profiles and testing profiles which all controllers need to deal with. The purpose of selected target profiles is first described. Then different control strategies are provided for comparison. For fuzzy control strategies, the approach to determine fuzzy rules is offered. Finally, the procedure of using GA for optimization is described, and a performance index is designed as the cost function for optimization.

4.1. Target Profiles

A training profile as shown in Figure 5, which is a series of BIS values required for an appropriate DOA, is employed for the training of the PID control strategies using GA. It defines the target BIS values in different periods that a PID control strategy has to achieve. After the training, a testing profile as shown in Figure 6 is employed to verify if the trained PID control strategies are able to control the BIS value to reach some targets subject to different operating conditions.

Clinically, the induction process starts with. As a result, the initial condition for the simulation is and the target profile begins by, where is the target BIS at sampling time. At seconds, the induction is executed and. Since the induction period is short compared with the maintenance period, there is not much difference during induction in terms of control performance between various control strategies. Thus, we use two linear PID controllers for two drugs to drive the BIS from 98 to 50. By trial and error, the PID gains are predefined, and it is guaranteed that is

around 50 and the maintenance process starts from seconds.

At seconds, designed controllers replace the predefined controller to achieve the maintenance process. For a suitable anesthesia, it is required that and an acceptable range is. Hence, the control objective is to maintain the BIS value within [40,60]. In order to make the control difficult enough to differentiate various control strategies, we add more targets into the profiles such as

and in Figure 5.

As for the recovery process, since it is not allowed take the drug out of humans, the cease of controllers is the only and fastest option. These will not exist overshoot either because the target BIS is the maximum value. Therefore, all controllers will perform the same at this stage. It is not necessary to design controllers for the recovery process which is thus not included in the target profiles. It is noted that although the induction process is not under comparison either, we cannot ignore it. The reason is that it is difficult to find all model parameters describing the state of. On the other hand, the parameters for are known and it is easy to start from this initial condition.

4.2. Control Strategies

In this paper, we aim to compare different control strategies to achieve a better control performance. All six control strategies are listed in Table 1. These six cases can be separated into two groups: two-controller Cases 1 to 3 and one-controller Cases 4 to 6. In two-controller cases, we control each drug by an independent PID con-

troller offering control outputs and. In one-controller cases, however, the control output of the

PID controller is divided into two control outputs by two scaling factors and. Apparently, parameters of these two divided controllers are constrained by the ratio of these two factors. Theoretically, this constraint leads to conservativeness and make the performance worse than the ones of two-controller cases. Nevertheless, one-controller cases have less number of parameters to be determined resulting in lower computational demand, faster convergence in the training process and lower implementation cost for the PID control strategies.

In each of these two groups, we have the following control strategies: linear PID controllers, T1 fuzzy PID controllers, and IT2 fuzzy PID controllers. For fuzzy PID controllers, the block diagram of BIS index regulation using two controllers and one controller with scaling factors are shown in Figure 7 and Figure 8, respectively. Referring to the figures, denotes the desired BIS target value, is the output BIS value from theanesthesia model, denotes the difference between the desired target and actual BIS value at time

. The PID gains and scaling factors

are decided by the fuzzy inference systems. The PID controller will generate the control signals

and for BIS regulation where and are in Figure 1.

4.3. Fuzzy Rules

T1 and IT2 fuzzy PID controllers are discussed in Subsections 3.2 and 3.3, respectively. While the number of fuzzy rules and the shape of membership function are predefined, the parameters of input and output membership functions are optimized by GA.

In this paper, the shape of membership functions is defined as triangular shape for simplicity. The triangular IT2 membership functions are shown in Figure 9. For each input IT2 membership functions, the lower and upper membership functions are characterized by seven points which are to be optimized by GA. As for T1 fuzzy PID controller, each input membership function is characterized by three points which are to be optimized by GA. For both T1 and IT2 fuzzy PID controllers, their output membership functions are T1 and IT2 singleton membership functions whose values are to be determined by GA.

In fact, from the determination of membership functions and fuzzy rules, we can find the relation between these control strategies: solutions from linear PID controllers can be implemented by T1 fuzzy PID controllers, and solutions from T1 fuzzy PID controllers can be implemented by IT2 fuzzy PID controllers. In other words, linear PID controller is a subset of T1 fuzzy PID controller, and T1 fuzzy PID controller is a subset of IT2 fuzzy PID controller. By adding some constraints on the membership functions and consequents in fuzzy rules, IT2 fuzzy PID controllers can be reduced to T1 fuzzy PID controllers, and T1 fuzzy PID controllers can be reduced to linear PID controllers.

Three rules are employed for each fuzzy inference system. The premise of each rule takes only one linguistic

variable, i.e., where is the target BIS value at time and

is the actual BIS value. Three linguistic terms, namely negative (N), zero (Z), and positive (P), are employed to characterize the premise variable. More rules can be used to partition the universe of discourse for a better result. However, it will lead to slower convergence of training and more computational burden. Additionally, the rate of change of can also be treated as another linguistic variable together with. Likewise, it will cause more fuzzy rules and parameters which increase the computational burden and defers the convergence.

From the above discussion, triangular membership functions and three rules are employed for both T1 and IT2 fuzzy PID controllers. In the GA optimization, the triangular shape and sequence of membership functions of three rules should be guaranteed. For the sequence, specifically, the same point in membership functions corresponding to linguistic terms N, Z and P should be in ascending order (except). For example, , where, , and are points for membership functions corresponding to linguistic terms N, Z and P, respectively.

Another condition needed to be ensured is that at least one rule is fired for. Between two adjacent membership functions, as long as they intersect each other, there exists one rule to be fired. On the left-hand side of membership function N and the right-hand side of membership function P, there are two ap-

proaches to ensure this condition. One is defining and, and another is defining

and. It can be found that solutions from the first approach can be imple-

mented by the second approach. In other words, the second approach is more general than the first approach. Hence, we adopt the second method. For that reason, rules N and P are two shoulder-shape membership functions which can be treated as a special case of triangular membership functions. In this approach, it is guaranteed that at least one rule is fired for, although only is concerned practically.

Note that under this method the equalities and are temporarily employed

in the optimization.

4.4. Parameters Optimization

The PID gains, scaling factors and membership functions are optimized by GA subject to a cost function reflecting the control performance. Due to the disparity of each training, GA is run for 10 times for each case of PID control strategies as shown in Table 1. The best set of parameters for each control strategy for each run of GA is recorded for further comparison, analysis and practical application. Statistical information including the worst, mean and best costs and the standard deviation for the 10 runs are collected. Among the 10 runs for each PID control strategy, the best set of parameters given by the best cost is used to implement the corresponding PID controller. During the optimization, the lower and upper bounds (LB and UB) of PID gains and scaling factors are listed in Table 2, which are determined by trial and error for good control performance.

Due to a large number of variables and the highly nonlinear cost function (defined in the following subsection 4.5), especially for fuzzy PID controllers, GA may not be able to reach the global optimal solution. To speed up the training process, we define the initial population for fuzzy PID controllers based on our knowledge on the PID control strategy. It is known that linear PID controller is a subset of T1 fuzzy PID controller, and T1 fuzzy PID controller is a subset of IT2 fuzzy PID controller. It is thus reasonable that the best set of PID gains for linear PID controllers obtained by GA is employed as the initial PID gains for T1 fuzzy PID controllers for all rules. As a result, initially, the T1 fuzzy controller is equivalent to a linear PID controller. Similarly, the best set of PID gains and membership function parameters obtained for T1 fuzzy PID controller is employed as the initial set of PID gains and membership function parameters for IT2 fuzzy PID controller. That is to say, we utilize the “knowledge” on the PID control strategy such that GA starts with a considerable well initial condition. In this way, although there is still the same number of variables to be trained, GA only needs to utilize the additional parameters and new structures (fuzzy rules and membership functions) provided by fuzzy PID controllers in order to obtain a better cost value.

4.5. Performance Index

The performance index is used to judge whether the control objective is achieved and show the merits of various control strategies. Different performance indices can be selected such as settling time, overshoot, steady error, mean absolute error (MAE), and mean square error (MSE). For the GA optimization, a well-defined performance index is required as the cost function. All parameters of controllers are optimized according to the cost value provided by the cost function. In this paper, we present the following index mainly based on MAE:

(26)

where is the sampling time, and is a positive integer; is the interval of sampling time;, is the target BIS value at time, and;

and are plasmatic concentration targets of propofol and remifentanil, respectively, and

,; are predefined weights.

In the cost function Equation (26), the first term is MAE which aims at minimizing the difference between the current BIS value and its target. Unlike settling time, overshoot and steady error, MAE and MSE record the er-

ror information through all simulation periods which reflect more comprehensive properties of performance. The reason we choose MAE instead of MSE is that MSE amplifies the effect of large errors and then the optimization tries to minimize large errors which results in oscillation. On the other hand, MAE keeps the original weights on large and small errors leading to a mild and smooth response.

The three terms after the first term of (26) are the rate of change of, , and, respectively, which are designed to reduce the oscillation for a gentle progress. The last term is to reduce the difference of the concentration of two drugs because they are equally important for different purposes during the anesthesia. Serious bias to either of these two drugs will not provide an adequate anesthesia.

5. Simulation Results

In this section, we implement the simulation on the control of anesthesia using the anesthesia model in Figure 1. Different control strategies in Table 1 are applied to regulate the BIS using the training profiles shown in Figure 5 for training and their control performance are verified by testing profiles shown in Figure 6. PID gains, parameters of scaling factors and membership functions are optimized by GA according to the cost function Equation (26). Comparisons of performance are made between the six cases. It should be noted that the simulation was carried out in discrete time. The PK model in Equation (1) was discretized by zero-order-hold (ZOH) assumption. The TCI system implemented has also been adjusted to perform under discrete PK model.

A real-coded GA available from Matlab global optimization toolbox is employed for the training. The control parameters of GA are listed in Table 3.

We define the parameters in Equation (26) as follows: seconds,. It is worth to mention that although the weight is not the largest, the MAE is still the main contribution for the total cost. The bounds of parameters in Table 2 are adopted for the training process of GA. By running GA 10 times for each case in Table 1, we obtain the statistical information of the cost as shown in Table 4.

Referring to Table 4, comparing with Cases 1 to 3, Case 3 offers the best cost and Case 2 comes to the second which coincides with the theory. The same rank can be found in Cases 4 to 6. Comparing Cases 1 and 4, two- controller case is better than one-controller case, which can be also proved by comparing between Cases 2 and 5, and Cases 3 and 6. The reason is that one-controller case has the constraint on the ratio of control signal between two drugs while two-controller case does not have such constraint. One-controller case is a subset of two-con- troller case. Although one-controller case performs worse, it has less number of parameters to be determined and thus reduces the computational demand and implementation cost of PID control strategy.

In spite of the “Best” cost, the “Std” is descending from linear PID controllers to IT2 fuzzy PID controller ex- cept Case 4. It provides the information that the convergence gets easier for fuzzy PID controllers. The reason is that we use the best set of parameters of previous cases as the initial populations of following training. It also in- dicates that the improvement becomes more and more difficult, and thus 10 times running offer similar results. The “Std” for two-controller cases is larger than the corresponding one-controller cases, which implies that the convergence for two-controller cases is harder due to more variables to be trained.

The best sets of parameters and membership functions are shown in Tables 5-14. Corresponding membership functions are exhibited in Figures 10-13. From the obtained membership functions, it can be summarized that points and are preferred to be around 20 related to the training profiles we define in Figure 5. Considering that the only linguistic variable is and the maximum change of is in the training profile, therefore, it is reasonable to choose as an estimated domain where the PID controllers work suggesting that the fuzzy blending should start from around. For the IT2 member- ship functions of IT2 fuzzy PID controllers in Figure 11 and Figure 13 after training, it was found that the lower and upper membership functions are very close to each other such as membership functions of N and Z in Figure 11. The values of points associated with the IT2 membership functions are shown in Table 9 and Table 14.

By applying the best sets of parameters to the training profile of BIS, the time responses of BIS and corresponding drug concentration information are shown in Figures 14-19 for Cases 1 to 6, respectively. By applying the best sets of parameters to the testing profile of BIS, the ones for the testing profile are exhibited in Figures 20-25. The green line is drawn by. From the practical clinical point of view, the acceptable

range of BIS is. Therefore, if BIS is stabilized within the region bounded by the

lower and upper green lines, the performance will be well acceptable.

It can be seen from the figures that all cases of PID control strategies succeed in tracking the BIS targets and maintaining BIS values. In terms of control performance, while the linear PID controllers have less steady error, fuzzy PID controllers have less overshoot and settling time. Since there is no rigorous requirement on small steady error, fuzzy PID controllers are preferable that they are able to keep the BIS values within the regions bounded by the green lines with less overshoot and settling time. Both training and testing profiles demonstrate

the same characteristics as discussed for the corresponding control strategies. Recalling a common controller with two linear PID controllers and predefined parameters is used to achieve the induction process for all cases, the proposed PID control strategies kick in at the start of the maintenance process. The switching point of induc- tion and maintenance process occurs at 1000 seconds which results in a small vibration.

The lower panels of Figures 20-25 show the behavior of drug concentration. It can be seen that both and of propofol and remifentanil are adjusted by controllers according to target profiles. More importantly, they are adjusted in advance to decrease the overshoot. This is mainly because all controllers are based on PID control strategy, whose abilities are inherited. In addition, the total amount of remifentanil is larger than propofol,

which illustrates that the remifentanil is more favorable to the stabilization of BIS than the propofol. This might be due to the faster action time of remifentanil comparatively to propofol. Although the total amount of remifentanil is larger, the cost function takes the difference of the amount of two drugs into consideration such that serious bias to remifentanil is avoided. Even though the slightly bias towards remifentanil can be noticed overall the propofol and remifentanil reached are within acceptable clinical boundaries even when the BIS target was set at 20, a value inferior to the safe interval of 40 to 60.

Overall from a clinical point of view the results and performance obtained are acceptable, considering the lack of significant overshoot, the stability of the controller response and the ability of achieving the desired targets on a short time with clinically acceptable. The controller also provided the expected response when the value of BIS target increases, and it stops infusion of anesthetic drugs but only restarts its infusion in other to avoid an overshoot.

When the PID control strategies with the best set of parameters (offering the best cost of) are employed to regulate the BIS value using the testing profile in Figure 6, their cost values are listed in Table 15. The rank of which is identical to the sequence of cost values from the training profile. It indicates that the advantage of fuzzy PID controller over linear PID controller exists not only in the training profile, but also in various target profiles.

6. Conclusion

In this paper, drug administration for anesthesia has been realized by various PID control strategies. We have first constructed a multivariable anesthesia model with propofol and remifentanil based on 42 patients’ clinical data. Simulation has been conducted to regulate the output BIS value governed by this model using six control strategies including linear PID controllers, T1 fuzzy PID controllers and IT2 fuzzy PID controllers. These six strategies are separated into two groups, namely two controllers and one controller with scaling factors, to handle the co-administration of two drugs. Parameters of these controllers have all been optimized by GA subject to a performance index quantitatively measuring the control performance. In order to make the control sophistica- ted, target profiles have been designed for BIS regulation. Different target profiles have been utilized to test and verify the performance of controllers. It has been demonstrated that the IT2 fuzzy PID controllers offer the best performance, and T1 fuzzy PID controllers come to the second. In the future, further investigation can be carried out to control a general and non-parametric multivariable anesthesia model rather than the model with parameters obtained from specific group of patients.

Acknowledgements

The work described in this paper was partly supported by King’s College London, Fundação para a Ciênciae Tecnologia (grant number SFHR/BD/44162/2008) and European Social Fund in the POPH framework.

References