Fractional-Order Control Strategy for Anesthesia–Hemodynamic Stabilization in Patients Undergoing Surgical Procedures

Abstract

Fractional calculus has been opening new doors in terms of better modeling and control of several phenomena and processes. Biomedical engineering has seen a lot of combined attention from clinicians, control engineers and researchers in their attempt to offer individualized treatment. A large number of medical procedures require anesthesia, which in turn requires a closely monitored and controlled level of hypnosis, analgesia and neuromuscular blockade, as well maintenance of hemodynamic variables in a safe range. Computer-controlled anesthesia has been given a tremendous amount of attention lately. Hemodynamic stabilization via computer-based control is also a hot topic. However, very few studies on automatic control of combined anesthesia–hemodynamic systems exist despite the fact that hemodynamics is strongly influenced by hypnotic drugs, while the depth of hypnosis is affected by drugs used in hemodynamic control. The very first multivariable fractional-order controller is developed in this paper for the combined anesthesia–hemodynamic system. Simulation studies on 24 patients show the effectiveness of the proposed approach.

1. Introduction

Fractional calculus has been opening new doors in terms of better modeling and control of several phenomena and processes [1]. Among these, biomedical engineering is considered a hot topic, as researchers and clinicians alike are in search for an improved understanding of how the human body works and an optimal automatic control of drug dosing [2,3]. Driven by technological advances, robot surgeons [4], targeted drug administration using nano-robots [5], computer-assisted drug dosing [6], etc., have emerged. The future of medicine is shaped to follow the Industry 4.0 paradigm and to adopt computer-based monitoring and control [7]. A large number of medical procedures require general anesthesia. This in turn implies an optimal balancing of its three components: hypnosis, analgesia and neuromuscular blockade [8]. At the same time, all other vital parameters, such as hemodynamic variables, have to be adequately monitored and maintained in a safe range [9,10].

Fractional calculus has received a lot of attention with respect to modeling [10,11] and control of the three components in anesthesia [12]. General anesthesia is a broad term encompassing the use of drugs to induce and maintain three states during surgery: hypnosis (depth of unconsciousness), analgesia (absence of pain) and areflexia (neuromuscular blockade). An extensive review of closed-loop control of anesthesia is included in [8].

Depth of hypnosis has been extensively studied and several control algorithms have been proposed to maintain it at a predefined level [13,14,15,16,17,18]. The bispectral (BIS) signal is used to estimate the depth of hypnosis and it is obtained by processing electroencephalogram (EEG) signals in patients. In a control loop, BIS is considered as the controlled output, while propofol, the corresponding drug used to induce hypnosis, is considered as the manipulated input. In addition to propofol, remifentanil and atracurium are used as manipulated variables to control analgesia and neuromuscular blockade (NMB), respectively. The Ramsay Agitation Score (RASS) and the electromyogram (EMG) are used as controlled outputs to monitor analgesia and NMB. In order to properly design a suitable control algorithm, a model of the human body and knowledge regarding how it reacts to anesthetic drugs are required. Pharmacokinetic and pharmacodynamic models are used to describe the patient’s response to a certain drug [8]. In this paper, an existing benchmark patient simulator is used [19]. The model of anesthesia is based on a compartmental approach [20]. The pharmacokinetic model consists in a combination of three compartments: a fast-acting compartment for blood and two slow compartments representing muscle and fat. This is generally represented by a set of differential equations [19]. The pharmacodynamic model reflects the dose–effect response and it is represented by a nonlinear Hill equation. However, research has shown that a synergic dose–effect response exists in the case of propofol and remifentanil upon BIS. In this case, the combined dose–effect relation is represented by complex Hill functions that result in a 3D nonlinear surface. Different control strategies have been proposed for regulating BIS levels in sedated patients, such as PIDs [15], internal model control [21], predictive control [16], event based PID control [14,17], fuzzy control [18] and others [8,22,23].

The neuromuscular blockade component in closed-loop anesthesia has also seen its share of research and simulation and clinical results. Several control methods have been proposed, ranging from various PID-type controllers [24] to advanced control algorithms [25,26]. Combined control of depth of hypnosis and NMB has been clinically evaluated [27]. However, evidence shows that a multivariable approach is not necessary since on one hand, atracurium has no effect upon BIS nor RASS, while on the other, propofol has no influence upon NMB. It is remifentanil that has a minor effect on the NMB, but this is negligeable compared to atracurium [8].

Most anesthesia-related control algorithms discuss simplified approaches, where interactions among the various drugs are not completely tackled. Research has shown that hemodynamic variables, such as mean arterial pressure (MAP) and cardiac output (CO), significantly depend on the anesthesia levels. An increase in the propofol dose decreases both MAP and CO, while an increase in the remifentanil dose decreases MAP and increases CO. These interactions should be tackled by a properly designed multivariable controller. Sodium nitroprusside (SNP) and dopamine (DP) are usually employed to maintain MAP and CO levels at safe values [28]. Other drugs may be used instead, with similar mathematical models. A hemodynamic system, having MAP and CO as the controlled outputs and SNP and DP as corresponding inputs, is modelled using a 2 × 2 multivariable system with first-order-plus-dead-time transfer functions [9], [19]. Control algorithms for the hemodynamic system have also been researched since the early 90s, ranging from simple PIDs to adaptive, optimal controllers, etc. [29,30,31,32]. Studies on the combined hemodynamic and anesthetic systems exist [33,34,35].

Fractional-order PID controllers have been previously considered for controlling biological systems [36,37,38], as well as possible solutions for automated drug dosing in anesthesia to regulate the depth of hypnosis [39,40,41]. Most of the existing control strategies deal with a single-input–single-output perspective of controlling BIS levels in patients using propofol or other drugs. Multiple-input–single-output approaches have also been researched in terms of co-administration of propofol and remifentanil to control BIS [17], [42,43], but not for fractional-order controllers. For regulating NMB, very few fractional-order controllers exist, and these solely consider the atracurium–NMB relation. In terms of maintaining the hemodynamic variables in a safe range, fractional-order controllers have been proposed [44,45], including multivariable ones [9,46].

