Economic Model Predictive Control of Bihormonal Artificial Pancreas System Based on Switching Control and Dynamic R-parameter

Abstract

Background:

Blood glucose (BG) regulation is a long-term task for people with diabetes. In recent years, more and more researchers have attempted to achieve automated regulation of BG using automatic control algorithms, called the artificial pancreas (AP) system. In clinical practice, it is equally important to guarantee the treatment effect and reduce the treatment costs. The main motivation of this study is to reduce the cure burden.

Methods:

The dynamic R-parameter economic model predictive control (R-EMPC) is chosen to regulate the delivery rates of exogenous hormones (insulin and glucagon). It uses particle swarm optimization (PSO) to optimize the economic cost function and the switching logic between insulin delivery and glucagon delivery is designed based on switching control theory.

Results:

The proposed method is first tested on the standard subject; the result is compared with the switching PID and the switching MPC. The effect of the dynamic R-parameter on improving the control performance is illustrated by comparing the results of the EMPC and the R-EMPC. Finally, the robustness tests on meal change (size and timing), hormone sensitivity (insulin and glucagon), and subject variability are performed. All results show that the proposed method can improve the control performance and reduce the economic costs.

Conclusions:

The simulation results verify the effectiveness of the proposed algorithm on improving the tracking performance, enhancing robustness, and reducing economic costs. The method proposed in this study owns great worth in practical application.

Keywords

artificial pancreas dynamic R-parameter economic cost economic model predictive control switching control

Type 1 diabetes (T1D) is one of many incurable metabolic diseases with the symptom that the body only secretes a negligible amount of insulin because of the destruction of the pancreatic beta-cells. The blood glucose (BG), for people with diabetes, always lies at abnormal levels, which causes a lot of complications, such as nephropathy, peripheral neuropathy, blindness, and so on. In 2015, the International Diabetes Federation (IDF) published the newest Diabetes Atlas, that is, the IDF Diabetes Atlas (7th ed.), which showed that a total of 5 million people died of diabetes by the end of 2015, and the world’s diabetes treatment costs exceeded $673 billion.¹ The rapid increase of people with diabetes and the serious results that causes in health care and economic costs all promote the development of diabetes therapies.

The most common therapy should be multiple daily injections of insulin to regulate the glucose level. In recent decades, more and more researchers, not in medical field, engage in diabetes treatment. For example, automatic control scholars attempt to combine continuous glucose monitoring, continuous subcutaneous insulin infusion and an automatic control algorithm to achieve closed-loop control of BG, that is, artificial pancreas (AP). The AP can be divided into two types: unihormonal AP and bihormonal AP. Several algorithms have been used in unihormonal AP, such as PID in Marchetti et al,² model predictive control (MPC) in Wang and Doyle³ and Wang et al,⁴ and fuzzy logic control in Grant.⁵ The previous studies showed that hypoglycemia can cause more direct damage than hyperglycemia and can even cause death.⁶ To avoid this risk, researchers add the glucagon infusion into the system to transfer unihormonal AP into bihormonal AP. There are also many algorithms using in this type AP, such as the switching PID control in Gao and Wang⁷ and the switching MPC in Ning and Wang.⁸

In clinical practice, the primary task is to achieve an effective treatment that guarantees the safety of the patients; however, the decrease of treatment costs is also important. This article fills the gap of previous studies that pay no attention to decreasing the treatment costs. In this article, to reduce the treatment costs, the EMPC is chosen to be the core algorithm. The EMPC originated from traditional two-layered real-time optimization with a similar framework to conventional MPC,⁹ which focused more on the economic consideration using a general cost function which differs from the quadratic cost function of MPC.¹⁰ Hence, the EMPC inherits the merits of the MPC, which can manage time delays, constraints, and multivariate effectively. These advantages all prove that the EMPC is a good choice to regulate the glucose in the AP because the subject model system is a nonlinear time-varying model with great time delays and strong disturbances. The design of switching rules also refers to the core idea of the EMPC, which can reduce the economic cost further.

