Complex-order PID controller design for enhanced blood-glucose regulation in Type-I diabetes patients

Abstract

Type-I Diabetes (TID) is a chronic autoimmune disease that elevates the glucose levels in the patient’s bloodstream. This paper formulates a fractional complex-order Proportional-Integral-Derivative (PID) control strategy for robust Blood Glucose (BG) regulation in TID patients. The glucose-insulin dynamics in blood plasma are modeled via the Bergman-Minimal-Model. The proposed control procedure employs the ubiquitous fractional order PID controller as the baseline BG regulator. The design flexibility of the baseline regulator to effectively normalize the BG levels is enhanced by assigning complex orders to the integral and differential operators instead. The resulting Complex Order PID (CO-PID) regulator strengthens the controller’s robustness against abrupt variations in the patient’s BG levels caused by meal disturbances or sensor noise. The controller parameters are numerically optimized offline. The aforesaid propositions are justified by performing credible simulations in which the proposed controller is tasked to effectively track a set point value of 80 mg/dL from an initial state of hyperglycemia under various disturbance factors. As compared to the FO-PID controller, the CO-PID controller improves the reference tracking-error, transient recovery-time, and control expenditure by 13.1%, 33.4%, and 28.1%, respectively. The simulation results validate the superior reference-tracking accuracy of the proposed CO-PID controller for BG regulation.

Introduction

Type-I Diabetes (TID) is a chronic organ-specific disease that destroys the insulin-generating pancreatic T-cells, thus, preventing them to generate sufficient insulin in the blood plasma. ¹ This inevitably creates an impairment in the generation and utilization of insulin, which elevates the glucose levels in the patient’s bloodstream. ² The insufficient insulin production, or hyperglycemia, in the patient’s body eventually leads to a life-long reliance on exogenous insulin therapy for the maintenance of normal glycemic control. ³ The closed loop glucose-insulin regulation systems, as illustrated in Figure 1, are used by TID patients to normalize the blood glucose (BG) levels. The closed-loop system contains a sensor that continually measures the BG levels and feeds them to a feedback control law that operates the insulin pump to appropriately alter the insulin infusion rate (IIR) in patient’s bloodstream. ⁴ Designing agile and robust glycemic control laws for TID patients is a complex control engineering problem.

Figure 1.

Closed-loop glucose-insulin regulation mechanism.

Extensive research has been conducted by scientists and researchers to devise reliable BG regulation control strategies to improve the external insulin administration process for the reliable glycemic control of TID patients.^5–7 The PID controller is a computationally simple model-free control strategy that uses the system’s state error dynamics to yield a reliable control effort. ⁸ However, it lacks the degree of freedom to effectively compensate for the non-vanishing disturbances and random noise. ⁹ Other notable control schemes used for robust BG regulatory control are sliding-mode controllers,^10,11 linear quadratic regulators, ¹² neural controllers, ¹³ etc. Each of these control schemes tends to make a trade-off between the controller’s computational complexity, robustness, time optimality, and control efficiency. Recently, a novel deep learned type-II fuzzy logic glycemic control approach for type-I diabetes patients has been proposed in Mosavi et al. ¹⁴ The said scheme employs the fuzzy immersion and invariance theorem to establish online adaptive tuning rules for unknown parameters that harness the controller’s ability to robustly compensate for the estimations errors and regulate the BG levels in a short time. However, the aforementioned scheme requires specialized computational resources for its execution.

The PID controller’s design flexibility is generally improved by restructuring the control law with the aid of fractional calculus. ¹⁵ This objective is achieved by replacing the PID controller’s integer-order integral and differential operators with fractional-order counterparts. ¹⁶ The fractional-order controller optimizes the control resource application to effectively reject external disturbances, chaotic behavior, and intrinsic nonlinearities. ¹⁷ Although, this scheme optimizes the controller’s design flexibility and increases its degrees of freedom to better address the exogenous disturbances; however, it also introduces two new parameters in the control law to be pre-specified by the control expert. The fractional order controllers have yielded promising results in regulating the BG levels in TID patients. ¹⁸

The ability of the conventional fractional order controller to address the nonlinear behavior of chaotic systems can be harnessed by extending it to complex controllers. ¹⁹ The complex controllers are realized by replacing the integer-order integral and differential operators with the complex-order counterparts of the form $z=a+\mathrm{jb}$ , where $a$ and $b$ are the real and imaginary parts of the complex number $z$ . ²⁰ The studies have shown that the assignment of complex orders to the aforementioned operators increases the controller’s degrees of freedom, which enhances its resilience against parametric uncertainties especially when the gain and phase characteristics are notable. ²¹ The complex-order controllers are extensively used to solve practical control engineering problems owing to their aforementioned attributes. ²²

The main contribution of this article is the systematic formulation and verification of the Complex-Order PID (CO-PID) controller for the normalization and robust regulation of BG levels in TID patients. The Bergman Minimal Model (BMM) is commonly used to dynamically model the glucose-insulin interaction in patients.^10,15 The FO-PID controller is used as the baseline BG regulator, since it has been observed to effectively normalize the BG levels in TID patients.^15,17 This observation sets a clear motivation to systematically evolve the FO-PID controller into the proposed CO-PID control law, an improved variant, to further robustify the system’s glycemic disturbance attenuation capacity. Moreover, the CO-PID scheme is relatively more flexible to deal with inherently chaotic and complex physical systems. Since the BG regulation problem manifests the aforementioned properties; therefore, the usage of CO-PID control law in this research is warranted. The innovative contributions of this article are thus presented as follows:

1. Formulation and calibration of a CO-PID control law for the BG regulation and normalization. This is done by augmenting the integral and differential operators of the baseline controller with complex orders instead.

2. Verification of the proposed control law via credible simulations to normalize the BG levels of a patient (modeled via BMM) to a set point of 80 mg/dL from an initial state of hyperglycemia, under the influence of meal disturbance and sensor noise.

