Metaheuristic Patient Estimation Based Patient-Specific Fuzzy

Article pubs.acs.org/IECR

Metaheuristic Patient Estimation Based Patient-Specific Fuzzy Aggregated Artificial Pancreas Design V. Kirubakaran,† T. K. Radhakrishnan,*,† and N. Sivakumaran‡ †

Department of Chemical Engineering and ‡Department of Instrumentation and Control Engineering, National Institute of Technology, Tiruchirappalli 620015, India S Supporting Information *

ABSTRACT: Patient-specific artificial pancreas design has been receiving increasing attention lately. In this article, using the chaotic bat algorithm (CBA), Hovorka−Wilinska (H−W) model parameters are estimated from nominal H−W virtual patient data. Using this identified H−W model for the virtual patient, multiple empirical second-order plus delay time (SOPDT) models representing glucose−insulin dynamics are derived for the range of blood glucose concentrations (BGCs) considered. Clustering of these models using the k-means algorithm yields three distinct clusters. Implicitly enumerated multiparametric model predictive controllers (mpMPCs) are designed using the cluster representatives. A fuzzy logic aggregation (FLA) of prediction and control improves the design parsimony. An insulin on board (IOB) safety trigger is designed using FLA of multiple full-order linearized CBA-estimated H−W models. The FLA-based mpMPC along with IOB and meal estimation are implemented on an embedded platform and by hardware-in-the-loop (HIL) simulation. In silico trials of the regulation of multiple meal disturbances are performed on the nominal H−W patient in MATLAB through serial communication. With meal estimation errors and varying insulin sensitivity, a very good low blood glucose index (LBGI) of 0, R ≥ 0

x(k + N ) = AN x(k) +

0⎤ ⎥ 0⎥ 0⎥ ⎥ 0⎥ 1⎥ −1⎥⎦

0 R 0 0

The control problem in eq 17 is transformed to the form provided by

Δu i(n)T R Δu i(n)

u ilb ≤ u i(n) ≤ u iub

··· ··· 0 ⋱ ··· ···

z = u i + H −1gx(k)

x T(n) Qx(n)

n=k+1 k + Nc − 1

+

··· ··· ⋱ 0 ··· ···

⎡R ⎢ 0 Γu = ⎢ ⎢0 ⎢⎣ 0

M relates the rate of insulin infused (Δui) to the rate constraints (Δuilb and Δuiub) using L. O represents the state constraint; ws. With no state constraints considered in this work, its value is maintained as in eq 17. Consider a new variable z that is a piecewise-affine function of state vector x in its vector space defined by

min J[x(k)] = x T(k + Np) Px(k + Np) ui

(17)

where

x(k + 1) = Ax(k) + Bu i(k) y(k + 1) = Cx(k + 1)

Mu i ≤ L + Ox(k)

(21)

The active constraints for the given state vector x are determined by the Lagrange multipliers obtained from eq 21.

(16) 15059

dx.doi.org/10.1021/ie5009647 | Ind. Eng. Chem. Res. 2014, 53, 15052−15070

Industrial & Engineering Chemistry Research

Article

Table 5. Fuzzy Weights for Aggregation model no.

weight 1

1 5 11 15 17 19 23 24 26

0.9900 0.9994 0.5480 0.0004 0 0.1361 0 0 0.0002

1 4 7 12 14 16 18 20 25

0.9930 0.9997 0.9507 0.0067 0.0006 0.0042 0.0084 0 0.0002

weight 2 Patient 1 0.0091 0.0006 0.4501 0.9996 1.0000 0.5740 0 0 0.0002 Patient 2 0.0066 0.0003 0.0490 0.9932 0.9994 0.9955 0.9908 0 0.0002

weight 3

cluster

0.0009 0 0 0 0 0.2899 1.0000 1.0000 0.9996

1 1 1 2 2 2 3 3 3

0.0004 0 0.0002 0.0001 0 0.0003 0.0008 0.9999 0.9996

1 1 1 2 2 2 2 3 3

Figure 6. Insulin kinetics of the H−W model.

For the three cluster representatives, the controller is designed as detailed above for both patients. Each row in matrices M, L, and X (defined in eq 17) represents a constraint. The active set of constraints for a given state vector is determined based on the values of the Lagrange multipliers. A positive value of λ indicates the presence of a constraint. The condition for a state vector to belong to an active constraint set is provided by D*(z , x) ≜ {c ∈ I |Mcz − Lc − Xcx ≤ 0}

5. PHYSICAL CONSTRAINTS The constraints for artificial pancreas design predominantly include insulin limits and also the rate of change of the insulin infusion rate, as shown in eq 15. Here, the lower bound is also inclusive of the basal insulin flow rate. A safety constraint based on the active insulin present in the body is described.16 The calculation is based on the insulin infused against the decay rate of the same, as shown in the equation

(22)

where Mc, Lc, and Xc are obtained by concatenating the rows of matrices M, L, and X whose corresponding constraints are active. The set of inactive constraints F is described by the equation F(z , x) ≜ I D*(z , x)

Nd

u iob(k) =

∑ u i(k − d)hd d=1

(23)

(26)

51