No research on fractional-order control of combined anesthesia–hemodynamic systems has been reported so far, except for our preliminary study regarding some simple fractional-order PI controllers [47]. In this paper, based on input–output analysis of an existing patient model [19], a computer-based multivariable fractional-order control algorithm is designed. According to interactions between the various loops [19], a decentralized control approach for hemodynamic stabilization is firstly considered [47] and briefly discussed. Next, the problem of maintaining the RASS score constant by manipulating the remifentanil dose is solved using a simple fractional-order PID controller. Once these loops have been closed, a double-switching fractional-order PID controller is developed for controlling the BIS level using propofol. The NMB can be easily controlled in a single-input–single-output framework, being quasi-decoupled from the rest of the anesthesia–hemodynamic system.

The main contribution of the manuscript consists in the very first multivariable fractional-order control strategy for combined anesthesia–hemodynamic systems. The approach for using a fractional-order controller for hypnosis induction and a different fractional-order controller for maintaining it are also original concepts. Additionally, a novel aspect and contribution of the manuscript consists in the design of a fractional-order function for bumpless switching between the induction and maintenance phases of hypnosis.

In clinical practice, an induction phase is considered first, where the anesthesiologist administers a propofol bolus that drives the patient from full awareness (BIS = 100%) to the required hypnotic state. In this case, the BIS index must be in a 40–60% range. The induction phase can also be achieved using an automated control system with the following performance specifications: a settling time of less than 4–5 min and an undershoot as small as possible, which might otherwise cause hypotension. Once the patient has reached a desired BIS level, usually 50%, the surgical procedures can be initiated. This corresponds to the maintenance phase. The anesthesiologist has to administer various boluses to account for the surgical stimulus by reading the patient’s vital signs. A properly designed control system could partially replace the anesthesiologist, whose remaining task would be that of supervising and readjusting the reference value for the BIS signal. The control system would have to ensure that a set of performance specifications is also met during the maintenance phase, such as: fast disturbance rejection due to nociception stimulation and a BIS value between 40–60%. In both the induction and maintenance phases, the control system has to exhibit robustness to intra- and inter-patient variability. Emergence is the third and final phase of hypnosis and consists in stopping the administration of the drugs. In this manuscript, only the first two phases are considered.

In this paper, a benchmark patient model [19] is considered and 24 patients with different characteristics are investigated, as indicated in

Table 1. Representative Patient Database for Propofol-to-BIS with pharmacokinetic Model Biometric Values and PD Model Sensitivity Values [19].

Index	Age (yrs)	Height (cm)	Weight (kg)	C50 (mg/ml)	γ (-)
1	74	164	88	2.5	3
2	67	161	69	4.6	2
3	75	176	101	5	1.6
4	69	173	97	1.8	2.5
5	45	171	64	6.8	1.78
6	57	182	80	2.7	2.8
7	74	155	55	1.7	3.5
8	71	172	78	7.8	2.9
9	65	176	77	2.9	1.88
10	72	192	73	3.9	3.1
11	69	168	84	2.3	3.1
12	60	190	92	4.8	2.1
13	61	177	81	2.5	3
14	54	173	86	2.5	3
15	71	172	83	4.3	1.9
16	53	186	114	2.7	1.6
17	72	162	87	4.5	2.9
18	61	182	93	2.7	1.78
19	70	167	77	6.8	3.1
20	69	168	82	9.8	1.6
21	69	158	81	3.2	2.1
22	60	165	85	5.1	2.51
23	70	173	69	3.67	3.1
24	56	186	99	5.8	2.3

. An open-loop analysis of the BIS levels in the 24 patients (Figure 1) due to a step change in the propofol dose suggests the need for robust fractional-order controllers to tackle patient intra- and inter-variability.

2. A Dedicated Double-FOPID Control System with Bumpless Transfer for Controlling the Depth of Hypnosis

In the literature, regarding general anesthesia control, various studies suggest separating the induction phase from the maintenance phase in terms of control logic. Controlling these independently is much more optimal and specific. During the induction phase, a combination of anesthetic and opioid agents involves a strong synergistic effect that has to be unified for many different patients, whereas in the maintenance phase, the controller should reject the surgical stimuli as fast as possible. Thus, we define the induction (a,b,c) and maintenance (d) phase control specifications as follows:

(a)

The intra- and inter-patient variability must be minimized;

(b)

The time to target must respect the four-minute mark;

(c)

Excessive undershoot and significant oscillations must be avoided;

(d)

Disturbances acting on the depth of anesthesia such as nociceptor stimuli must be rejected, keeping BIS in the range of 40% to 60%.

In order to fulfil these conditions, the general fractional-order PID (FOPID) controller is considered:

$${\mathrm{H}}_{\mathrm{c}}{\left(\mathrm{s}\right)=\mathrm{k}}_{\mathrm{p}}\left(1+\frac{{\mathrm{k}}_{\mathrm{i}}}{{\mathrm{s}}^{\mathsf{\lambda }}}{+\mathrm{k}}_{\mathrm{d}}{\mathrm{s}}^{\mathsf{\mu }}\right)$$

(1)

where ${\mathrm{k}}_{\mathrm{p}}$, ${\mathrm{k}}_{\mathrm{i}}$, ${\mathrm{k}}_{\mathrm{d}}$ are the proportional, integral and derivative gains, respectively, and $\mathsf{\lambda }$, $\mathsf{\mu }$ are the integral and derivative orders. The estimation of the FOPID parameters needs to take into account the design specifications mentioned in a–d. Since tuning FOPID controllers can be simplified by considering the frequency domain representation of (1), the time-domain performance specifications are converted to their frequency-domain counterparts.