The formulation of the CO-PID controllers introduces multiple parameters in the control law that significantly increase the controller’s degrees of freedom, which in turn enables it to better address the parametric disturbances and effectively normalize the BG levels. All controller parameters are numerically selected offline by using a well-established numerical optimization method. The CO-PID control law can be easily realized by using modern embedded processors. The computational realization as well as the simulated testing of the proposed BG regulator is done in the MATLAB SIMULINK environment. The design and verification of the proposed CO-PID glycemic control law, to optimize the BG regulation in TID patients under exogenous disturbances, has never been attempted previously in the scientific literature available. Thus, this article primarily focuses on this innovative idea.

The remaining paper is organized as follows. The dynamics of the glucose-insulin regulation mechanism and the associated baseline FO-PID control law are discussed in Section 2. The formulation of the proposed CO-PID control law is presented in Section 3. The parameter tuning procedure is discussed in Section 4. The simulation results are analyzed in Section 5. Finally, the article is concluded in Section 6.

Glucose-insulin regulation mechanism

This article synthesizes an optimal insulin delivery system for diabetics to effectively reduce their high BG levels. ¹⁰ The BG levels within the patient’s body tend to get disturbed by the excessive or insufficient infusion of insulin, which adversely affects the patient’s health. The closed-loop glucose-insulin regulation system is an automated BG control scheme that administers the precise amount of glucose. The subcutaneous injection technique imitates the biological insulin discharge, and thus, it is regarded as the most secure method of administering an insulin dose to the body. The close-loop insulin infusion technique is already depicted in Figure 1. The system continually monitors the BG level during insulin infusion via a dedicated sensor. The sensor feedback is compared with the reference set-point. The corresponding deviations in the desired and actual BG levels are fed to a feedback controller that appropriately manipulates the exogenous insulin infusion rate to normalize the BG levels. The controller stops the infusion when the actual and desired BG levels become equal.

Mathematical model for glucose-insulin dynamics

The BMM is the most commonly used mathematical model for representing the kinetics of glucose-insulin interaction. ¹⁵ The model offers minimal biological complexities which makes it a viable option for evaluating the effectiveness of the artificial pancreas. The BMM represents the dynamics of glucose and insulin via the following set of first-order differential equations. ²³

$$\begin{array}{c}\stackrel{·}{G}\left(t\right)=-p_{1}G\left(t\right)-X\left(t\right)(G\left(t\right)-G_b)+\frac{G_m\left(t\right)}{V_1}\hfill \\ \stackrel{·}{X}\left(t\right)=-p_{2}X\left(t\right)+p_{3}I\left(t\right)\hfill \\ \stackrel{·}{I}\left(t\right)=-n(I\left(t\right)-I_b)+\frac{u\left(t\right)}{V_1}\hfill \end{array}$$

(1)

where, $G\left(t\right)$ is the BG concentration variable, $X\left(t\right)$ is the insulin concentration in a “remote” compartment variable, $I\left(t\right)$ is the blood-insulin concentration variable, $G_m\left(t\right)$ is the meal disturbance input variable, $u\left(t\right)$ is the manipulated insulin infusion rate (IIR) variable, $G_b$ is the basal BG concentration, $I_b$ is the basal blood-insulin concentration, $n$ is the first-order decay rate of plasma insulin, $V_1$ represents the blood volume, and $p_1$ , $p_2$ , $p_3$ are pre-calibrated parameters that are specific to the blood specimen being used. A linear system is generally represented in state-space as shown below.

$$\stackrel{·}{x}\left(t\right)=\mathit{A}x\left(t\right)+\mathit{B}v\left(t\right),y\left(t\right)=\mathit{C}x\left(t\right)+\mathit{D}v\left(t\right)$$

(2)

where, $x\left(t\right)$ is the state vector, $y\left(t\right)$ is the output vector, $v\left(t\right)$ is the input vector, $\mathit{A}$ is the system matrix, $\mathit{B}$ is the input matrix, $\mathit{C}$ is the output matrix, and $\mathit{D}$ is the feed-forward matrix. The state vector and input vector of the glucose-insulin dynamical system are provided in equation (3). ²³

$$x\left(t\right)={\left[\begin{array}{ccc}G\left(t\right)\hfill & X\left(t\right)\hfill & I\left(t\right)\hfill \end{array}\right]}^T,v\left(t\right)={\left[\begin{array}{cc}{G}_m\left(t\right)\hfill & u\left(t\right)\hfill \end{array}\right]}^T$$

(3)

The nominal state-space model representing the glucose-insulin dynamics are given by equation (4). ²³