where hd is the decay vector. Based on the IOB, using the correction factor (cF, the variation in BGC for 1 U of insulin infused) of the individual patient, it is possible to calculate the amount of BGC that can be regulated at a given instance using the insulin already infused. It has been reported that this uiob value can be included as an explicit constraint in controller formation by appending it to the reducedorder system state vector space.16 In such a case, the dynamics of uiob is constrained by a generic decay rate of insulin infused. Inclusion of this patient-specific approach in IOB calculations (as an explicit constraint) requires that the entire set of internal states of the H−W model be appended to the reduced-order system state space used in controller design. Instead, in this study, the constraint is implemented as an external safety instrument during closed-loop BGC regulation, preserving customization. Based on the cF and accumulation of infused insulin, the amount of BGC that can be compensated is calculated. Using the deviation of the current BGC measurement, a decision to suspend any further infusion of insulin is made according to eq 31 (below). To customize this safety constraint for each patient, in this study, an FLA-based multiple-model composite design is proposed. The models are based on the Jacobian of the CBA model given by

For the active set of constraints observed, a linear independence constraint qualification (LICQ) test is performed. This involves evaluation of the linear independence of constraints using the rank of the active set D*. A full-rank active set D* indicates qualification of LICQ for the given set of constraints. A strict complementary slackness (SCS) test ensures that no λ value in the inactive set is equal to zero (which indicates the presence of weakly active constraints), the presence of which will make the feasible quadratic programming problem degenerate. Upon satisfying these two conditions, the constraint set is considered to be an optimal set. Based on the set of inactive constraints for a given state vector, by concatenating the rows of these matrices (MF, LF, and XF) based on inactive constraints, the polyhedral active region for a constraint set is obtained using the equation CR D ≜ Δ{x ∈ Θ: MFz(x) − L F − XFx ≤ 0, λM(x) > 0} (24)

The multiparametric control action to suppress active constraints is provided by u i(k) = H −1Mc T(McH −1Mc T)−1Xcx(k) − H −1gx(k) + H −1Mc T(McH −1Mc T)−1Lc

(25) 15060



AHWn

Article

⎡− x1ss k12 0 0 0 ⎢ ⎢ x1ss − (k12 + x 2ss) 0 0 0 ⎢ ⎢ ⎛ ⎞ Vmax,LD ⎢ 0 ⎟⎟ − ⎜⎜ka1 + 0 0 0 ⎢ kM,LD + S1ass ⎠ ⎝ ⎢ ⎢ ⎛ Vmax,LD ⎞ ⎢ 0 ⎟⎟ − ⎜⎜ka2 + 0 0 0 ⎢ kM,LD + S1bss ⎠ ⎝ ⎢ − ka1 * ka1 * 0 0 =⎢ 0 ⎢ ka2 * ka1 * 0 0 ⎢ 0 ⎢ ⎢ 0 0 0 0 0 ⎢ ⎢ ⎢ 0 0 0 0 ⎢ 0 ⎢ ⎢ ⎢ 0 0 0 0 0 ⎢⎣

0

− Q 1ss

0

0

Q 1ss

− Q 2ss

0

0

0

0

0

0

0

0

0

− ke

0

0

k b1 Vi

− ka1

0

k b2 Vi

0

− ka2

k b3 Vi

0

0

− EGP0 ⎤ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ − ka3 ⎥ ⎥⎦

BHWn = [0 0 k ins (1 − k ins) 0 0 0 0 0]T ⎡ 18.018 ⎤ 0 0 0 0 0 0 0 0⎥ C HWn = ⎢ V g ⎣ ⎦ n = 1, 2, ..., Nn x HW = [Q 1 Q 2 S1a S1b S2 S3 x1 x 2 x3 ]T

(27)

The steady-state values of state variables for each region’s cluster representative are obtained in MATLAB by using fsolve application peripheral interphase (API). A fuzzy aggregation of the linear models based on the current region of operation is performed as follows Nn

AHW =

6. EMBEDDED IMPLEMENTATION A practical approach to developing an AP as a product has been reported. Issues related to CGM sensor calibration, lost data during communication, and patient uncertainties pose great challenges in such a hardware implementation. The patient is the mathematical nominal H−W model implemented in MATLAB. The controller designed in previous sections is implemented on a hardware platform. HIL trials are necessary to obtain approvals for clinical trials.9 The data flow diagram of the HIL scheme is shown in Figure 7. To implement the control strategy, a 32-bit microcontroller (PIC32mx795H512L) that operates using an 8 MHz external clock was chosen. The embedded development tool chain includes a C 32 compiler, MPLAB IDE V 8.84, and MPLAB REAL ICE in-circuit emulator (ICE), all from MICROCHIP. Because the control strategy is a simple state feedback gain-based approach and also because the FLA and IOB are based on simple arithmetic operations, it is not mandatory to implement this strategy on a 32bit platform. As shown in Figure 7, the modular design ensures easier migration of software into other families of microcontrollers and related tool chains with minor changes in the hardware application layer (HAL).

Nn

∑ wnmAHWn ,

BHW =

n=1

∑ wnmBHWn , n=1

Nn

C HW =

∑ wnmC HWn n=1

(28)

Absorption of insulin infused in the H−W model is classified into two categories: (i) local degradation (LD), which includes a slow channel (LDa) and a fast channel (LDb), and (ii) plasma insulin absorption. This is shown in Figure 6. Under steady state, insulin absorption can be defined as u iss = LDass + LDbss + (ka1sS2ss + ka2S1bss)

(29)

The accumulation of insulin (uacc) can hence be calculated as uacc = uacc + {u i(k) + ub(k) − [LDa (k) + LDb(k) + ka1S2(k) + ka2S1b(k)]}Ts

(30)