For the induction phase, the time to target is addressed via a gain crossover frequency. The larger this is, the faster the settling time is and, implicitly, the time to target. Excessive undershoot is then addressed via a large phase margin, which ensures that the overall control systems cannot be easily destabilized. A constraint regarding the robustness of the control system to ensure as little as possible undershoot in the case of inter- and intra-patient variability is addressed via the iso-damping property of the designed controller. In this way, the phase margin is kept constant despite gain variations. Recent studies [19,47] regarding the BIS variability for 24 patients suggests that gain variations occur frequently (Figure 1) and should be addressed in a controlled anesthesia framework.

Having

$${\mathrm{s}}^{\mathsf{\mu }}={\left(\mathrm{j}\mathsf{\omega }\right)}^{\mathsf{\mu }}{=\mathsf{\omega }}^{\mathsf{\mu }}\left(\cos \frac{\mathsf{\mu }\mathsf{\pi }}{2}+\mathrm{jsin}\frac{\mathsf{\mu }\mathsf{\pi }}{2}\right){, \mathrm{s}}^{-\mathsf{\lambda }}={\left(\mathrm{j}\mathsf{\omega }\right)}^{-\mathsf{\lambda }}\left(\cos \frac{\mathsf{\lambda }\mathsf{\pi }}{2}-\mathrm{jsin}\frac{\mathsf{\lambda }\mathsf{\pi }}{2}\right)$$

(2)

gives the real, Re, and imaginary, Im, parts of the complex FOPID by substituting (2) in (1):

$${\mathrm{Im}\left(\mathsf{\omega }\right)=\mathrm{k}}_{\mathrm{p}}\left({\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}^{\mathsf{\mu }}\sin \frac{\mathsf{\mu }\mathsf{\pi }}{2}{-\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}^{-\mathsf{\lambda }}\sin \frac{\mathsf{\lambda }\mathsf{\pi }}{2}\right){, \mathrm{Re}\left(\mathsf{\omega }\right)=\mathrm{k}}_{\mathrm{p}}\left({1+\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}^{-\mathsf{\lambda }}\cos \frac{\mathsf{\lambda }\mathsf{\pi }}{2}{+\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}^{\mathsf{\mu }}\cos \frac{\mathsf{\mu }\mathsf{\pi }}{2}\right).$$

(3)

The controller magnitude, phase and derivative of phase at the gain crossover frequency are determined according to:

$$\left|{\mathrm{H}}_{\mathrm{c}}(\mathrm{j}{\mathsf{\omega }}_{\mathrm{c}})\right|=\sqrt{\mathrm{Re}({\mathsf{\omega }}_{\mathrm{c}}{)}^2+\mathrm{Im}({\mathsf{\omega }}_{\mathrm{c}}{)}^2}, \angle {\mathrm{H}}_{\mathrm{c}}\left(\mathrm{j}{\mathsf{\omega }}_{\mathrm{c}}\right)=\mathrm{atan}\left(\frac{\mathrm{Im}({\mathsf{\omega }}_{\mathrm{c}})}{\mathrm{Re}({\mathsf{\omega }}_{\mathrm{c}})}\right), \mathsf{\varphi }={\frac{\mathrm{d}\mathsf{\angle }{\mathrm{H}}_{\mathrm{c}}\left(\mathrm{j}\mathsf{\omega }\right)}{\mathrm{d}\mathsf{\omega }}|}_{\mathsf{\omega }={\mathsf{\omega }}_{\mathrm{c}}}$$

(4)

The frequency-domain performance specifications yield the following system of non-linear equations:

$$\{\begin{array}{c}{\mathrm{k}}_{\mathrm{p}}\sqrt{{\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\mu }}\left(\sin \frac{\mathsf{\mu }\mathsf{\pi }}{2}+\cos \frac{\mathsf{\mu }\mathsf{\pi }}{2}\right)+{\left({\mathsf{\omega }}_{\mathrm{c}}{}^{-\mathsf{\lambda }}{\mathrm{k}}_{\mathrm{i}}\cos \frac{\mathsf{\lambda }\mathsf{\pi }}{2}+1\right)}^2-{\left({\mathsf{\omega }}_{\mathrm{c}}{}^{-\mathsf{\lambda }}{\mathrm{k}}_{\mathrm{i}}\sin \frac{\mathsf{\lambda }\mathsf{\pi }}{2}\right)}^2}-\frac{1}{\left|\mathrm{P}\left({\mathrm{j}\mathsf{\omega }}_{\mathrm{c}}\right)\right|}=0\\ \frac{{\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\mu }}\sin \frac{\mathsf{\mu }\mathsf{\pi }}{2}{-\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}_{\mathrm{c}}{}^{-\mathsf{\lambda }}\sin \frac{\mathsf{\lambda }\mathsf{\pi }}{2}}{{1+\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}_{\mathrm{c}}{}^{-\mathsf{\lambda }}\cos \frac{\mathsf{\lambda }\mathsf{\pi }}{2}{+\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}_{\mathrm{c}}{}^{-\mathsf{\mu }}\cos \frac{\mathsf{\mu }\mathsf{\pi }}{2}}-\tan \left(-90°+PM-\angle \mathrm{P}\left({\mathrm{j}\mathsf{\omega }}_{\mathrm{c}}\right)\right)=0\\ \frac{{\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\lambda }-1}X}{Y}+{\frac{\mathrm{d}\angle \mathrm{P}\left(\mathrm{j}\mathsf{\omega }\right)}{\mathrm{d}\mathsf{\omega }}|}_{{\mathsf{\omega }=\mathsf{\omega }}_{\mathrm{c}}}=0\end{array}$$

(5)

with

$${\mathrm{X}=\mathrm{k}}_{\mathrm{i}}\mathsf{\lambda }\sin \frac{\mathsf{\pi }\mathsf{\lambda }}{2}{+\mathrm{k}}_{\mathrm{d}}{\mathrm{k}}_{\mathrm{i}}{\mathsf{\lambda }\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\mu }}\sin \frac{\mathsf{\pi }\left(\mathsf{\lambda }+\mathsf{\mu }\right)}{2}{+\mathrm{k}}_{\mathrm{d}}{\mathrm{k}}_{\mathrm{i}}{\mathsf{\mu }\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\mu }}\sin \frac{\mathsf{\pi }\left(\mathsf{\lambda }+\mathsf{\mu }\right)}{2}{+\mathrm{k}}_{\mathrm{d}}{\mathsf{\mu }\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\lambda }+\mathsf{\mu }}\sin \frac{\mathsf{\pi }\mathsf{\mu }}{2}$$