$$\begin{array}{c}\mathit{A}=\left[\begin{array}{ccc}-p_1\hfill & G_b\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -p_2\hfill & p_3\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -n\hfill \end{array}\right],\mathit{B}=\left[\begin{array}{cc}0\hfill & \frac{1}{V_1}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{\frac{1}{V_1}}& 0\hfill \end{array}\right],\hfill \\ \mathit{C}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right],\mathit{D}=\left[\begin{array}{cc}0\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

(4)

As per the matrix, the system has only one output variable, $y\left(t\right)=G\left(t\right)$ . The model parameters for a healthy subject and three different patients, used in this research, are listed in Table 1. ¹⁰

Table 1.

BMM parameters for three different patients. ¹⁰

Parameter	Normal	Patient 1	Patient 2	Patient 3	Units
$p_1$	0.0317	0	0	0	$\mathrm{mi}{n}^{-1}$
$p_2$	12.3 × 10−3	20 × 10−3	7.2 × 10−3	14.2 × ′10−3	$\mathrm{mi}{n}^{-1}$
$p_3$	4.92 × 10−6	5.3 × 10−6	2.16 × 10−6	99.4 × ′10−6	$\mathrm{mi}{n}^{-1}$
$n$	0.2659	0.3	0.2465	0.2814	$\mathrm{mi}{n}^{-1}$
$I_b$	7	7	7	7	$\mathrm{mU}/L$
$G_b$	80	80	80	80	$\mathrm{mg}/\mathrm{dL}$
$V_1$	12	12	12	12	$L$
$G\left(0\right)$	200	200	200	200	$\mathrm{mg}/\mathrm{dL}$
$X\left(0\right)$	0	0	0	0	$\mathrm{mi}{n}^{-1}$
$I\left(0\right)$	50	50	50	50	$\mathrm{mU}/L$

Fractional-order PID control

The PID control law is formulated as a weighted sum of the system’s classical state error, error-integral, and error-derivative. ²⁴ The proportional term rejects the instantaneous error, the derivative term speeds up the transient response and predicts the system’s future state variations, and the error-integral term optimizes the system’s reference-tracking behavior and damps the steady-state fluctuations. The conventional integer-order PID control law formulation for BG regulation is expressed below. ²⁵

$$u\left(t\right)=k_{p}e\left(t\right)+k_i\underset{0}{\overset{t}{\int }}e\left(\tau \right)d\tau +k_d\stackrel{·}{e}\left(t\right)$$

(5)

$$\mathrm{where},e\left(t\right)=G_{\mathrm{ref}}-G\left(t\right)$$

where, $G_{\mathrm{ref}}$ and $G\left(t\right)$ represent the desired and actual levels of BG concentration in the patient’s body, respectively, and $e\left(t\right)$ represents the error between the desired and actual BG concentration levels. In this research, the value of $G_{\mathrm{ref}}$ is set at 80 mg/dL.

The conventional PID control law lacks the necessary degrees of freedom to effectively compensate for bounded disturbances in BG regulation systems owing to its integer-order integral and differential operators. Hence, in this work, the PID control law is augmented with fractional calculus to improve its design flexibility.

In fractional calculus, the integral and differential operators involved in the control law are assigned pre-calibrated fractional powers rather than the usual integer ones. ²⁶ The generalized fractional order of the operator is denoted by the symbol $D^{\lambda }$ ; where, $\lambda$ represents the operator’s fractional order. Fractional calculus is frequently chosen to represent and regulate practical engineering systems with extremely chaotic and nonlinear dynamics. The studies show that the innate properties of the fractional-order operators allow the corresponding fractional control law to recognize and effectively reject the exogenous disturbances encountered by the BG regulation systems. ²⁷ The following well-known formulations offered by Riemann-Liouville, Gruunwald-Letnikov, and Caputo ²⁸ serve as the mathematical definitions of fractional operators.

$$D^{\lambda }f\left(t\right)=\frac{1}{\mathrm{\Gamma }(m-\lambda )}\frac{d^m}{d{t}^m}\underset{a}{\overset{t}{\int }}\frac{f\left(\tau \right)}{{(t-\tau )}^{\lambda -m+1}}d\tau$$

(6)

where, $\mathrm{\Gamma }\left(x\right)$ is the Euler gamma function, $m$ is an integer such that $m-1<\lambda <m$ .

$$D^{\lambda }f\left(t\right)=\underset{h\to 0}{lim}\frac{1}{h^i}\sum _{i=0}^{(t-a)/h}{(-1)}^i\left(\begin{array}{c}\lambda \hfill \\ \multicolumn{1}{c}{i}\end{array}\right)f(t-\mathrm{ih})$$

(7)

where, $\left(\begin{array}{c}\lambda \hfill \\ \multicolumn{1}{c}{i}\end{array}\right)=\mathrm{\Gamma }(\lambda +1)/\mathrm{\Gamma }(i+1)\mathrm{\Gamma }(\lambda -i+1)$ , and $h$ is the step size.

$$D^{\lambda }f\left(t\right)=\frac{1}{\mathrm{\Gamma }(\lambda -m)}\underset{a}{\overset{t}{\int }}\frac{f^m\left(\tau \right)}{{(t-\tau )}^{\lambda -m+1}}d\tau$$

(8)

The FO-PID control law is formulated by supplementing the conventional PID controller with predetermined fractional orders integral and differential operators as indicated below. ²⁹

$$u\left(t\right)=k_{p}e\left(t\right)+k_i\left(D^{-\alpha }e\left(t\right)\right)+k_d\left(D^{\gamma }e\left(t\right)\right)$$

(9)

The aforementioned FO-PID control law possesses five unique parameters that enhance the controller’s design flexibility; namely, $k_p$ , $k_i$ , $k_d$ , $\alpha$ , and $\gamma$ . Together these five parameters increase the controller’s flexibility to manipulate the damping control effort of the procedure. The parameters $\alpha$ and $\gamma$ represent the real-numbered fractional orders allocated to the integral and differential operators of $e\left(t\right)$ , respectively. The following expression identifies the transfer function of the FO-PID control law.

$$C\left(s\right)=\frac{U\left(s\right)}{E\left(s\right)}=k_p+\frac{k_i}{s^{\alpha }}+k_d{s}^{\gamma }$$

(10)

where, $s$ represents the Laplace operator. It is quite challenging to computationally realize the operators $s^{\gamma }$ and $s^{\alpha }$ due to their fractional powers. To simplify their implementation, the fractional operators are approximated via the Oustaloup recursive filters. ²⁹ The Oustaloup approach is used to implement the fractional operator as demonstrated below.

$$s^{\lambda }=C\underset{i=1}{\overset{M}{\Pi }}\frac{1+(\mathrm{s}/\_\mathrm{z},\mathrm{i})}{1+(\mathrm{s}/\_\mathrm{p},\mathrm{i})}$$

(11)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{such}\mathrm{that},{\omega }_{z,i}={\omega }_L{(\_\mathrm{H}/\_\mathrm{L})}^{2\mathrm{i}-1-/2\mathrm{M}},}\\ {\omega }_{p,i}={\omega }_L{(\_\mathrm{H}/\_\mathrm{L})}^{2\mathrm{i}-1+/\mathrm{2M}}\hfill \end{array}$$