7. RESULTS AND DISCUSSION The controller design is deployed on the embedded microcontroller, whereas the patient is deployed in MATLAB. Serial communication using a master (embedded board)−slave (MATLAB) configuration transfers insulin and BGC data periodically. As described in section 6, CGM sensor error, data loss in the CGM sensor (simulated using MATLAB), and patient uncertainties are considered during simulation and HIL experiments. The CGM sensor error [θ(k)], given by ⎡ e(k) − γ ⎤ θ(k) = ε + ζ sinh⎢ ⎥ ⎣ ⎦ δ (32)

where Ts is the sampling time in minutes and ub is the bolus insulin. Based on the switching function defined by ⎧ if uacccF < BGC − reference ⎪ u i (k ) u i (k ) = ⎨ ⎪ ⎩ ubas if uacccF > BGC − reference

(31)

insulin infusion is manipulated if IOB poses a risk of hypoglycemia. Based on the given condition, the insulin infusion is reduced to the basal flow rate (ubas) during possible hypoglycemic risk. 15061



Article

Figure 7. Data flow diagram of hardware-in-the-loop testing of an artificial pancreas.

where ε, δ, γ, and ζ represent the Johnson transformation parameters,26 is designed to be a non-Gaussian, nonwhite error with a non-normal distribution. The ARMA model e(k) = 0.7[e(k − 1) + φ(k)]

mean square error (MSE) of control for each scenario is included. The estimated average glucose (eAG) can be used to find A1c according to

(33)

eAG = 28.7A1c − 46.7

constitutes a white Gaussian noise [φ(k)] assumed to have unity covariance and zero mean. Loss of data during closed-loop control poses a significant threat, with the possibility of hypoglycemic episodes if undetected. Data losses for time periods of 10 min30 and up to 20 min31 were considered in previous controller designs. In this study, a data-loss period for a longer duration of 30 min is considered. During data loss, the insulin infusion rate reverts back to the basal rate, as indicated by

or A1c =

eAG + 46.7 28.7

(36)

LBGI and HBGI evaluate the frequency and extent of risks by means of BGC readings. f(BGC) provides risk functions for low and high BGC levels as follows f (BGC) = 1.509[ln(BGC)1.084 − 5.381]

(37)

The quadratic risk function is defined by

u i(k) = ubas ub(k) = 0

(35)

s(BGC) = 10f (BGC)2

(34)

(38)

The individual quadratic risk functions for both low and high BGC are thus given by

To maintain uniform test conditions, the CGM sensor noise and data-loss period remain unchanged between different scenarios. To quantify the performance of controllers for specific applications, it is imperative that measures beyond error-related indices be included. Apart from the visualization of BGC regulation, safety indices used by physicians to quantify controller performance are employed here. These measures include mean BGC, maximum BGC, minimum BGC, A1c, low blood glucose index (LBGI), and high blood glucose index (HBGI).1 Also, the

⎧ if f (BGC) < 0 ⎪ s(BGC) s l(BGC) = ⎨ ⎪ if f (BGC) ≥ 0 ⎩ 0 ⎧ if f (BGC) > 0 ⎪ s(BGC) s h(BGC) = ⎨ ⎪ if f (BGC) ≤ 0 ⎩ 0 15062

(39)



Article

Table 6. Testing Scenario Selectiona meal-estimation error

a

insulin sensitivity

scenario

mpMPC

GS

FLA

FF

bolus

IOB

+25%

−25%

+20%

−20%

varying

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

× × × × × × × × × × × × × × × ×

− × − − − − − − − − − − − − − −

− − × × × × × × × × × × × × × ×

− − − × × − × × × × × × × × × ×

− − − − × × − × × × × × × × × ×

− − − − − − × × × × × × × × × ×

− − − − − − − − × × × × − − − −

− − − − − − − − − − − − × × × ×

− − − − − − − − × − − − × − − −

− − − − − − − − − − × − − − × −

− − − − − − − − − − − × − − − ×

×, engaged; −, disengaged.

Table 7. Quantitative Controller Performance: Patient 1 index scenario

maximum BGC (mg/dL)

minimum BGC (mg/dL)

eAG (mg/dL)

A1C (%)

LBGI

HBGI

MSE

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

240.317 228.357 237.280 222.107 203.984 215.478 231.511 220.725 192.262 208.061 233.027 223.595 185.919 206.720 232.920 207.487

60.369 64.195 63.032 58.881 63.358 65.694 74.690 64.825 75.983 88.315 92.729 73.296 74.444 88.903 88.191 74.757

119.957 118.293 119.783 112.096 107.985 114.899 132.269 129.265 120.607 132.084 145.670 132.897 116.523 131.136 144.618 129.432

5.807 5.749 5.801 5.533 5.390 5.631 6.236 6.131 5.830 6.229 6.703 6.258 5.687 6.196 6.666 6.137

1.694 1.546 1.202 1.816 1.691 1.020 0.357 0.281 0.528 0.107 0.041 0.263 0.512 0.099 0.044 0.166

2.297 1.933 1.834 1.318 0.864 1.161 2.749 2.052 1.260 2.004 3.512 2.181 0.903 1.921 3.372 1.707

2042.333 1699.986 1576.828 1265.967 896.066 1013.749 2130.928 1567.276 997.750 1477.737 2633.158 1711.036 733.377 1419.432 2533.610 1317.425