and

$${\mathrm{Y}=\mathsf{\omega }}_{\mathrm{c}}{}^{2\mathsf{\lambda }}{+\mathrm{k}}_{\mathrm{i}}^2{+\mathrm{k}}_{\mathrm{d}}^2{\mathsf{\omega }}_{\mathrm{c}}{}^{2(\mathsf{\lambda }+\mathsf{\mu })}{+2\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\lambda }}\cos \frac{\mathsf{\pi }\mathsf{\lambda }}{2}{+2\mathrm{k}}_{\mathrm{d}}{\mathsf{\omega }}_{\mathrm{c}}{}^{2\mathsf{\lambda }+\mathsf{\mu }}\cos \frac{\mathsf{\pi }\mathsf{\mu }}{2}{+2\mathrm{k}}_{\mathrm{d}}{\mathrm{k}}_{\mathrm{i}}{\mathsf{\omega }}_{\mathrm{c}}{}^{\mathsf{\lambda }+\mathsf{\mu }}\cos \frac{\mathsf{\pi }\left(\mathsf{\lambda }+\mathsf{\mu }\right)}{2}$$

$\mathrm{P}\left({\mathrm{j}\mathsf{\omega }}_{\mathrm{c}}\right)$ is the process transfer function at the imposed gain crossover frequency ${\mathsf{\omega }}_{\mathrm{c}}$ and PM is the desired phase margin.

In the maintenance phase, keeping the patient’s depth of hypnosis in the recommended range of 40% to 60% of BIS prioritizes the speed of the controller, the attenuation of the potential disturbances and keeping the patients’ variability minimal. Three main factors compromise the state of well-being of individuals: the intra-patient variability, the feedback signal corrupted by noise, and the surgical stimuli. In this case, to achieve the best control over the system, each term in the FOPID transfer function is tuned based on the minimization of the integral of absolute error (IAE). Unlike in the induction phase, where the controller design was based on an imposed phase margin and gain crossover frequency, in this case, an optimization problem is defined and solved to determine the FOPID parameters. Consider the following optimization problem:

$$\min \left(\mathrm{IAE}\right),°\mathrm{IAE}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}\left|\mathrm{e}\left(\mathrm{t}\right)\right|\mathrm{dt}$$

(6)

where $\mathrm{e}\left(\mathrm{t}\right)=\mathrm{r}\left(\mathrm{t}\right)-\mathrm{BIS}\left(\mathrm{t}\right)$ is the error signal and r(t) is the reference value, usually taken as 50%.

However, in addition to the IAE performance indicator, other parameters are checked based on special functions such as the sensitivity and complementary sensitivity function. In this regard, for the noise-corrupted signal, the complementary sensitivity function described in

$${\left|\mathrm{T}\left(\mathrm{j}\mathsf{\omega }\right)=\frac{{\mathrm{H}}_{\mathrm{open}-\mathrm{loop}}\left(\mathrm{j}\mathsf{\omega }\right)}{{1+\mathrm{H}}_{\mathrm{open}-\mathrm{loop}}\left(\mathrm{j}\mathsf{\omega }\right)}\right|}_{\mathrm{dB}}\le { \mathrm{A}}_{\mathrm{dB}}$$

(7)

can be considered as a performance indicator regarding the attenuation of high-frequency noises, where ${\mathrm{H}}_{\mathrm{open}-\mathrm{loop}}$ denotes the loop transfer function.

Regarding the disturbance rejection performance, the sensitivity function described below

$${\left|\mathrm{S}\left(\mathrm{j}\mathsf{\omega }\right)=\frac{1}{{1+\mathrm{H}}_{\mathrm{open}-\mathrm{loop}}\left(\mathrm{j}\mathsf{\omega }\right)}\right|}_{\mathrm{dB}}\le { \mathrm{B}}_{\mathrm{dB}}$$

(8)

is a relevant component in quantifying the ability of the designed FOPID to reject or diminish the effect of load disturbances. The two performance specifications in (7) and (8) are used as design constraints in the overall optimization routine based on (6). Handling surgical stimuli is significantly more complex compared to the induction-phase specifications. Studies are still researching the perception of pain in humans [11]. For this manuscript, a surgical stimulus as indicated in Figure 2 is used.

Equally important in terms of the obtained closed-loop performance is the ability to ensure a smooth transition of the control system from an induction to a maintenance phase. Switching of two or more controllers implies the synchronization of the command signals so that the bumpless transfer can be achieved. Let ${\mathsf{\epsilon }}_1,{\mathsf{\epsilon }}_2$ be the errors of the induction ${\mathrm{H}}_{{\mathrm{C}}_{\mathrm{I}\mathrm{B}\mathrm{I}\mathrm{S}}}\left(\mathrm{s}\right)$- and maintenance ${\mathrm{H}}_{{\mathrm{C}}_{\mathrm{M}\mathrm{B}\mathrm{I}\mathrm{S}}}\left(\mathrm{s}\right)$-phase controllers, respectively, denoted here by $\mathrm{f}\left(\mathrm{s}\right),\mathrm{g}\left(\mathrm{s}\right)$ having their associated control signals ${\mathrm{c}}_1,{\mathrm{c}}_2$. We simply define the synchronization problem, that is, in the induction phase, the maintenance controller must follow the induction control signal. Therefore, a function $\mathrm{h}\left(\mathrm{s}\right)$ must be found such that the following equality is satisfied:

$${\mathsf{\epsilon }}_1 · \mathrm{h}\left(\mathrm{s}\right){=\mathsf{\epsilon }}_2$$