where, $M$ is the filter’s order, and ${\omega }_L$ and ${\omega }_H$ represent the filter’s lower and the upper translational frequencies, respectively. The value of $C$ is chosen such that ${\left(j\omega \right)}^{\lambda }=1$ at unity frequency. In this research, a fifth-order Oustaloup’s recursive filter is employed with ${\omega }_L={10}^{-5}$ rad/s and ${\omega }_H={10}^5$ rad/s to approximate the fractional operators. The analytical computation of the fractional orders is quite difficult. Therefore, all five parameters of the FO-PID controller are optimized offline by using the tuning procedure discussed in Section 4.

Proposed complex order PID control methodology

The CO-PID controller is an improved variant of the fractional-order controllers that further enhance the robustness and flexibility of the control law by retrofitting the integral and differential operators with complex orders rather than real-numbered integer or fractional orders. The complex order control schemes are derived from the third generation CRONE control. ³⁰ As discussed earlier, the CO-PID control procedure has been proposed in the available literature for DC motor’s speed control, inverted pendulum’s balancing control, and surface roughness control in machining applications.^21,22,31 However, the dynamics of the glucose-insulin interactions in human body are quite different than that of the aforementioned electro-mechanical systems. The CO-PID controller, formulated in this section, introduces additional controller parameters that increases the controller’s degrees of freedom and improves its design flexibility, which enables it to effectively track the reference under exogenous disturbances. Hence, the idea of formulating a CO-PID controller for glycemic control application is being investigated in this article.

The formulation of the CO-PID control law for the BG regulation system is presented in equation (12).

$$u\left(t\right)=k_{p}e\left(t\right)+k_i\left(D^{-f}e\left(t\right)\right)+k_d\left(D^{g}e\left(t\right)\right)$$

(12)

$$\mathrm{such}\mathrm{that},f=\alpha +j\beta ,g=\gamma +j\delta$$

where, $\alpha$ and $\gamma$ are coefficients of the real parts of the complex order of the integral/differential operator, and $\beta$ and $\delta$ are the coefficients of the imaginary parts of the complex order of the integral/differential operator, respectively. As compared to the FO-PID law, this scheme introduces two additional parameters in the control law. Together, these seven parameters ( $k_p$ , $k_i$ , $k_d$ , $\alpha$ , $\gamma$ , $\beta$ , and $\delta$ ) significantly enhance the adaptability of the control law to efficiently reject the disturbances while ensuring a minimum-time transient recovery. The following expression represents the modified transfer function of the control law. ³¹

$$C\left(s\right)=\frac{U\left(s\right)}{E\left(s\right)}=k_p+k_i{\left(\frac{1}{s}\right)}^{\alpha +j\beta }+k_d{\left(s\right)}^{\gamma +j\delta }$$

(13)

The differential operator in the modified transfer function can be simplified as shown below. ²²

$$\begin{array}{c}\multicolumn{1}{c}{{\left(s\right)}^{\gamma +j\delta }={\left(s\right)}^{\gamma }{\left(s\right)}^{j\delta }}\\ ={\left(s\right)}^{\gamma }{e}^{\ln {\left(s\right)}^{j\delta }}\hfill \\ ={\left(s\right)}^{\gamma }{e}^{j\delta \ln \left(s\right)}\hfill \\ ={\left(s\right)}^{\gamma }[\cos \left(\delta \ln \left(s\right)\right)+j\sin \left(\delta \ln \left(s\right)\right)]\hfill \end{array}$$

(14)

In real-time applications, a feedback control law is required to deliver real output for a real input signal. Hence, the differential operator is approximated by considering only its real part as shown below. ²²

$${\left(s\right)}^{\gamma +j\delta }={\left(s\right)}^{\gamma }\cos (\delta \times \ln \left(s\right))$$

(15)

The integral operator in the modified transfer function can also be simplified as shown below. ²²

$${\left(\frac{1}{s}\right)}^{\alpha +j\beta }={\left(\frac{1}{s}\right)}^{\alpha }{\left(\frac{1}{s}\right)}^{j\beta }$$

(16)

By considering only the real part, the integral operator can be approximated as shown below. ²²

$${\left(\frac{1}{s}\right)}^{\alpha +j\beta }={\left(\frac{1}{s}\right)}^{\alpha }\cos \left(\beta \ln \left(\frac{1}{s}\right)\right)$$

(17)

The simplified transfer function of the control law is thus formulated as shown below.

$$C\left(s\right)=\frac{U\left(s\right)}{E\left(s\right)}=k_p+\frac{k_i}{s^{\alpha }}\cos \left(\beta \ln \left(\frac{1}{s}\right)\right)+k_d{s}^{\gamma }\cos \left(\delta \ln \left(s\right)\right)$$

(18)

The CO-PID transfer function in (18) is quite similar to the FO-PID transfer function shown in (10). The only modification comes in the form of the two new weights assigned to the integral and differential terms. These weights are formulated as the cosines of natural logarithms of integer-order integrals and derivatives of the state error variable. It is to be noted that the coefficient of the complex order’s imaginary part modulates each $\ln (.)$ term in the expression. The CO-PID control law is expressed in equation (19).

$$\begin{array}{c}\multicolumn{1}{c}{u\left(t\right)=k_{p}e\left(t\right)+k_i\left(D^{-\alpha }e\left(t\right)\right)\cos \left(\beta \ln \left|\underset{0}{\overset{t}{\int }}e\left(\tau \right)d\tau \right|\right)}\\ +k_d\left(D^{\gamma }e\left(t\right)\right)\cos \left(\delta \ln \left|\stackrel{·}{e}\left(t\right)\right|\right)\hfill \end{array}$$