comparisons of the best strategies from the first stage inclusive of the IOB and bolus modules. Each control strategy involves the mpMPC alongside other modules. From Figure 8, it is evident that, instead of a single linear mpMPC (scenario 1) for the entire BGC range, a zonebased gain scheduling approach (scenario 2, mpMPC + GS) marginally improves the performance of the controller during multiple meal disturbances by reducing the maximum BGC and also provides a lower HBGI (see Tables 7 and 8) . An FLAbased strategy (scenario 3, mpMPC + FLA) improves the LGBI while further reducing the HBGI. It is observed that, as the model selection for control deployment becomes more region-specific, a further improvement in closed-loop control results. Meal announcement (scenario 4, mpMPC + FLA + FF), which is added in the future prediction of the BGC, increases the aggressiveness of the control. This increases the risk of hypoglycemia, reflected in the high LGBI value. Although this combination has the lowest MSE, it can be seen that the risk factors are increased as a result of the increase in LBGI and the low BGC metrics.

From the BGC information, using eqs 37−39, the LBGI and HBGI are calculated as LBGI =

HBGI =

Ne

1 Ne

∑ sl(BGCs)

1 Ne

∑ sh(BGCs)

s=1 Ne s=1

(40)

Different scenarios to investigate the performance of the controller were executed. The various modules involved in each scenario are pisted in Table 6. The values of these metrics obtained for patients 1 and 2 are available in Tables 7 and 8, respectively. The modules consist of a linear mpMPC, gain scheduling (GS), FLA, meal announcement as a feed-forward (FF) disturbance, bolus insulin (BOL) delivered 30 min prior to meal time, and an IOB-based safety trigger. The combination of these modules is separated into two stages. The first stage consists of comparisons of combinations that do not include the bolus and IOB combination. The second stage consists of 15063



Article

Table 8. Quantitative Controller Performance: Patient 2 index scenario

maximum BGC (mg/dL)

minimum BGC (mg/dL)

eAG (mg/dL)

A1C (%)

LBGI

HBGI

MSE

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

248.910 235.215 234.452 247.038 213.622 207.383 253.225 214.341 202.333 216.206 220.436 214.975 194.536 219.126 239.400 239.755

64.970 65.782 66.191 67.043 71.884 71.294 83.189 80.473 59.678 82.602 86.981 75.047 61.617 80.943 104.202 64.826

125.721 123.990 124.097 131.133 126.435 119.049 141.581 133.800 121.878 135.249 137.587 125.379 119.312 134.559 149.788 132.697

6.008 5.947 5.951 6.196 6.033 5.775 6.560 6.289 5.874 6.340 6.421 5.996 5.784 6.316 6.846 6.251

1.396 1.264 0.984 0.809 0.588 0.772 0.095 0.013 0.406 0.015 0.008 0.171 0.439 0.019 0.002 0.321

2.757 2.354 2.130 2.874 1.938 1.326 3.565 2.028 1.379 2.175 2.380 1.510 1.154 2.139 3.769 2.884

2220.436 1855.150 1634.714 2147.636 1371.425 1000.932 2513.892 1291.023 945.913 1391.914 1533.738 1000.144 804.653 1375.939 2500.588 2071.016

Figure 8. Stage 1 comparison of control combinations.

scenario 4 (scenario 5, mpMPC + FLA + FF + BOL) and scenario 3 (scenario 6, mpMPC + FLA + BOL). IOB provides a means to prevent excessive insulin infusion during the regulation of meal-ingested BGC release. The impact of IOB on scenario 4 (scenario 7, mpMPC + FLA + FF + IOB) in the absence of bolus insulin was studied and compared with its combination with scenario 5 (scenario 8, mpMPC + FLA + FF + BOL + IOB). The responses of scenarios 5−8 are presented in

The impact of bolus insulin on controlling high BGC exposure was studied. In general, a preemptive action by bolus insulin on the influence of meal ingestion on the BGC lowers the need for a reactive insulin infusion by the feedback controller. This improves the LBGI measure and also improves the HBGI measure because the maximum levels of BGC are reduced significantly. The impact of bolus insulin on scenarios 3 and 4 was also studied. Overall, it improves the LBGI in combination with both 15064



Article

Figure 9. Stage 2 comparison of control combinations.

A meal estimation error of ±25%, an insulin sensitivity (IS) bias of ±20%, and a varying circadian-rhythm-based insulin sensitivity (one of the three) were used to create scenarios 9−16, the combinations of which are provided in Table 6. The meal sizes and insulin boluses and the timings of the two were maintained same in all scenarios. The data loss and CGM sensor noise were also left unchanged for these scenarios. BGC regulation by the embedded controller under a +25% meal estimation error and insulin sensitivity biases for both patients is provided in Figure 10a. Also, a similar visualization of the impact of a −25% meal estimation error in the presence of insulin sensitivity biases in shown in Figure 10b. From a metrics stand point, the LBGI value signifies the ability of this control combination (scenario 8) to tackle meal estimation errors as well as IS uncertainties. Under all scenarios, a healthy minimum BGC that is above the 60 mg/dL lower bound is maintained. Also, a healthy average BGC level is maintained under tested scenarios 9−16. An HBGI measure of less than 5 ensures the avoidance of hyperglycemic risks. CVGA, accepted as a valuable tool for evaluating the effectiveness of BGC regulation, was performed.52 On the basis of the minimum and maximum BGC levels obtained during multiple meal disturbances for scenarios 8−16 for both virtual patients, a CVGA plot is shown in Figure 11. As determined by the LBGI and minimum BGC metrics, CVGA also confirms the absence of any hypoglycemic events.