(9)

Thus, we have

$$\{\begin{array}{c}{\mathsf{\epsilon }}_1 · \mathrm{f}\left(\mathrm{s}\right){=\mathrm{c}}_1\\ {\mathsf{\epsilon }}_2 · \mathrm{g}\left(\mathrm{s}\right){=\mathrm{c}}_1\end{array}$$

(10)

and their division gives

$$\frac{{\mathsf{\epsilon }}_1}{{\mathsf{\epsilon }}_2}·\frac{\mathrm{f}\left(\mathrm{s}\right)}{\mathrm{g}\left(\mathrm{s}\right)}=1$$

(11)

And, equivalently,

$${\mathsf{\epsilon }}_1·\frac{\mathrm{f}\left(\mathrm{s}\right)}{\mathrm{g}\left(\mathrm{s}\right)}{=\mathsf{\epsilon }}_2\Rightarrow \mathrm{h}\left(\mathrm{s}\right)=\mathrm{f}\left(\mathrm{s}\right) ·{ \mathrm{g}}^{-1}\left(\mathrm{s}\right){=\mathrm{H}}_{{\mathrm{C}}_{\mathrm{IBIS}}}\left(\mathrm{s}\right) ·{ \mathrm{H}}_{{\mathrm{C}}_{\mathrm{MBIS}}}^{-1}$$

(12)

where h(s) denotes the fractional-order switching function.

Although for the standard PID controller, the inversion of its transfer function would render (12) impossible, for the FOPID controller, this poses no mathematical problem. The inverse of the maintenance controller’s transfer function is achievable in this particular context, as the approximation of fractional-order controllers given by

$${\mathrm{s}}^{\mathsf{\mu }}=\mathrm{K}{\displaystyle \prod }_{\mathrm{k}=-\mathrm{N}}^{\mathrm{N}}\frac{{\mathrm{s}+\mathsf{\omega }}_{\mathrm{k}}^{\prime }}{{\mathrm{s}+\mathsf{\omega }}_{\mathrm{k}}}$$

(13)

known as Oustaloup’s Recursive Approximation [48] features the necessary properties to invert the transfer function, i.e., the transfer function is always causal and stable, having all the poles inside the unit circle even when inverting the transfer function.

3. Results

In this section, the design of the fractional-order controllers for BIS, RASS and NMB is presented, along with that for the hemodynamic variables CO and MAP. As it has been previously stated, this manuscript presents the first fractional-order control strategy for the combined anesthesia–hemodynamic system. A block diagram representation of the control system is depicted in Figure 3.

The design of the fractional-order controllers indicated in (14)–(18) has been achieved by solving the system of nonlinear equations in (5) with the Matlab “fsolve” function for (14), (17) and (18) and with the graphical method according to [49] for (15) and (16). For the fractional-order controllers used in the maintenance phase, as indicated in Section 3.2, the design is based on a Matlab optimization routine, where the IAE cost function is minimized.

The neuromuscular blockade is completely decoupled from the BIS and RASS signals, which allows for a simple tuning of the FO-PI in this case. The controller for the NMB signal was tuned to meet a phase margin PM = 85°, a gain crossover frequency of ${\mathsf{\omega }}_{\mathrm{c}}=0.01$ rad/s, as well as the iso-damping property [47]. The system of nonlinear equations in (5) was solved, with k_d = 0 and $\mathsf{\mu }$ = 0, yielding the following transfer function of the FO-PI controller:

$${\mathrm{H}}_{\mathrm{NMB}}=0.08\left(1+\frac{0.0227}{{\mathrm{s}}^{1.05}}\right)$$

(14)

To further simplify the overall control strategy, the anesthesia and hemodynamic systems are considered separately. Thus, for the control of the hemodynamic variables, cardiac output CO and mean arterial pressure MAP, a decentralized fractional-order control strategy is used [47]. Two FO-PI controllers are preferred in this case as well. The controller parameters are determined by solving the system in (5), considering again k_d = 0 and $\mathsf{\mu }$= 0. For the CO controller, the following design criteria were imposed: PM = 65°, ${\mathsf{\omega }}_{\mathrm{c}}=0.005$ rad/s and the iso-damping property. For the MAP controller, a PM = 65°, ${\mathsf{\omega }}_{\mathrm{c}}=0.012$ rad/s and the iso-damping property were imposed. Solving (5) using the method proposed in [49] leads to the controller transfer functions as given below:

$${\mathrm{H}}_{\mathrm{CO}}\left(\mathrm{s}\right)=0.38\left(1+\frac{0.0012}{{\mathrm{s}}^{1.2}}\right)$$

(15)

$${\mathrm{H}}_{\mathrm{MAP}}\left(\mathrm{s}\right)=0.07\left(1+\frac{0.0034}{{\mathrm{s}}^{1.27}}\right).$$

(16)

The novelty of this manuscript consists in the double-FOPID design for the induction and maintenance phase of the bispectral index, along with a simple yet efficient synchronization function that ensures a bumpless transfer between the two phases. Since there is a strong interaction between the drugs used to control BIS and RASS, the design of these controllers is considered in a unitary framework. The design procedure is presented hereafter.

3.1. Design of the FOPIDs for the Induction Phase

Focusing on the induction controller, the design steps are considered based on the level of interactions between agents. Since RASS is affected solely by changes in remifentanil, this loop is firstly tuned, followed by the tuning of the BIS controller, as the BIS signal is greatly influenced by the combination of propofol and remifentanil. Empirical knowledge leads to the observation that faster settling times of RASS results in the best behavior for the combined control of BIS and RASS in terms of undershoot and time to target. It is important to mention that (1) describes the general fractional-order PID equation. Since, for the tuning of the induction controllers, only three performance specifications were used and the resulting set of equations is (5), the following simplification of the FOPID transfer function was considered: $\mathsf{\lambda }=\mathsf{\mu }$ and T_i = 4T_d [50]. To tune the parameters, for both the RASS and BIS controllers, the iso-damping requirement was considered along with an imposed phase margin of ${\mathsf{\phi }}_{\mathrm{m}}=78°$ and gain crossover frequency ${\mathsf{\omega }}_{\mathrm{c}}=0.021 \mathrm{rad}/\mathrm{s}$ for the RASS controller, and a ${\mathsf{\phi }}_{\mathrm{m}}{=84° \mathrm{and} \mathsf{\omega }}_{\mathrm{c}}=0.017\mathrm{rad}/\mathrm{s}$ for the BIS controller. The resulting controllers are shown below:

$${\mathrm{H}}_{\mathrm{BIS}}\left(\mathrm{s}\right)=0.0033\left(1+\frac{0.0537}{{\mathrm{s}}^{0.97}}+0{.0134\mathrm{s}}^{0.97}\right)$$

(17)

$${\mathrm{H}}_{\mathrm{RASS}}\left(\mathrm{s}\right)=0.0753\left(1+\frac{0.1553}{{\mathrm{s}}^{0.94}}+0{.0388\mathrm{s}}^{0.94}\right)$$

(18)

The robustness of the designed controllers must be validated by considering the inter-patient variability. Thus, the same tests are performed for the nominal patient 1 in Table 1, as well as for the remaining 23 patients from Table 1. The same controllers are used both for the induction and maintenance phases. The closed-loop performance is indicated in Figure 4 for the output signal BIS and RASS and Figure 5 for the corresponding input drugs, propofol and remifentanil. Figure 6 presents a detailed response of the BIS signal, with a focus on the induction phase. Performances are satisfactory and balanced from speed to undershoot as specified in

Table 2. Performance indexes for each patient’s response.

Patient	TT (s)	Undershoot (%)	BIS-NADIR
1	140.56	1.78	40.82
2	104.41	1.93	40.77
3	219.36	1.63	40.86
4	166.37	1.64	40.85
5	112.25	2.68	40.83
6	129.51	1.66	40.83
7	96.73	2.83	40.74
8	117.32	1.95	40.79
9	120.20	1.83	40.81
10	131.79	1.86	40.81
11	129.36	1.78	40.81
12	145.74	1.66	40.85
13	124.43	1.77	40.82
14	126.91	1.69	40.83
15	124.84	1.91	40.80
16	230.84	1.03	40.95
17	134.84	1.76	40.82
18	142.40	1.67	40.84
19	114.15	1.93	40.79
20	122.36	1.86	40.80
21	124.04	1.75	40.81
22	124.34	1.73	40.82
23	109.79	1.94	40.78
24	152.03	1.56	40.86
Mean	135.19	1.83	40.82
Std. Dev.	31.65	0.34	0.04
Min	96.73	1.03	40.74
Max	230.84	2.83	40.95

, having the control signals within their limits as shown in Figure 5. The time to target ranges from 95 to 230 seconds, with most patients averaging around 135 seconds. Inter-patient variability is significantly reduced compared to the open-loop responses. Further, the maximum undershoot is 2.8% and no oscillations can be observed on the output.

3.2. Design of the FOPIDs for the Maintenance Phase

The design of the FOPIDs in this section is performed via an optimization routine considering both FOPID and FO-PI control structures and studies the optimal maintenance performances as well as the relevance of the iso-damping property in the maintenance phase.

The FO-PI ensuring the iso-damping property optimizes the IAE index in (6) and has the constraints as stated in (5) reduced to the specific FO-PI parameters k_p, k_i, μ. The performance specifications involve a phase margin of 60 degrees, gain crossover frequency of 0.0305 rad/s and the iso-damping property. Two other controllers, namely the FO-PI and the FOPID, optimize the IAE index in (6) and have the constraints as stated in (7) and (8). The attenuation parameters A_dB, B_dB were set to be as high as possible in the range of 0.05–0.15 rad/s, as this is the average span of the surgical stimuli. The Matlab function “fmincon” has been used to implement the optimization routine for each controller, yielding the control parameters in (19), (20) and (21).

The resulting controllers are shown below:

$${\mathrm{H}}_{\mathrm{FOPI}-\mathrm{ISO}}\left(\mathrm{s}\right)=0.0072\left(1+\frac{0.0188}{{\mathrm{s}}^{1.24}}\right)$$

(19)

$${\mathrm{H}}_{\mathrm{FOPI}-\mathrm{IAE}}\left(\mathrm{s}\right)=0.0053\left(1+\frac{0.0532}{{\mathrm{s}}^{1.04}}\right)$$

(20)

$${\mathrm{H}}_{\mathrm{FOPID}-\mathrm{IAE}}\left(\mathrm{s}\right)=0.0065\left(1+\frac{0.0549}{{\mathrm{s}}^{1.02}}+4{.99\mathrm{s}}^{1.12}\right)$$

(21)

The comparison between the BIS controller previously used in the induction phase described in (17) and the specifically optimized controllers for the maintenance phase in (19)–(21) is shown in Figure 7. The surgical stimulus profile used as the disturbance is shown in Figure 2.

Due to synergistic effects, each controller struggles in some manner to counteract the sudden changes, undershooting at the time the peak surgical stimulus happens (around t = 550 s). In the case of the FO-PI controller with the iso-damping property, the speed performances are decent, but this comes with the cost of an increased oscillation tendency, as well as a larger IAE of 45.38. A moderate number of oscillations are present in the other two cases, especially for the FOPID case. These two controllers exhibit a better IAE score of IAE = 30.58 for the FO-PI controller and an IAE = 33.23 for the FOPID. Despite setting the attenuation thresholds A_dB, B_dB high in the range of the occurrence of the stimuli, the disturbance attenuation B_dB failed to satisfy the constraint, as the nominal sensitivity peaks are positive values and are present in the range of our interest considering the span of surgical stimuli ranging from 40 to 120 seconds (0.05–0.15 rad/s). However, their values have still been minimized to reduce their amplifying influence. The corresponding error signals associated to the results in Figure 7 are included in Figure 8. This clearly shows that the smallest error is achieved using the proposed FOPID controller. The numerical performance results are summarized in

Table 3. Numerical performance comparison of the optimized controllers.

Controller Type	Iso-Damping	Phase Margin PM (deg)	Gain Crossover Frequency ${\mathsf{\omega }}_{\mathrm{c}}$ (rad/s)	Nominal Sensitivity Peak (dB)	T Function Value at ${\mathsf{\omega }}_{\mathsf{t}}={10}^2\left(\mathbf{dB}\right)$	Integral of Absolute Error (IAE) (-)
Induction FOPID	YES	84	0.017	+3.28	−75.0	66.31
FO-PI with iso-damping	YES	60	0.0305	+9.06	−65.5	45.38
FO-PI with minimum IAE	NO	61.39	0.0308	+2.7	−101.1	33.23
FOPID with minimum IAE	NO	55.73	0.0400	+5.23	−56.2	30.58

below.

Looking at the sensitivity function of these two best performing controllers reveals that the nominal sensitivity peak is significantly larger in the FOPID case (Figure 9), leading to the conclusion that the FO-PI is less sensitive to disturbances as shown in Figure 10. The T function describing the sensitivity to noise suggests differences as well, but satisfactory performances from both are obtained, as noises with a frequency larger than ${\mathsf{\omega }}_{\mathrm{t}}=100 \mathrm{rad}/\mathrm{s}$are being significantly attenuated.

Therefore, the controller in (20) has been tested on all 24 patients. The closed-loop simulation results from Figure 11 reveal satisfactory performances in the maintenance phase. From the point of view of intra- and inter-patient variability, all patients reacted similarly to the input drugs. Additionally, pronounced oscillations are absent and the surgical stimuli is rejected successfully in every case. Detailed performance results are shown in

Table 4. Disturbance rejection performances for each patient.

Patient	TTn (seconds)	BIS-NADIRn (%)	TTp (seconds)	BIS-NADIRp (%)
1	35.11	41.21	0	54.46
2	36.95	41.02	0	54.67
3	32.05	41.28	0	54.46
4	34.69	41.27	0	54.44
5	36.79	41.20	0	54.57
6	35.69	41.26	0	54.43
7	40.24	40.93	4.56	55.30
8	35.95	41.12	0	54.46
9	35.99	41.17	0	54.44
10	35.36	41.18	0	54.46
11	35.55	41.19	0	54.45
12	35.11	41.29	0	54.43
13	35.86	41.21	0	54.44
14	35.83	41.27	0	54.43
15	35.61	41.15	0	54.46
16	31.08	41.43	0	54.41
17	35.32	41.20	0	54.45
18	35.19	41.28	0	54.43
19	36.17	41.10	0	54.49
20	35.81	41.16	0	54.45
21	35.81	41.17	0	54.44
22	35.90	41.22	0	54.43
23	36.44	41.07	0	54.59
24	34.98	41.33	0	54.43
Mean	35.56	41.20	0.19	54.50
Std. Dev.	1.63	0.10	0.93	0.18
Min	31.08	40.93	0	54.41
Max	40.24	41.43	4.56	55.30

.

3.3. Design of the Synchronization Function for Bumpless Transfer between the Induction- and Maintenance-Phase Controllers

The implementation of the synchronization function h(s) that ensures a bumpless transfer between induction- and maintenance-phase controllers is conducted as shown in Figure 12, depicting the function h(s), the two controllers and the switching elements.

The comparator “<=50”, featuring the logic “OR” block and a memory block, controls the switching logic. The blocks are configured such that once the induction controller reaches the 50% feedback mark, the switching of controllers occurs. This step only occurs once. The synchronized control signals, therefore, are shown in Figure 13.

Lastly, the results of synchronizing the induction-phase controller and the maintenance controller are showcased in Figure 14 and Figure 15, with their associated control signals in Figure 16 and Figure 17.

Table 5. Final induction-phase performances for each patient.

Patient	TT (s)	Undershoot (%)	BIS-NADIR
1	140.56	1.37	41.21
2	104.41	1.26	41.02
3	219.36	1.41	41.28
4	166.37	1.35	41.27
5	112.25	2.53	41.20
6	129.51	1.24	41.26
7	96.73	2.57	40.93
8	117.32	1.35	41.12
9	120.20	1.29	41.17
10	131.79	1.38	41.18
11	129.36	1.32	41.19
12	145.74	1.33	41.29
13	124.43	1.28	41.21
14	126.91	1.24	41.27
15	124.84	1.37	41.15
16	230.84	1.00	41.43
17	134.84	1.33	41.20
18	142.40	1.32	41.28
19	114.15	1.31	41.10
20	122.36	1.33	41.16
21	124.04	1.26	41.17
22	124.34	1.25	41.22
23	109.79	1.28	41.07
24	152.03	1.29	41.33
Mean	135.19	1.40	41.20
Std. Dev.	31.65	0.36	0.10
Min	96.73	1.00	40.93
Max	230.84	2.57	41.43

presents the final induction-phase performances and

Table 6. Final disturbance rejection performances for each patient.

Patient	TTn (s)	BIS-NADIRn (%)	TTp (s)	BIS-NADIRp (%)
1	35.11	41.21	0	54.46
2	36.95	41.02	0	54.67
3	32.05	41.28	0	54.46
4	34.69	41.27	0	54.44
5	36.79	41.20	0	54.57
6	35.69	41.26	0	54.43
7	40.24	40.93	4.56	55.30
8	35.95	41.12	0	54.46
9	35.99	41.17	0	54.44
10	35.36	41.18	0	54.46
11	35.55	41.19	0	54.45
12	35.11	41.29	0	54.43
13	35.86	41.21	0	54.44
14	35.83	41.27	0	54.43
15	35.61	41.15	0	54.46
16	31.08	41.43	0	54.41
17	35.32	41.20	0	54.45
18	35.19	41.28	0	54.43
19	36.17	41.10	0	54.49
20	35.81	41.16	0	54.45
21	35.81	41.17	0	54.44
22	35.90	41.22	0	54.43
23	36.44	41.07	0	54.59
24	34.98	41.33	0	54.43
Mean	35.56	41.20	0.19	54.50
Std. Dev.	1.63	0.10	0.93	0.18
Min	31.08	40.93	0	54.41
Max	40.24	41.43	4.56	55.30

shows the final maintenance-phase performances in the synchronized context. Further improvements have been achieved due to the use of dedicated fractional-order controllers (18) and (20) for the induction and maintenance phases, synchronized using the bumpless transfer function procedure. For instance, a slight improvement in the maximum undershoot can be observed (2.5%) compared to the last scenario (3%) having the same controllers. The time to target remained unchanged, satisfying the specifications (95 to 230 seconds). Redundant oscillations are absent and the control signals in Figure 16 are within safety limits. From the maintenance phase’s point of view, most surgical stimuli are rejected in time for all patients without being affected by intra-patient changes. Sudden nociception stimuli rejection (700 to 750 s) and large bolus-type stimuli rejection (540 to 560 s) are successful in this context, considering that BIS is maintained within safe limits (40% to 60%).

4. Overdose Analysis and Controller Validation for Variable Measurement Delays

The maximum allowable values of medical substances in general anesthesia for patients under 55 years of age are shown in

Table 7. Maximum dosage of anesthetic substances.

Substance	Maximum Allowed Dosage for Adult Patients up to 55 Years Old
	Induction	Maintenance
Propofol	2.5 mg/kg IV titrated to 40 mg every 10 seconds until onset of induction	150 to 200 ug/kg/min IV for the first 10 to 15 min, then decreased by 30% to 50% during the first 30 min of maintenance
Remifentanil	1 ug/kg/min for 30 to 60 seconds	0.4 ug/kg/min (range 0.1 to 2 micrograms/kg/min)
Dopamine	10 ug/kg/min	50 ug/kg/min
SNP	0.3 ug/kg/min, but the dose can be adjusted up to 10 ug/kg/min in major cases	0.3 ug/kg/min, but the dose can be adjusted up to 10 ug/kg/min in major cases
Atracurium	0.6 mg/kg (bolus injection)	13 ug/kg/min (can be adjusted up to 29.5 ug/kg/min if patient does not respond to initial dose—major cases only)

. Only intravenous (IV) administration using volumetric pumps are considered, except for atracurium, which is injected in a single dose in the induction phase. Values represent maximum doses according to existing real-life surgical protocols, considering all synergic and antagonistic substance interactions, according to IBM^® Watson^® Micromedex™, Cerner Multum INC and the American Society of Health-System Pharmacists (ASHP).

Any delay in measuring BIS values has a negative effect on controller performance. For very long delays, the controller may even lose feasibility (and stability). In real-life BIS measurements, there are various delays associated with BIS measurements due to the filtering algorithms used by the monitors, as well as in adapting the pre-processing steps of artifact rejection [51,52]. This section aims at investigating the effect of delays on controller performance for a wide delay range: between 1 to 60 seconds. According to [53], the 95% confidence interval for delays is 12.7–27.6 seconds, with a nominal value of 19.7 seconds, as presented in Figure 1.

Remifentanil doses are unaffected by the variation in dead time, being constant in the maintenance phase at approx. 1.07 ug/kg/min, as illustrated in Figure 18 The absence of overdose can be concluded from Figure 19 and Figure 20, where the anesthetic agents are below the maximum doses shown in Table 7.

Lower dead time values cause a slower BIS response, because the propofol control signal is also slower. Intermediate values show a high induction time for BIS; however, the trigger time is affected because the delay is greater than the 30 second margin. In this case, the response begins to oscillate as it approaches the stability limit (for 50 to 60 seconds dead time).

A further investigation of the maximum dead-time tolerance of the controllers for all 24 patients concludes that some patients have a higher tolerance. However, most of them converge to 60 seconds of maximum dead time, as shown in Figure 20, where the dead time is set to 60 seconds for all patients.

For small delays, responses are, on average, slower than for a constant average delay of 19.7 seconds. However, the risk of overdose is significantly reduced. The test case illustrated in Figure 21 represents the responses of all 24 patients who have an almost instantaneous response time to anesthetic agents for a fixed 1-second delay.

5. Conclusions

Automatically controlled general anesthesia comes as an assistant tool to anesthesiologists so that safe surgical operations can be offered for every patient in need. In this study, several specifications were considered in elaborating the final result, such as strictly imposed time in driving every patient from a fully awake state to a moderate hypnotic state, eliminating the concerns regarding the inter- and intra-patient variability by imposing the iso-damping property and avoiding excessive undershoot in the induction phase for all patients, and handling surgical stimuli in the maintenance phase such that the hypnotic states of the patients are minimally affected. The design attempt of a fractional-order PI results in decent performances in the induction phase as described in [3], but not in the maintenance phase, as these two phases represent two drastically different operating points. Regarding the control design, the stability of the overall system was acceptable around a range of parameters; however, the working range of the imposed gain crossover frequency and phase margin was extremely limited. Crossing this range often resulted in an unstable system. Additionally, the fractional elements λ, μ in this context happened to be considerably sensitive, perhaps due to the high non-linearity and the synergistic nature of the system. Despite the severe limiting factors, a combination of a fractional-order PID controller and a fractional-order PI controller, specifically tuned for induction-phase specifications and maintenance-phase specifications, respectively, has been successfully validated on all patients within satisfactory performances. A bumpless transfer procedure was considered in the synchronization process of the two controllers acting on the stages of general anesthesia. A final overdose risk analysis has been performed, taking into consideration variable measurement delays, showing the efficiency of the controllers for delays up to 60 s.

Remark: although designed for the specific BIS control system, the theoretical design of the fractional-order function ensuring bumpless transfer can be extended to various other applications.

Further research includes the design of integer-order PIDs and comparisons with the proposed fractional-order control system presented in this paper.