(19)

The magnitude of the $\cos (.)$ terms varies between zero and unity. In the transient (or disturbed) state, the magnitude of the error-derivative and error-integral enlarges which reduces the $\cos (.)$ term to zero allowing the response to quickly converge to the set point with minimal control cost. In the steady state, the magnitude of the error-derivative and error-integral enlarges which inflates the $\cos (.)$ term to unity contributing stronger damping against stead-state fluctuations and overshoots. This arrangement significantly improves the response speed and disturbance attenuation capacity of the control law while effectively economizing the overall control energy expenditure and preventing large control input requirements that eventually saturate the actuator. From a functional viewpoint, the $\cos (.)$ term acts as a superior self-regulator that adapts the integral and differential gains of the control law as a nonlinear function of error variables, as shown below.

$$u\left(t\right)=k_{p}e\left(t\right)+k_i\left(t\right)\left(D^{-\alpha }e\left(t\right)\right)+k_d\left(t\right)\left(D^{\gamma }e\left(t\right)\right)$$

(20)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{where},k_i\left(t\right)=k_i\cos \left(\beta \ln \left|\underset{0}{\overset{t}{\int }}e\left(\tau \right)d\tau \right|\right),}\\ k_d\left(t\right)=k_d\cos \left(\delta \ln \left|\stackrel{·}{e}\left(t\right)\right|\right)\hfill \end{array}$$

To maintain a stable control behavior, the variations in the internal constituents of the $\cos (.)$ term are saturated between 0 rad and $\pi /2$ rad. The block diagram of the CO-PID structure is shown in Figure 2. The CO-PID controller parameters are optimized by using the tuning procedure discussed in Section 4.

Figure 2.

Block diagram of the CO-PID control law.

Parameter tuning procedure

Optimizing the controller parameters is very crucial to generate an optimal control yield. However, selecting a unique set of parameters that yield optimum time-domain performance and control efficiency under every operating condition is an ill-posed problem. ³⁰ Owing to the contradictory nature of these factors, generally, a compromise is made between the controller’s time optimality, and control input economy. The analytical tuning methods are hard to track. The empirical settings of the PID gains are constrained by the designer’s experience, and thus, may not always produce precise position regulation or transitory recovery behavior. Hence, in this work, the following quadratic cost function is employed that iteratively minimizes the variations in the control input as well as the classical error. ¹⁵

$$J=\underset{0}{\overset{T}{\int }}(100{\left|e\left(\tau \right)\right|}^2+{\left|u\left(\tau \right)\right|}^2)d\tau$$

(21)

To ensure optimal tuning, a weight of 10 is assigned to error-minimization criteria, so that the transgressions in the BG concentration levels are penalized more than that of the control-minimization criteria. ¹⁵ The controller parameters are optimized via simulations as discussed in Section 5.1. The flow chart in Figure 3 depicts the parameter tuning algorithm. The BMM parameters of patient 1 are used as the reference for the parameter tuning procedure. The offline selection process is initiated by picking a random set of controller parameters (PID gains and fractional/complex orders) from the pre-defined search space. In every simulation trial, the controller parameters are empirically updated and the control system is tasked to regulate the virtual patient’s BG levels at 80 mg/dL from an initial state of hyperglycemia for 500 min. The resulting cost $J_n$ for that trial is then computed; where, $n$ represents the number of trial. The search space is investigated in the direction of the cost function’s steepest gradient descent. ³² If the cost of the present trial ( $J_n$ ) is found to be lesser than the cost of the previous trial ( $J_{n-1}$ ), the local minimum-cost variable $J_{\min }$ is updated. This arrangement ensures that the exploration proceeds in the direction of the descending gradient of $J$ . The search for the best-fit parameter values is terminated if either $J_{\min }$ achieves a predefined threshold cost or the algorithm has completed the maximum number of trials ( $n_{\max }$ ) allowed. In this work, the threshold for $J_{e,\min }$ is preset at 1 × 10⁷ based on the expert’s experience. Owing to the long time-duration of each simulation, the value of $n_{\max }$ is set at 12. Finally, the acquired set of parameters values are manually fine-tuned to further refine the selections. The selected parameter values are shown in Table 2. The same procedure is employed to tune all the parameters offline.

Figure 3.

Flow chart of the parameter tuning algorithm.

Table 2.

Optimized controller parameters.

Parameters	Selection range	Initial value	Optimized values
			FO-PID	CO-PID
$k_p$	[0, 0.1]	1 × 10−3	3.55 × 10−3	3.71 × 10−3
$k_i$	[0, 0.1]	1 × 10−4	1.19 × 10−4	1.24 × 10−4
$k_d$	[0, 0.1]	1 × 10−3	0.0478	0.0511
$\alpha$	[0, 1]	0.1	0.675	0.538
$\gamma$	[0, 1]	0.1	0.395	0.322
$\beta$	[0, 1]	0.01	0	0.0519
$\delta$	[0, 1]	0.01	0	0.0762

Simulations and analysis

This section presents the simulation-based testing of the designed BG regulation control schemes along with a comprehensive analysis of the acquired results.

Simulation setup