Figure 9. The use of the IOB safety trigger is depicted in the insulin infusion shown in Figure 9. It reduces the current infusion rate to the basal inflow and ensures no further increase until the IOB criterion is overcome. This has a huge impact on the LBGI values, as shown in Tables 7 and 8. Although this results in an increase in the eAG value, it is still maintained within the safe BGC range (80−160 mg/dL). Hence, IOB allows the controllers to respond to a probable hypoglycemic condition and improves the LBGI measures even in the presence of uncertainties. Considering the safety factors, based on the metrics obtained, for scenarios 1−8, there is no risk of hyperglycemic exposure because all HBGI values fall below 5.1 Different combinations of modules are able to improve the performance of the controller with respect to the LBGI measure. Considering the lowest LBGI alongside improved controller performance with lower MSE, the combination used in scenario 8 (mpMPC + FLA + FF + BOL + IOB) is concluded to be the best possible combination of modules in the current study. All of the results obtained from scenarios 1−8 were obtained from the entire simulation run in MATLAB. The best combination, namely, scenario 8 was converted to C code, and further experiments based on uncertainties in insulin sensitivity and meal estimation, along with scenario 8, were carried out with implementation in an embedded device whose scheme is described in section 6. The meal sizes and timings are assumed to be informed by the patient without fail. The embedded microcontroller-based implementation was tested on both patients. 15065



Article

Figure 10. Embedded controller performance under meal estimation errors of (a) +25% and (b) −25%.

8. CONCLUSIONS

(CBA) is used in patient-specific model parameter estimation using MOGTT data. An input−output segregation of the blood glucose concentration (BGC) range followed by a gap-metric-based k-means clustering of multiple linear models is demonstrated. The gap metrics of linear models from cluster representatives are used in

In this work, mpMPC, an explicit MPC, and different combinations of the two are investigated. Two patients are derived from a nominal H−W type 1 diabetic by varying the insulin sensitivity of the nominal patient. The chaotic bat algorithm 15066



Article

The gut absorption rate of glucose is expressed as UG =

DGA Gt e−t / tmax,G tmax ,G 2

(A.7)

where the glucose appearance rate is highest at tmax,G minutes. The insulin absorption submodel is provided by the equations

Figure 11. Control variability grid analysis of minimum/maximum BGC for both virtual patients.

■

S2̇ = ka1 *S1a − ka2 *S2

(A.10)

S3̇ = ka1 *S2 + ka2 *S1b − keS3

(A.11)

S1a k m,LD + S1a

(A.12)

S1b k m,LD + S1b

(A.13)

ASSOCIATED CONTENT

S Supporting Information *

Details on the chaotic bat algorithm and k-means algorithm; comparison of existing insulin on board estimation and fuzzy aggregated linearized full-order H−W model based insulin on board estimation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +91431-2503104. Fax: +91431-2500133. E-mail: radha@ nitt.edu.

⎛ FC ⎞ Q̇ 1 = −⎜ 01 + x1⎟Q 1 + k12Q 2 − FR + UG ⎝ VGG ⎠

Notes

The authors declare no competing financial interest.

+ EGP0(1 − x3)

(A.1)

Q̇ 2 = x1Q 1 − (k12 + x 2)Q 2

(A.2)

G = Q 1/VG

(A.3)

■

ACKNOWLEDGMENTS V.K. is supported by the Technical Education Quality Improvement Programme (TEQIP) Phase II, a World Bank initiative. The authors also acknowledge the support of TEQIP-II and NIT, Tiruchirappalli, India, in this research work. Finally, the authors acknowledge the valuable feedback provided by anonymous reviewers that helped improve this article greatly.

where EGP0 represents endogenous glucose production and FC01 is the total noninsulin dependent glucose flux, given by

■

⎪

LIST OF SYMBOLS

General Symbols

(A.4)

A = discretized system matrix b = bat b B = discretized input matrix c = constraint set C = output matrix cF = correction factor CR = critical region D = optimal constraint set Dg = meal intake E = loudness f b = frequency of bat b f max = maximum frequency f min = minimum frequency F = inactive constraint set

Also, the renal clearance (FR) above the glucose threshold is given by

(A.5)

The three compartmental insulin action subsystems are given by f x1̇ = −ka1x1 + SIT ka1(S3/Vi ) f x 2̇ = −ka2x 2 + SID ka2(S3/Vi ) f x3̇ = −ka3x3 + SIE ka3(S3/Vi )

(A.9)

LDb = Vmax,LD

APPENDIX The H−W model equations are provided in this appendix. Glucose in accessible and nonaccessible compartments is given by

⎧ 0.003(G − 9)VG if G > 9 mmol/L FR = ⎨ ⎩0 otherwise

̇ = (1 − k ins)u − ka2 *S1b − LDb S1b

LDa = Vmax,LD

■

⎪

(A.8)

where

obtaining FLA weights. Linearized full-order models of CBA-estimated virtual patients in cluster representative regions are derived. A fuzzy composite of these models using FLA weights is used in calculating the IOB criterion. The CGM sensor errors and CGM sensor data loss are included during closed-loop experiments to reflect practical scenarios. The best combination of these modules developed (mpMPC + FLA + FF + IOB + BOL) is implemented in an embedded platform, and HIL testing of the controller under meal estimation errors and IS uncertainties is performed for both virtual patients. It is concluded that FLA greatly increases the transient performance of the artificial pancreas. Also, insulin on board (IOB) improves the safety of this artificial pancreas greatly by regulating insulin infusion as indicated by low LBGI values.