The dynamic R-parameter¹¹ is utilized to improve the control performance, make zero-error tracking possible, which is designed based on the tracking error of the last sampling moment. For the insulin subsystem and the glucagon subsystem, the design of (the R-parameter) is different because they have opposite effect on glucose regulation. Several penalty parameters are used to adjust the value of r.

Simulation results are shown in figures and the statistical data are recorded in tables. After the simulation on the standard subject, the comparison for tracking performance and the tests for robustness, we can draw the conclusion that the proposed method can improve tracking performance and reduce treatment costs at the same time. For the standard subject, BG can be maintained within the normal range (70-180 mg/dl), the average tracking error (ATE) and BG risk index (BGRI) are 14.04 and 1.01, respectively, which is sufficient. In addition, the economic costs of insulin and glucagon are reduced by 41.74% and 70.40% compared with the switching PID,⁷ 28.47% and 53.92% compared with the switching MPC.⁸

Simulation Model

A meal model, an insulin model, a glucagon model and a glucose model form the whole simulation model. In this study, the meal model and the glucose model are selected from Hovorka et al,¹² the insulin model can be seen in Wilinska et al,¹³ and the glucagon model is proposed in Gao and Wang.⁷ All these models are authoritative and have been used widely in AP study.

We name the virtual subject with the same parameters of Gao and Wang,⁷ Hovorka et al,¹² and Wilinska et al¹³ as the standard subject. Another nine virtual subjects (numbered 1-9) with different parameters obtained by the same method in Ning and Wang⁸ will be used to make the robustness tests on subject variability.

The Switching Dynamic R-parameter Economic Model Predictive Control

The framework of the switching R-EMPC is shown in Figure 1.

Figure 1.

The structure chart of the switching R-EMPC.

Dynamic R-parameter Economic Model Predictive Control

The detailed overview of the EMPC can be found in Ellis et al.¹⁴ The key procedures include the identification of the prediction model, the design of the cost function, and the receding horizon optimization of inputs.

The prediction models are used to approximate the relationship between insulin and glucose and glucagon and glucose, respectively. In this study, the autoregressive exogenous (ARX) model is used, which is described as:

\begin{array}{l} A (z^{- 1}) G (t + 1) = B (z^{- 1}) I (t - d) + w (t) \\ A (z^{- 1}) = 1 + a_{1} z^{- 1} + \dots + a_{n} z^{- n} \\ B (z^{- 1}) = b_{0} + b_{1} z^{- 1} + \dots + b_{m} z^{- m} \end{array}

where $z^{- 1}$ is the backward shift operator; d is the time delay; w(t) denotes disturbances. All parameters can be obtained after a 24-hour open-loop step response experiment.

For the design of the cost function, we add the idea of the dynamic R-parameter¹⁵ into the general EMPC to improve the tracking performance; hence, the cost function of each subsystem can be designed as follows:

\begin{array}{l} Ω \hat{=} \sum_{i = 1}^{N} α_{1} {(\hat{G} (t + i | t) - y_{r} (t + i))}^{2} \\ + r \times \sum_{j = 0}^{M - 1} (α_{2} U (t + j | t) + α_{3} | Δ U (t + j | t) |) \end{array}

where the second term aims to reduce the usage of hormone directly; N and M (N>M) are the prediction horizon and the control horizon, chosen as 30 and 3 respectively; $\hat{G}$ denotes the predictive output; $y_{r}$ is the set point, selected as 110 mg/dl; $U (t + j | t) |_{j = 0}^{M - 1}$ is the predictive input (insulin or glucagon); $Δ U (t + j | t) |_{j = 0}^{M - 1}$ is the change rate of input; $α_{i} (i = 1, 2, 3)$ is the penalty coefficient, chosen as {1, 22.5, 15} for the insulin subsystem and {1, 75, 15} for the glucagon subsystem. And r is the dynamic R-parameter, which is designed based on the tracking error of the last sampling moment. r of the insulin subsystem and the glucagon subsystem are designed as (3) and (4), respectively:

r = {\begin{cases} τ_{1} (1 - \exp (τ_{2} e (t - 1))), e (t - 1) > 0 \\ \exp (τ_{3} e (t - 1)), e (t - 1) < 0 \end{cases}

r = {\begin{cases} \exp (τ_{4} e (t - 1)), e (t - 1) > 0 \\ τ_{5} (1 - \exp (τ_{6} e (t - 1))), e (t - 1) < 0 \end{cases}

where $e (t - 1) = y_{r} - G (t - 1)$ is the tracking error of the last sampling moment and $τ_{i} (i = 1, 2, 3, 4, 5, 6)$ are the penalty parameters, chosen as {15, -0.005, 0.03, -0.0005, 1.03, 15}.

In this study, particle swarm optimization (PSO) is chosen to solve the receding horizon optimization problem, the predictive inputs signal can be obtained by solving the optimization problem (5).

U^{*} (t + j | t) |_{j = 0}^{M - 1} = \underset{U (t + j | t)}{\arg \min} Ω

Then, only the first input $U (t | t)$ is used in the control, the optimization will be repeated in the next sampling moment.

The Switching Rules

In this article, the whole system is divided into three modes: the insulin infusion mode, the glucagon infusion mode, and the zero-input mode. The switching rules are used to simulate the switching logic between them.

The switching cost function is designed like (6) considering the tracking error and the economic cost comprehensively:

\begin{array}{l} J_{δ} \hat{=} \sum_{i = 1}^{N_{1}} {(\hat{G} (t + i | t) - y_{r} (t + i))}^{2} + \sum_{j = 0}^{M_{1} - 1} α_{δ} {(U (t + j | t))}^{2} {}_{δ}\in {1, 2} \\ J_{δ} \hat{=} \sum_{i = 1}^{N_{1}} {(\hat{G} (t + i | t) - y_{r} (t + i))}^{2} {}_{δ}\in {3} \end{array}

where the same notations has the same meanings as equation (2), and M₁ and N₁ represent the control horizon and prediction horizon, which are chosen as 3 and 10, respectively, of the switching rules (M₁ ≤ M, N₁ ≤ N). $(M_{1} \leq M, N_{1} \leq N) .$

At every sample t, the appropriate mode can be obtained by optimizing (7).

{}_{δ}{(t)} = \underset{{}_{δ}\in {1, 2, 3}}{\arg \min J_{δ}}

The Simulation Results

Performance Assessment Indices

In this study, the safe range of BG concentration is defined as 70-180 mg/dl, BG greater than 180 mg/dl is called hyperglycemia, and BG lower than 70 mg/dl is called hypoglycemia.

The ATE is used to show the tracking performance of the controller, which is introduced as follows:

A T E (k) \hat{=} \sum_{t = (k - 1) T + 1}^{k T} | G (t) - y_{r} | / T

where k represents the k-th day, $G (t)$ is the real output glucose, $y_{r}$ is the set point, and $T$ is the number of samplings for one day, chosen as 288.

The BGRI¹⁶ are chosen to reflect the risk of BG.

Regarding the economic costs, the amount of insulin and glucagon used in simulations are recorded and compared.

The control variability grid analysis (CVGA)¹⁷ is used to visualize the quality of the control algorithm on different subjects.

Simulation on the Standard Subject

The duration of simulation is set to 2 days, the basal and bolus open-loop therapy is used in the first day, while the closed-loop control is implemented in the second day. During the simulation, patients eat 40 g, 60 g, and 85 g carbohydrates at 7:00, 12:00, and 18:00, respectively, which correspond to the patients’ three meals every day. These parameter settings accord with the real cases. Simulation does its best to imitate clinical test. In this simulation, the superiority of the switching R-EMPC is compared with the other two closed-loop control methods, that is, the switching PID⁷ and the switching MPC.⁸ The simulation results are shown in Figures 2 to 4 and the statistical data are recorded in Tables 1 and 2.

Figure 2.

Simulation results for the switching R-EMPC.

Figure 3.

Simulation results for the switching PID.⁷

Figure 4.

Simulation results for the switching MPC.⁸

Table 1.

Control Performance Statistics.

Therapy	% of BG	% of BG	ATE	BGRI
Therapy	< 70 mg/dl	> 180 mg/dl	(mg/dl)	BGRI
Switching PID⁷	0	0	31.89	2.13
Switching MPC⁸	0	0	22.94	1.44
Switching R-EMPC	0	0	14.04	1.01

Table 2.

Economic Cost Comparisons.

Therapy	Insulin amount (U)	Glucagon amount (U)
Switching PID⁷	74.46	13.31
Switching MPC⁸	60.65	8.55
Switching R-EMPC	43.38	3.94
Decrease than PID/MPC	41.74%/28.47%	70.40%/53.92%

The simulation results in Figures 2 to 4 and the statistical data in Table 1 show that the control performance of the switching R-EMPC is the best among these three methods. In addition, from the data in Table 2, compared with the switching PID,⁷ the switching R-EMPC decreases the economic costs of insulin and glucagon by 41.74% and 70.40%, respectively; compared with the switching MPC,⁸ it decreases the economic costs of insulin and glucagon 28.47% and 53.92%, respectively. These results prove that the proposed switching R-EMPC can improve control performance and decrease economic costs simultaneously.

Tracking Performance Improvement of Dynamic R-parameter

The simulation results of the switching R-EMPC are compared with the switching EMPC to illustrate the effectiveness on improving the tracking performance of the dynamic R-parameter, shown in Figure 5.

Figure 5.

Comparison of simulation results between the switching EMPC and the switching R-EMPC.

From Figure 5, the glucose control curve of the switching R-EMPC is smoother than that of the switching EMPC. Besides, the ATE and the BGRI of the switching R-EMPC are smaller than these of the switching EMPC (14.04 and 1.01 vs 18.25 and 1.09). Thus, the comparison results show that the dynamic R-parameter can improve the tracking performance.

Robustness Tests

In this section, simulation is designed to test the robustness of the proposed method with respect to inevitable variations in the real cases. First, the change on meal size and meal time occur frequently, in real life. It is necessary to test the robustness on meal change, we add a uniform distributed random variable to meal size, meal time, or both to realize meal changes in this case. Second, the hormone sensitivity of the same subject could be time-varying due to exercise, stress, and other factors. It is significant to test the robustness on hormone sensitivity; we multiply a uniform distributed sensitivity gain on the input of the insulin subsystem, the glucagon subsystem, or both to realize hormone sensitivity changes in this case. Finally, superior performance must be shown across a range of subjects for the results to be generally applied. So, the robustness tests on subject variability are performed.

The robustness tests on the meal change consists of three parts: the robustness tests on meal size change only, with the results shown in Figure 6; the robustness tests on meal time change only, with the results shown in Figure 7; and the robustness tests on combined changes, with the results shown in Figure 8. In addition, the statistical data about control performance and economic costs are recorded in Tables 3 and 4.

Figure 6.

Simulation results for the standard subject with ±50% meal size change.

Figure 7.

Simulation results for the standard subject with ±40 minutes meal time change.

Figure 8.

Simulation results on the standard subject with ±50% and ±40 minutes meal change.

Table 3.

Statistical Results for Control Performance With Meal Changes.

	Therapy	% of BG	% of BG	ATE	BGRI
	Therapy	<70 mg/dl	>180 mg/dl	(mg/dl)	BGRI
Meal size ±50%	Switching PID⁷	1.04	0	34.88	2.53
	Switching MPC⁸	2.43	0	23.37	1.63
	Switching R-EMPC	0	0	13.34	1.14
Meal time ±40 minutes	Switching PID⁷	0	1.39	29.64	2.10
	Switching MPC⁸	3.82	0	22.88	1.69
	Switching R-EMPC	0	0	13.62	1.22
Meal size and time ±50% and ±40 minutes	Switching PID⁷	0	0	28.91	1.80
	Switching MPC⁸	2.08	0	21.85	1.44
	Switching R-EMPC	0	0	12.25	0.83

Table 4.

Statistical Results for Economic Costs With Meal Changes.

	Therapy	Amount insulin (U)	Amount glucagon (U)
Meal size ±50%	Switching PID⁷	80.33	14.40
	Switching MPC⁸	63.05	8.60
	Switching R-EMPC	44.69	4.08
Meal time ±40 minutes	Switching PID⁷	70.12	12.91
	Switching MPC⁸	58.27	7.90
	Switching R-EMPC	51.37	4.88
Meal size and time ±50% and ±40 minutes	Switching PID⁷	68.68	11.94
	Switching MPC⁸	55.11	6.53
	Switching R-EMPC	54.05	5.50

According to Figures 6 to 8 and Tables 3 and 4, we can conclude that the robustness of the switching R-EMPC on meal change is better than the other two methods, and this method can significantly reduce the economic cost.

Similar to the robustness tests on meal change, the robustness of the R-EMPC on insulin sensitivity and glucagon sensitivity are considered separately at first, then in combination. The results are shown in Figures 9 to 11, with the statistical data on control performance and economic costs recorded in Tables 5 and 6.

Figure 9.

Simulation results for the standard subject with ±20% insulin sensitivity.

Figure 10.

Simulation results for the standard subject with ±20% glucagon sensitivity.

Figure 11.

Simulation results for the standard subject with ±20% insulin and glucagon sensitivity.

Table 5.

Statistical Results for Control Performance With Hormone Sensitivity Change.

	Therapy	% of BG	% of BG	ATE	BGRI
	Therapy	< 70 mg/dl	>180 mg/dl	(mg/dl)	BGRI
Insulin sensitivity ±20%	Switching PID⁷	0	2.78	27.12	2.06
	Switching MPC⁸	6.94	0	19.47	2.16
	Switching R-EMPC	0	0	9.82	0.63
Glucagon sensitivity ±20%	Switching PID⁷	0	0	31.41	2.10
	Switching MPC⁸	0	0	22.18	1.34
	Switching R-EMPC	0	0	10.85	0.71
Insulin and glucagon ±20%	Switching PID⁷	0	2.78	30.98	2.33
	Switching MPC⁸	0	0	18.00	1.54
	Switching R-EMPC	0	0	13.05	0.95

Table 6.

Statistical Results for the Economic Costs With Hormone Sensitivity Change.

	Therapy	Amount insulin (U)	Amount glucagon (U)
Insulin sensitivity ±20%	Switching PID⁷	85.97	17.22
	Switching MPC⁸	86.09	13.41
	Switching R-EMPC	41.62	3.78
Glucagon sensitivity ±20%	Switching PID⁷	73.09	12.98
	Switching MPC⁸	60.11	8.27
	Switching R-EMPC	57.51	5.75
Insulin and glucagon ±20%	Switching PID⁷	85.16	16.25
	Switching MPC⁸	85.61	13.14
	Switching R-EMPC	67.16	8.80

From Figures 9 to 11 and Tables 5 and 6, the robustness of the switching R-EMPC on hormone sensitivity is superior to the other two methods and consumes less insulin and glucagon.

Regarding the robustness test on subject variability, the switching R-EMPC proposed in this study is tested on 10 different subjects, previously mentioned in the Section of Simulation Model. The glucose control curves of the 10 subjects are shown in Figure 12, where the subfigure (a) corresponds to the standard subject and the subfigures (b)-(j) correspond to subjects 1 to 9, respectively. Figure 13 is the CVGA of them.

Figure 12.

Simulation results for the 10 subjects.

Figure 13.

The control variability grid analysis (CVGA) for the 10 subjects.

In Figure 13, as all points fall in zone B; combined with the results shown in Figure 12, we have reason to believe that the method proposed in this study has excellent robustness on subject variability.

Conclusion

The economic model predictive control, dynamic R-parameter, and switching theory, are combined to design the bihormonal AP system in this study. Simulation results and statistical data show that the proposed method has excellent control performance and robustness. The dynamic R-parameter can improve the tracking performance effectively. In addition, the switching R-EMPC performs well in reducing economic costs. Simulation scenarios in the Simulation Results are designed to imitate the clinical trials. Simulation results demonstrate that the proposed method owns the potential of clinical insight. The proposed method will be clinically tested in the near future.

Footnotes

Abbreviations

AP, artificial pancreas; ARX, autoregressive exogenous; ATE, average tracking error; BG, blood glucose; BGRI, blood glucose risk index; CVGA, control variability grid analysis; EMPC, economic model predictive control; IDF, International Diabetes Federation; MPC, model predictive control; PID, proportional integral derivative control; PSO, particle swarm optimization; R-EMPC, dynamic R-parameter economic model predictive control; T1D, type 1 diabetes.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: National Natural Science Foundation of China (61374099) and Research Fund for the Taishan Scholar Project of Shandong Province of China.

References

International Diabetes Federation. IDF Diabetes Atlas. 7th ed. 2015. Available at: http://www.idf.org.

Marchetti

Barolo

Jovanovič

Zisser

Seborg

DE.

An improved PID switching control strategy for type 1 diabetes. IEEE Trans Biomed Eng. 2008;55:857-865.

Wang

Doyle

III . Indirect iterative learning control: application on artificial pancreatic β-cell. Presentation at: 2009 IEEE Control and Decision Conference. Guilin, China; June 17-19, 2009:1739-1744.

Wang

Zisser

Dassau

Jovanovič

Doyle

III . Model predictive control with learning-type set-point: application to artificial pancreatic β-cell. AIChE J. 2010;56:1510-1518.

Grant

A new approach to diabetic control: fuzzy logic and insulin pump technology. Med Eng Phys. 2007;29:824-827.

Bloomgarden

ZT.

Diabetes complications. Diabetes Care. 2004;27:1506-1514.

Gao

Wang

. Closed-loop blood glucose control using dual subcutaneous infusion of insulin and glucagon based on switching PID controller. In: Proceedings of the 10th World Congress on Intelligent Control and Automation. Beijing, China; July 6-8, 2012:5023-5029.

Ning

Wang

. Bihormonal artificial pancreas system based on switching model predictive control. In: Proceedings of the 34th Chinese Control Conference (CCC). Hangzhou, China; July 28-30, 2015:4156-4161.

Jing

Liu

Economic model predictive control with triggered evaluations: state and output feedback. J Process Contr. 2014;24:1197-1206.

10.

Ellis

Christofides

PD.

On finite-time and infinite-time cost improvement of economic model predictive control for nonlinear systems. Automatica. 2014;50:2561-2569.

11.

Jia

Yang

Chiu

An integrated iterative learning control strategy with model identification and dynamic R-parameter for batch processes. J Process Contr. 2013;23:1332-1341.

12.

Hovorka

Canonico

Chassin

et al . Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiol Meas. 2004;25:905-920.

13.

Wilinska

Chassin

Schaller

Schaupp

Pieber

Hovorka

Insulin kinetics in type-1 diabetes: continuous and bolus delivery of rapid acting insulin. IEEE Trans Biomed Eng. 2005;52:3-12.

14.

Ellis

Durand

Christofides

PD.

A tutorial review of economic model predictive control methods. J Process Contr. 2014;24:1156-1178.

15.

Jia

Han

Chiu

MS.

Dynamic R-parameter based integrated model predictive iterative learning control for batch processes. J Process Contr. 2017;49:26-35.

16.

Clarke

Kovatchev

Statistical tools to analyze continuous glucose monitor data. Diabetes Technol Ther. 2009;11(suppl 1):S45-S54.

17.

Magni

Raimondo

Man

et al . Evaluating the efficacy of closed-loop glucose regulation via control-variability grid analysis. J Diabetes Sci Technol. 2008;2:630-635.