The time-domain performance of the FO-PID and CO-PID control schemes is comparatively analyzed via credible simulations. The customized control application is implemented and the simulations are carried out using the MATLAB/Simulink R2018b software environment. ³³ The software is operated on a personal computer that is equipped with a 64 -bit, 2.4 GHz CPU, and 8.0 GB RAM. The fractional operators in the FO-PID and CO-PID control laws are implemented by using the built-in functions in MATLAB’s FOMCON toolbox. The sampling time is set at 1.0 min. As discussed in Section 2, the dynamics of the virtual diabetic patients are simulated using three different BMMs to validate the BG regulation performance of the control laws to ascertain the controller’s performance under parametric variations. The BMM parameters for three different virtual patients are also considered in the presence of sensor noise and meal disturbance. The three different sets of model parameters are already listed in Table 1. The IIR signal $u\left(t\right)$ is restricted within the limits of 0 mU/min and 100 mU/min via a saturation function of the form, $100\mathrm{sat}\left(u\left(t\right)\right)$ , to prevent the BG levels from plunging significantly, which may put the patient in a hypoglycemic state. ¹⁰ The control saturation function is described as follows.

$$100\mathrm{sat}\left(u\left(t\right)\right)=\{\begin{array}{c}100,u\left(t\right)>100\hfill \\ \multicolumn{1}{c}{u\left(t\right),0\le u\left(t\right)\le 100}\\ \multicolumn{1}{c}{0,u\left(t\right)<100}\end{array}$$

(22)

Simulation results

The efficacy of the proposed CO-PID controller is analyzed against the FO-PID controllers via the following two test cases. The tests are conducted on the three virtual patients whose BMM parameters are identified in Table 1. In each test case, the closed-loop system is tasked to track the BG level set-point of 80 mg/dL. White Gaussian noise is introduced in $G\left(t\right)$ in each test case to emulate the effects of sensor noise on the reference tracking behavior.

• A. BG regulation under normal conditions: This simulation serves to examine the capability of the designed control law to regulate the BG levels of the three patients to 80 mg/dL from an initial state of hyperglycemia (200 mg/dL). To emulate the effects of measurement noise contributed by the glucose sensor, a white Gaussian noise signal having a mean of zero and variance of 0.2 is introduced in $G\left(t\right)$ at the beginning of the trial. The time-domain profiles of BG level and IIR (control input) under FO-PID and CO-PID controllers for the Patients 1, 2, and 3 are depicted in Figures 4 to 6, respectively.

• B. BG regulation under meal disturbance: This simulation serves to analyze the controller’s resilience against the bounded exogenous disturbances that are normally caused by the meal intake. The effects of meal disturbance and sensor noise are emulated by administering a simulated impulse signal in $G_m\left(t\right)$ and a white Gaussian noise signal in $G\left(t\right)$ , respectively. The closed-loop system is tasked to recover from the transient disturbance and normalize the BG levels of each patient to 80 mg/dL from an initial state of hyperglycemia (200 mg/dL). The simulated impulse signal of amplitude 80 mg/dL is injected in the system at t ≈ 500 min, whereas, the white noise signal having a mean of zero and variance of 0.2 is introduced in the system at the beginning of the trial. The time-domain profiles of BG level and IIR under FO-PID and CO-PID controllers for the Patients 1, 2, and 3 are depicted in Figures 7 to 9, respectively.

Figure 4.

(a) BG levels of Patient 1 under normal conditions and (b) IIR (control input) for Patient 1 under normal conditions.

Figure 5.

(a) BG levels of Patient 2 under normal conditions and (b) IIR (control input) for Patient 2 under normal conditions.

Figure 6.

(a) BG levels of Patient 3 under normal conditions and (b) IIR (control input) for Patient 3 under normal conditions.

Figure 7.

(a) BG levels of Patient 1 under meal disturbance and (b) IIR (control input) for Patient 1 under meal disturbance.

Figure 8.

(a) BG levels of Patient 2 under meal disturbance and (b) IIR (control input) for Patient 2 under meal disturbance.

Figure 9.

(a) BG levels of Patient 3 under meal disturbance and (b) IIR (control input) for Patient 3 under meal disturbance.

Analytical discussion

The simulation results are examined as per the following Key-Performance-Metrics (KPMs).

RMSE	: Root-mean-squared value of error $e\left(t\right)$ in the BG levels.
IAE	: Integral of the absolute value of error $e\left(t\right)$ in the BG levels.
Tfall	: Time taken by the BG level to fall to +10% of $G_{\mathrm{ref}}$ from an initial state of hyperglycemia.
Tset	: Time taken by the BG level to normalize and settle within ±2% of $G_{\mathrm{ref}}$ .
Trec	: Time taken by the BG level to recover from a transient disturbance.
IIRMS	: Mean-squared value of the IIR.
IIRp,Start	: Peak value of the IIR during initial start-up
IIRp,Dist	: Peak value of the IIR during transient disturbances.

The RMSE is calculated as shown below.

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}\sum _{n=0}^N{\left|e\left(n\right)\right|}^2}$$

(23)

where, $N$ is the total number of samples and $n$ is number of measurement. The value of $N$ is 500 for simulation A and 1000 for simulation B in this article. The IAE is calculated as shown below.

$$\mathrm{IAE}=\int \left|e\left(\tau \right)\right|d\tau$$

(24)

The IIR_MS is calculated as shown below.

$$\mathrm{II}{\mathrm{R}}_{\mathrm{MS}}=\frac{1}{N}\sum _{n=0}^N{\left|u\left(n\right)\right|}^2$$

(25)

These KPMs are used to quantitatively analyze the performance and robustness of the designed BG regulation control procedures. The simulation results are summarized in Table 3. The quantitative analysis of the simulation results validates the superior time optimality and disturbance-response speed of the CO-PID controller as compared to the FO-PID controller.

Table 3.

Summary of simulation results.

Simulation	KPM		Patient 1		Patient 2		Patient 3
	Metric	Unit	FO-PID	CO-PID	FO-PID	CO-PID	FO-PID	CO-PID
A	RMSE	mg/dL	27.7	26.7	35.3	29.7	30.5	24.9
	IAE	mg/dL	1586.3	788.7	1719.1	811.3	1433.0	855.4
	Tfall	min.	136	128	191	148	186	104
	Tset	min.	228	201	302	230	328	162
	IIRMS	(mU/min)2	7.01	2.89	6.14	2.92	5.20	2.12
	IIRp,Start	mU/min	7.72	5.27	7.67	5.34	7.80	4.98