⎧ F01 if G > 4.5mmol/L C F01 =⎨ ⎩ F01G /4.5 otherwise

S1ȧ = k insu − ka1 *S1a − LDa

(A.6) 15067



Article

θ = CGM sensor error κ = gap metric λ = Lagrange multiplier μ = fuzzification function ξ = random number ρ = fuzzification factor σ = pulse decay rate τI,1, τI,2 = input−output transfer function time constants τd,I = input−output delay time τd,M = disturbance delay time τM,1, τM,2 = disturbance time constants φ = white Gaussian noise χ = predicted system states ω = frequency

g = gradient matrix G = measured blood glucose concentration G̃ = predicted blood glucose concentration G̅ = mean of experimental measurements h = insulin decay rate H = Hessian matrix i = iteration count I = constraint set k = sampling instant k12 = transfer rate ka1, ka2, ka3 = deactivation rates ka1* = slow-channel transfer rate ka2* = fast-channel transfer rate kb1, kb2, kb3 = activation rates kd = disturbance gain ke = insulin elimination rate kio = input−output gain kM,LD = half-concentration constant lb = lower bound L = matrix of constraint bounds LDa = local degradation of insulin in the slow channel LDb = local degradation of insulin in the fast channel M = matrix of constraint activation Nc = control horizon Ne = total number of samples Nn = total number of clusters Np = prediction horizon O = matrix of state constraints p = prediction count P = terminal cost pb = position of bat b Q = state weighting factor Q1 = mass of glucose in the accessible compartment Q2 = mass of glucose in the inaccessible compartment r = pulse rate R = input weighting factor s = risk function S = inverse of the fitness score S1a, S2 = slow-channel insulin S1b = fast-channel insulin S3 = plasma insulin SfID = insulin sensitivity of disposal SfIE = insulin sensitivity of endogenous glucose production SfIT = insulin sensitivity of transportation Ts = sampling time uacc = insulin accumulation ui = insulin infusion rate uiob = insulin on board ub = upper bound vb = velocity of bat b Vg = glucose distribution volume Vi = insulin distribution volume Vmax,LD = saturation level w = defuzzification weight X = mpMPC-transformed state constraint matrix x1 = action on glucose transport x2 = action on glucose elimination x3 = action on endogenous glucose production

Acronyms

■

AP = artificial pancreas API = application peripheral interface BA = bat algorithm BGC = blood glucose concentration BOL = bolus insulin CBA = chaotic bat algorithm CEG = Clarke’s error grid CGM = continuous glucose monitoring CHO = carbohydrate CVGA = control variability grid analysis DOE = design of experiments FF = feed forward FIR = finite impulse response FLA = fuzzy logic aggregation FOPDTPI = first-order plus delay time plus integrator GS = gain scheduling HAL = hardware application layer HBGI = high blood glucose index HIL = hardware in the loop H−W = Hovorka−Wilinska ICE = in-circuit emulator IOB = insulin on board IS = insulin sensitivity LBGI = low blood glucose index LICQ = linear inequality constraint qualification mpMPC = multiparametric model predictive controller MOGTT = modified oral glucose tolerance test MPC = model predictive controller MSE = mean square of error OGTT = oral glucose tolerance test PWA = piecewise affine SCS = strict complementary slackness SOPDT = second-order plus delay time T1DM = type 1 diabetes mellitus USB = universal serial bus

REFERENCES

(1) Liu, S. W.; Huang, H. P.; Lin, C. H.; Chien, I. L. A Hybrid Neural Network Model Predictive Control with Zone Penalty Weights for Type 1 Diabetes Mellitus. Ind. Eng. Chem. Res. 2012, 51, 9041. (2) Percival, M. W.; Bevier, W. C.; Wang, Y.; Dassau, E.; Zisser, H.; Jovanovič, L.; Doyle, F. J., III. Modeling the effects of subcutaneous insulin administration and carbohydrate consumption on blood glucose. J. Diabetes Sci. Technol. 2010, 4 (5), 1214. (3) Sorensen, J. A physiologic model of glucose metabolism in man and its use to design and assess improved insulin therapies for diabetes. Ph.D. Dissertation, Massachusetts Institute of Technology, Cambridge, MA, 1985.

Greek Symbols

α = loudness decay rate β = random vector γ, δ, ε, ζ = Johnson transformation parameters η = rate constant 15068



Article

(4) Hovorka, R.; Canonico, V.; Chassin, I.; Haueter, U.; MassiBenedetti, M.; Federici, M. O.; Pieber, T. R.; Schaller, H. C.; Schaupp, L.; Vering, T.; Wilinska, M. E. Nonlinear model predictive control of glucose concentration in subjects with type 1diabetes. Physiol. Meas. 2004, 25 (4), 905. (5) Wilinska, M. E.; Chassin, L. J.; Schaller, H. C.; Schaupp, L.; Pieber, T. R.; Hovorka, R. Insulin Kinetics in Type-1 Diabetes: Continuous and Bolus Delivery of Rapid Acting Insulin. IEEE Trans. Biomed. Eng. 2005, 52 (1), 3. (6) Dalla Man, C.; Rizza, R. A.; Cobelli, C. Meal simulation model of the glucose−insulin system. IEEE Trans. Biomed. Eng. 2007, 54 (10), 1740. (7) Bequette, B. W. Challenges and recent progress in the development of a closed-loop artificial pancreas. Annu. Rev. Control 2012, 36, 255. (8) Rmileh, A. A.; Gabin, W. G. Feedforward−feedback multiple predictive controllers for glucose regulation in type 1 diabetes. Comput. Methods Programs Biomed. 2010, 99, 113. (9) Percival, M. W.; Dassau, E.; Zisser, H.; Jovanovič, L.; Doyle, F. J., III. Practical Approach to Design and Implementation of a Control Algorithm in an Artificial Pancreatic Beta Cell. Ind. Eng. Chem. Res. 2009, 48, 6059. (10) Farmer, T. G., Jr.; Edgar, T. F.; Peppas, N. A. Effectiveness of Intravenous Infusion Algorithms for Glucose Control in Diabetic Patients Using Different Simulation Models. Ind. Eng. Chem. Res. 2009, 48, 4402. (11) Galvanin, F.; Barolo, M.; Macchietto, S.; Bezzo, F. Optimal Design of Clinical Tests for Identification of Physiological Models of Type 1 Diabetes Mellitus. Ind. Eng. Chem. Res. 2009, 48, 1989. (12) Yang, X. S. Nature-Inspired Metaheuristic Algorithms; Luniver Press: Bristol, U.K., 2008. (13) Yang, X. S. Bat algorithm: Literature review and applications. Int. J. Bio-Inspired Comput. 2013, 5 (3), 141. (14) Yang, X. S. A New Metaheuristic Bat-Inspired Algorithm. In Nature Inspired Cooperative Strategies for Optimization (NICSO 2010); Cruz, C., González, J. R., Krasnogor, N., Pelta, D. A., Terrazas, G., Eds.; Studies in Computational Intelligence Series; Springer: Berlin, 2010; Vol. 284, pp 65−74. (15) Gandomi, A. H.; Yang, X. S. Chaotic Bat Algorithm. J. Comput. Sci. 2014, 5, 224. (16) Percival, M. W.; Wang, Y.; Dassau, E.; Zisser, H. C.; Jovanovič, L.; Doyle, F. J., III. Development of a multi-parametric model predictive control algorithm for insulin delivery in type 1 diabetes mellitus using clinical parameters. J. Process Control 2011, 21, 391. (17) Oruklu, M. E.; Quinn, A. C. L.; Smith, D. Adaptive control strategy for regulation of blood glucose levels in patients with type 1 diabetes. J. Process Control 2009, 19, 1336. (18) Campetelli, G.; Lombarte, M.; Basualdo, M. S.; Rigalli, A. Extended adaptive predictive controller with robust filter to enhance blood glucose regulation in type I diabetic subjects. Comput. Chem. Eng. 2013, 59, 243. (19) Kirubakaran. V.; Radhakrishnan, T. K.; Sivakumaran, N. In Blood glucose concentration regulation in Type 1 diabetics using multi model multi parametric model predictive control: An empirical approach. In Proceedings of the 12th IFAC Symposium on Computer Applications in Biotechnology; Modak, D., Ed.; Elsevier Science Ltd.: Amsterdam, The Netherlands, 2013; pp 291−296. (20) Canete, J. F.; Perez, S. G.; Diaz, J. C. R. Artificial neural networks for closed loop control of in silico and ad hoc type 1 diabetes. Comput. Methods Programs Biomed. 2012, 106, 55. (21) Dua, P.; Doyle, F. J., III; Pistikopoulos, E. N. Model-Based Blood Glucose Control for Type 1 Diabetes via Parametric Programming. IEEE Trans. Biomed. Eng. 2006, 53 (8), 1478. (22) Pistikopoulos, E. N.; Dua, V.; Bozinis, N. A.; Bemporad, A.; Morari, M. On-line optimization via off-line parametric optimization tools. Comput. Chem. Eng. 2002, 26, 175. (23) Venkat, A.; Gudi, R. Fuzzy segregation based identification and control of nonlinear dynamic systems. Ind. Eng. Chem. Res. 2002, 41, 538.