B	RMSE	mg/dL	30.59	24.38	32.60	26.25	30.78	24.62
	IAE	mg/dL	1962.0	1252.9	2028.1	1277.3	1961.4	1591.8
	Trec	min.	218	154	352	180	229	198
	IIRMS	(mU/min)2	4.04	2.25	4.78	2.30	4.48	2.48
	IIRp,Dist	mU/min	4.83	2.40	4.75	2.83	4.53	2.68

In Simulation A, each controller exhibits a different transient response to normalize the BG levels of each patient at the desired level. The FO-PID controller demonstrates a mediocre reference tracking accuracy and a relatively higher settling time to converge to $G_{\mathrm{ref}}$ . As compared to FO-PID, the CO-PID controller exhibits a considerable improvement in the reference-tracking behavior by improving the response speed and the regulation accuracy. The CO-PID displays the fastest transient response and an accurate reference tracking behavior. It quickly settles the BG levels at the desired set-point and effectively suppresses the steady-state fluctuations while maintaining an economical application of the control input (IIR).

In Simulation B, each controller applies the necessary control resources to normalize the BG levels from a state of hyperglycemia caused by the simulated meal disturbance. The FO-PID controller exhibits weak immunity against meal disturbance and sensor noise. The response shows relatively higher chattering and slower transient recovery time. It also yields a highly discontinuous control activity with large insulin infusion requirements. Furthermore, the FO-PID controller injects the insulin for a longer time duration, which is therapeutically unsuitable for the patients from a practical point of view. The CO-PID controller shows the most time optimal behavior with a significant improvement in robustness against meal disturbance and sensor noise. It exhibits the relatively faster referencing response speed, minimum-time transient recovery, and minimal tracking error. It also yields a smooth and economical control yield, which in turn minimizes the chattering content in the state response.

As compared to the FO-PID controller, the CO-PID controller renders a mean improvement of 13.1%, 56.8%, 25.9%, 33.4%, 30.9%, 28.1%, and 26.2% in the RMSE, IAE, T_fall, T_rec, T_set, IIR_MS, and IIR_p, respectively. The significant improvement in the transient response speed and the control input efficiency of the CO-PID controller is credited to the enhancement in the controller’s degree of freedom contributed by the allocation of complex orders to integral and differential operators. This arrangement indirectly modifies the derivative and integral gains online to improve the controller’s adaptability to execute better reference tracking accuracy and transient recovery response, even under exogenous disturbances. Furthermore, it achieves the aforementioned feats without compromising the control input (IIR) economy which is yet another remarkable milestone. Although the CO-PID control procedure requires the offline tuning of seven distinct parameters which is indeed a laborious task; however, the enhanced time optimality offered by the said procedure exceeds this disadvantage.

The proposed CO-PID controller’s performance is also compared with the state-of-the-art controllers presented in Delavari et al. ¹⁰ to further validate its effectiveness in addressing the BG regulation problem. These state-of-the-art control schemes include the super-twisting sliding-mode-controller (STSMC), higher-order SMC (HOSMC), fractional-order SMC (FOSMC), and adaptive fractional-order SMC (AFOSMC). The aforementioned controllers were applied to the same set of three patient models, as described in Table 1, to carry out the simulations A and B in Delavari et al. ¹⁰ The reader is referred to Delavari et al. ¹⁰ for visualizing the simulation results of the aforementioned controllers. In this article, the performance comparison is done on the basis of the IAE metric expressed in Guras et al. ²⁴ . The summary of performance comparison is quantified in Table 4. The quantitative analysis also validates the enhanced BG regulation capability of the CO-PID controller as compared to other state-of-the-art controllers. The AFOSMC is declared as the proposed control law in Delavari et al. ¹⁰ As compared to the AFOSMC, the CO-PID controller demonstrates an improvement of 10.0%, −6.8%, and 15.7% in the IAE of Patients 1, 2, and 3, respectively, for Simulation A. For Simulation B, the CO-PID controller exhibits a relative improvement of 11.6%, −5.9%, and 3.8% in the IAE of Patients 1, 2, and 3, respectively, as compared to AFOSMC.

Table 4.

Performance comparison with other controllers.

Simulation	Metric	Unit	Controllers	Patient 1	Patient 2	Patient 3
A	IAE	mg/dL	STSMC 10	1956	2110	1880
			HOSMC 10	1456	1357	1609
			FOSMC 10	1135	916	1542
			AFOSMC 10	877	759	1015
			CO-PID (proposed)	789	811	855
B	IAE	mg/dL	STSMC 10	2720	2986	3560
			HOSMC 10	2014	1895	2398
			FOSMC 10	1766	1660	2228
			AFOSMC 10	1417	1206	1655
			CO-PID (proposed)	1253	1277	1592

A special comparison case

To better ascertain their effectiveness and time-optimality of the FO-PID and CO-PID controller in regulating the BG levels, the performance of these controllers with various orders of the differential operator is being analyzed in this section. The analysis is done by carrying out additional simulations of test cases A and B.

Firstly, two new FO-PID controller variants are introduced and their behavior is benchmarked against the original FO-PID controller (See Table 2). The original control law is referred to as FO-PID-1 controller for this comparison. The value of $\gamma$ in the FO-PID-2 controller variant is incremented by 10% while keeping all other controller parameters affixed at their original values. Conversely, the value of $\gamma$ in the CO-PID-3 controller variant is decremented by 10% while keeping all other parameters affixed at their original values. The changes in the orders of the fractional derivative terms are highlighted in Table 5. The simulations A and B are carried out for Patient 2 under the influence of the aforementioned FO-PID controller variants. The simulation results are shown in Figures 10 and 11, respectively. The simulation results are quantified in Table 6. The results show that increasing the magnitude of the differential operator’s fractional order minimizes the settling time and transient recovery time at the cost of generating a highly discontinuous control activity which introduces chattering in the BG levels, which is not desirable. On the contrary, reducing the said fractional order slows down the BG regulation response speed while maintaining a relatively smoother control activity. In terms of BG regulation and IIR input behavior, the FO-PID-1 controller outperforms the FO-PID-2 and FO-PID-3 controller.