(24) Hariprasad, K.; Bhartiya, S.; Gudi, R. A gap metric based multiple model approach for nonlinear switched systems. J. Process Control 2012, 22, 1743. (25) Kovatchev, B. P.; Breton, M.; Dalla Man, C.; Cobelli, C. In Silico Preclinical Trials: A Proof of Concept in Closed-Loop Control of Type 1 Diabetes. J. Diabetes Sci. Technol. 2009, 3 (1), 44. (26) Breton, M.; Kovatchev, B. Analysis, modeling, and simulation of accuracy of continuous glucose sensors. J. Diabetes Sci. Technol. 2008, 2 (5), 1122. (27) Russell, S. J.; El-Khatib, F. H.; Nathan, D. M.; Magyar, K. L.; Jiang, J.; Damiano, E. R. Blood Glucose Control in Type 1 Diabetes with a Bihormonal Bionic Endocrine Pancreas. Diabetes Care 2012, 35, 2148. (28) Hovorka, R. The future of continuous glucose monitoring: Closed loop. Curr. Diabetes Rev. 2008, 4 (3), 269. (29) Turksoy, K.; Bayrak, E. S.; Quinn, L.; Littlejohn, E.; Cinar, A. Multivariable Adaptive Closed-Loop Control of an Artificial Pancreas without Meal and Activity Announcement. Diabetes Technol. Therapeut. 2013, 15 (5), 386. (30) Patek, S. D.; Magni, L.; Dassau, E.; Karvetski, C. H.; Toffanin, C.; De Nicolao, G.; Favero, S. D.; Breton, M.; Dalla Man, C.; Renard, E.; Zisser, H.; Doyle, F. J., III; Cobelli, C.; Kovetchev, B. P. International Artificial Pancreas Study Group. Modular Closed-Loop Control of Diabetes. IEEE Trans. Biomed. Eng. 2012, 59 (11), 2986. (31) Kovatchev, B. P.; Renard, E.; Cobelli, C.; Zisser, H. C.; Hynes, P. K.; Anderson, S. M.; Brown, S. A.; Chernavvsky, D. R.; Breton, M.; Farret, A.; Pelletier, M. J.; Place, J.; Bruttomesso, D.; Favero, S. D.; Visentin, R.; Filippi, A.; Scotton, R.; Avogaro, A.; Doyle, F. J., III. Feasibility of Outpatient Fully Integrated Closed-Loop Control. Diabetes Care 2013, 36, 1851. (32) Aye, T.; Block, J. Toward Closing the Loop: An Update on Insulin Pumps and Continuous Glucose Monitoring Systems. Endocrinol. Metab. Clin. North Am. 2010, 39 (3), 609. (33) O’Connell, M. A.; Donath, S.; O’Neal, D. N.; Colman, P. G.; Ambler, G. R.; Jones, T. W.; Davis, E. A.; Cameron, F. J. Glycaemic impact of patient-led use of sensor-guided pump therapy in type 1 diabetes: A randomised controlled trial. Diabetologia 2009, 52, 1250. (34) Tamborlane, W.; Beck, R.; Laffel, L. Continuous Glucose Monitoring and Type 1 Diabetes. N. Engl. J. Med. 2009, 360, 191. (35) Parker, R. S.; Doyle, F. J., III; Ward, J. H.; Peppas, N. A. Robust H∞ glucose control in diabetes using a physiological model. AIChE J. 2000, 46, 2537. (36) Ramprasad, Y.; Rangaiah, G. P.; Lakshminarayanan, S. Enhanced IMC for Glucose Control in Type 1 Diabetics Using a Detailed Physiological Model. Food Bioprod. Proc. 2006, 84 (C3), 227. (37) Sambariya, D. K.; Prasad, R. Robust tuning of power system stabilizer for small signal stability enhancement using metaheuristic bat algorithm. Int. J. Electr. Power Energy Syst. 2014, 61, 229. (38) Rodrigues, D.; Pereira, L. A. M.; Nakamura, R. Y. M.; Costa, K. A. P.; Yang, X. S.; Souza, A. N.; Papa, J. P. A wrapper approach for feature selection based on Bat Algorithm and Optimum-Path Forest. Expert Syst. Appl. 2014, 41, 2250. (39) Hasançebi, O.; Teke, T.; Pekcan, O. A bat-inspired algorithm for structural optimization. Comput. Struct. 2013, 128, 77. (40) Clarke, W. L.; Cox, D.; Gonderfrederick, L. A.; Carter, W.; Pohl, S. L. Evaluating clinical accuracy of systems for self-monitoring of bloodglucose. Diabetes Care 1987, 10 (5), 622. (41) Rugh, W.; Shamma, J. Research on gain scheduling. Automatica 2000, 36 (10), 1401. (42) Chischi, L.; Falugi, P.; Zappa, G. Gain scheduling MPC for nonlinear systems. Int. J. Robust Nonlinear Control 2003, 13 (3), 295. (43) Kuipers, B.; Astrom, K. The composition and validation of heterogeneous control laws. Automatica 1994, 30, 233. (44) Narendra, K. S.; Balakrishnan, J.; Ciliz, M. K. Adaptation and learning using multiple models, switching and tuning. IEEE Control Syst. Mag. 1995, 15 (3), 37. (45) El-Sakkary, A. The gap metric: Robustness of stabilization of feedback systems. IEEE Trans. Automat. Control 1985, 30 (3), 240. (46) Zhou, K.; Doyle, J. Essentials of Robust Control; Prentice Hall PTR: Upper Saddle River, NJ, 1998. 15069



Article

(47) Venkat, A.; Vijaisai, P.; Gudi, R. Identification of complex nonlinear processes based on fuzzy decomposition of the steady state space. J. Process Control 2003, 13, 473. (48) Tondel, P.; Johansen, T. A.; Bemporad, A. An algorithm for multiparametric quadratic programming and explicit MPC solutions. Automatica 2003, 39, 489. (49) Kirubakaran, V.; Radhakrishnan, T. K.; Sivakumaran, N. Distributed multiparametric model predictive control design for a quadruple tank process. Measurement 2014, 47, 841. (50) Gupta, A.; Bhartiya, S.; Nataraj, P. S. V. A novel approach to multiparametric quadratic programming. Automatica 2011, 47, 2112. (51) Walsh, J.; Roberts, R. Pumping Insulin; Torrey Pines Press: San Diego, CA, 2006. (52) Magni, L.; Raimondo, D. M.; Dalla Man, C.; Breton, M.; Patek, S.; Nicolao, G. D.; Cobelli, C.; Kovatchev, B. P. Evaluating the efficacy of closed-loop glucose regulation via control-variability grid analysis. J. Diabetes Sci. Technol. 2008, 2 (4), 630.

15070


Metaheuristic Patient Estimation Based Patient-Specific Fuzzy

Recommend Documents