Table 5.

Modified FO-PID controller variants.

Parameters	Modified derivative’s fractional order
	FO-PID-1	FO-PID-2	FO-PID-3
$\gamma$	0.395	0.435	0.355

Figure 10.

(a) BG levels of Patient 2 with different FO-PID controller variants under normal conditions, (b) IIR (control input) for Patient 2 with different FO-PID controller variants under normal conditions.

Figure 11.

(a) BG levels of Patient 2 with different FO-PID controller variants under meal disturbance, (b) IIR (control input) for Patient 2 with different FO-PID controller variants under meal disturbance.

Table 6.

Simulation results with different FO-PID controller variants.

Simulation	KPM		FO-PID-1	FO-PID-2	FO-PID-3
	Metric	Unit
A	RMSE	mg/dL	35.3	32.4	38.2
	IAE	mg/dL	1719.1	1588.5	1862.2
	Tfall	min.	191	185	208
	Tset	min.	302	295	318
	IIRMS	(mU/min)2	2.65	3.05	2.44
	IIRp,Start	mU/min	7.67	9.11	6.32
B	RMSE	mg/dL	32.60	30.1	35.8
	IAE	mg/dL	2028.1	1902.4	2182.9
	Trec	min.	352	295	365
	IIRMS	(mU/min)2	4.78	5.15	4.57
	IIRp,Dist	mU/min	4.75	5.68	4.09

Similar to the process discussed above, two new CO-PID controller variants are also introduced and their behavior is benchmarked against the original CO-PID (CO-PID-1) controller (optimized in Table 2). The values of both parameters $\gamma$ and $\delta$ in the CO-PID-2 controller variant are incremented by 10% while keeping all other parameters fixed at their original values. Conversely, the values of both parameters $\gamma$ and $\delta$ in the CO-PID-3 controller variant are decremented by 10% while keeping all other parameters fixed at their original values. The changes introduced in the orders of the fractional derivative terms are highlighted in Table 7. The simulations A and B are carried out for Patient 2 under the influence of the aforementioned CO-PID controller variants. The simulation results are shown in Figures 12 and 13, respectively. The simulation results are quantified in Table 8. The results show that increasing the magnitude of the complex order of the differential operator improves the transitional times of the BG regulation response under normal as well as disturbance conditions while generating a disrupted control activity and large control energy requirements. The chattering in the response is also inevitable in this case. Conversely, reducing the said complex order makes the BG regulation response quite sluggish speed; however, it also economizes the control energy expenditure. The overall performance of the CO-PID-1 controller surpasses the CO-PID-2 and CO-PID-3 controller as per the BG regulation and IIR input.

Table 7.

Modified CO-PID controller variants.

Parameters	Modified derivative’s fractional order
	CO-PID-1	CO-PID-2	CO-PID-3
$\gamma$	0.322	0.354	0.290
$\delta$	0.0762	0.0838	0.0686

Figure 12.

(a) BG levels of Patient 2 with different CO-PID controller variants under normal conditions, (b) IIR (control input) for Patient 2 with different CO-PID controller variants under normal conditions.

Figure 13.

(a) BG levels of Patient 2 with different CO-PID controller variants under meal disturbance, (b) IIR (control input) for Patient 2 with different CO-PID controller variants under meal disturbance.

Table 8.

Simulation results with different CO-PID controller variants.

Simulation	KPM		CO-PID-1	CO-PID-2	CO-PID-3
	Metric	Unit
A	RMSE	mg/dL	29.7	28.5	31.6
	IAE	mg/dL	811.3	788.2	896.5
	Tfall	min.	148	129	177
	Tset	min.	230	221	268
	IIRMS	(mU/min)2	2.17	3.84	1.98
	IIRp,Start	mU/min	5.34	6.22	4.34
B	RMSE	mg/dL	26.3	25.0	28.1
	IAE	mg/dL	1277.3	1068.4	1398.6
	Trec	min.	180	168	198
	IIRMS	(mU/min)2	2.30	4.18	2.05
	IIRp,Dist	mU/min	2.83	2.06	4.18

Conclusion

This article methodically formulates and verifies the efficacy of a flexible CO-PID controller to improve the BG normalization behavior in TID patients, especially under the influence of meal disturbances and exogenous noise sources. The BG regulation in TID patients via closed-loop feedback control procedure is a difficult task. The proposed control procedure is realized by augmenting the integral and differential operators of the ubiquitous PID baseline controller with complex order. The controller’s transformation into complex order controllers increases its design flexibility and time optimality, thus equipping it to exhibit a faster tracking speed and stronger immunity against disturbances.

The simulation results validate the superior robustness and time optimality of the CO-PID controller. Additional simulations are also conducted to further analyze the time-optimality of the FO-PID and CO-PID controllers under the influence of different orders of differential operators. These simulations also validate the optimum BG regulation behavior of the originally prescribed CO-PID control law. In the future, the CO-PID controller can be augmented with intelligent systems, online supervised learning laws, or model-free expert adaptive systems to dynamically self-tune the controller parameters to further robustify its performance against hyperglycemia and exogenous disturbances. Meta-heuristic optimization algorithms can be investigated to improve the controller’s parameter tuning process. Furthermore, the type-3 fuzzy systems and controllers should also be investigated in the future to analyze their ability to further improve the performance of the closed-loop BG regulators.