孫 超，孔雪華，高 軍，陸凱旋，杜鵬翔

(燕山大學(xué) 河北省測試計量技術(shù)及儀器重點實驗室，河北 秦皇島 066004)

1 引 言

隨著物質(zhì)生活水平的提高和環(huán)境的變化，糖尿病日益成為嚴重危害人類健康的多發(fā)疾病。據(jù)統(tǒng)計，我國糖尿病患者人數(shù)達到1.144億，已成為世界第一糖尿病大國，發(fā)病形式相當嚴峻。目前醫(yī)學(xué)界還沒有根治糖尿病的手段，胰島素治療是控制高血糖的重要手段，更是維持患者生命的唯一手段，人工胰臟系統(tǒng)被認為是治療糖尿病最有前景的方法之一[1]。

雙激素人工胰臟模擬了人體胰島對血糖的生理調(diào)節(jié)方式，是最有希望控制血糖水平的辦法。Damiano[2]提出了“雙皮下注射”的理念，由模型預(yù)測控制(model predictive ontrol,MPC)和比例-積分-微分(proportional integral derivative,PID)分別控制胰島素和胰高血糖素的注射速度[3]，減少低血糖的發(fā)生; Castle[4]用兩個衰減記憶比例-微分(fading memor proportional differential,FMPD)控制器分別控制胰島素和胰高血糖素的輸注速度，證明了雙激素血糖調(diào)節(jié)方法的臨床可行性。但上述兩項研究均采用兩個獨立的控制器，易造成胰島素和胰高血糖素的同時輸注?；谀Ｐ皖A(yù)測控制(MPC)的雙激素切換控制方法[5]、基于內(nèi)模的PID切換控制器(IMC-PID)[6,7]和雙激素人工胰臟的經(jīng)濟模型預(yù)測控制[8]，設(shè)計了輸入胰島素、輸入胰高血糖素和零輸入(兩種激素都不輸入)3種子模式，利用博弈論確定切換律，在3種模式中切換以達到控制效果的最優(yōu)，避免了兩種激素同時輸注的可能，但現(xiàn)有研究在運行效率和群體控制質(zhì)量上仍有待提高。

目前，關(guān)于雙激素人工胰臟子系統(tǒng)切換機制的研究大多存在運行效率和群體控制質(zhì)量較低的問題。為了進一步提高雙激素人工胰臟的控制效果和整體效率，本文在三支決策理論和模型預(yù)測控制算法的基礎(chǔ)上，提出一種基于三支決策的雙激素人工胰臟模型預(yù)測控制算法。首先，該算法采用現(xiàn)有的人工胰臟預(yù)測模型并對其組合方式進行了一定的改進，3個子系統(tǒng)的預(yù)測模型均能獨立得到各自的血糖預(yù)測值，也避免了進行大量的開環(huán)響應(yīng)實驗; 其次，該算法以三支決策理論為依托設(shè)計3個子系統(tǒng)間的切換規(guī)則，提供了一種切實可行的雙激素切換規(guī)則，以解決子系統(tǒng)切換的決策機制問題;最后利用美國弗吉尼亞大學(xué)和意大利帕多瓦大學(xué)聯(lián)合開發(fā)的UVA仿真平臺完成仿真，在血糖調(diào)節(jié)效果、跟蹤誤差、群體控制質(zhì)量和算法運行效率上取得了較好的結(jié)果。

2 基于三支決策的雙激素人工胰臟模型預(yù)測控制算法研究

基于三支決策的雙激素人工胰臟模型預(yù)測控制算法主要有模型預(yù)測控制和三支決策兩大理論基礎(chǔ)。首先，分別建立胰島素子系統(tǒng)、零輸注子系統(tǒng)和胰高血糖素子系統(tǒng)的預(yù)測模型，并用ode45算法進行求解得到各自的血糖預(yù)測值；其次，設(shè)計各個子系統(tǒng)的目標函數(shù)，并用粒子群算法進行優(yōu)化得到最優(yōu)控制輸入序列，選擇第1項作為當前最優(yōu)激素輸注速率；最后，根據(jù)最優(yōu)控制輸入、血糖預(yù)測值和輸注激素的價格因素，設(shè)計三支決策切換規(guī)則，即設(shè)計代價矩陣，分別計算切換到每個子系統(tǒng)面臨的風(fēng)險值，選擇切換到風(fēng)險值最小的子系統(tǒng)，將其激素輸注速率作為雙激素人工胰臟模型預(yù)測控制算法的輸出。本算法提供了一種切實可行的雙激素人工胰臟子系統(tǒng)切換機制，算法結(jié)構(gòu)框圖如圖1所示。

圖1 基于三支決策的雙激素人工胰臟模型預(yù)測控制算法結(jié)構(gòu)框圖Fig.1 Structural block diagram of model predictive control algorithm for double-hormone artificial pancreas based on three-branch decision theory

2.1 建立預(yù)測模型

本算法包含胰島素子系統(tǒng)、零輸注子系統(tǒng)和胰高血糖素子系統(tǒng)，分別建立3個子系統(tǒng)的預(yù)測模型并利用ode45算法求解得到各個子系統(tǒng)的血糖預(yù)測值。假設(shè)u1代表胰島素輸注速率；u2代表零輸注子系統(tǒng)的輸入，值恒為0；u3代表胰高血糖素輸注速率；y1代表胰島素子系統(tǒng)的血糖預(yù)測值；y2代表零輸注子系統(tǒng)的血糖預(yù)測值；y3代表胰高血糖素子系統(tǒng)的血糖預(yù)測值。

胰島素-血糖模型選用應(yīng)用較為廣泛的Bergman最小化模型[9～11]。Bergman最小化模型綜合考慮了血糖的降解、胰島素的降解、胰島素對機體吸收血糖速率、肝糖原輸出速率等因素對血糖值的影響，由3個微分方程構(gòu)成的微分方程組表示，公式為：

(1)

式中：t為時間；Gb、Ib分別為基礎(chǔ)狀態(tài)下的血糖、血漿胰島素的濃度；I(t)為血漿胰島素濃度與其基礎(chǔ)值的差值，μU/ml；X(t)為細胞組織間隙胰島素濃度；P1為在基礎(chǔ)胰島素水平下，葡萄糖自身代謝和抑制肝糖輸出的效能，min-1；P2為細胞組織間隙胰島素自身降解系數(shù)，min-1；P3為組織中胰島素代謝調(diào)控作用的系數(shù)，min-2·(μU/ml)-1；n為血漿胰島素自身降解系數(shù)，min-1；V1為體積，L。

零輸注子系統(tǒng)的模型采用改進的最小葡萄糖動力模型，公式如下：

(2)

式中：t為時間；Gb為基礎(chǔ)狀態(tài)下的血糖的濃度；SG為部分葡萄糖效果，min-1。

關(guān)于胰高血糖素-血糖模型的研究相對較少，胰高血糖素子系統(tǒng)的預(yù)測模型是將關(guān)于胰高血糖素輸注速率-血漿胰高血糖素濃度的吸收模型和關(guān)于血漿胰高血糖素濃度-血糖值的葡萄糖模型進行組合構(gòu)成復(fù)合模型[12,13]得到的。該復(fù)合模型是由5個微分方程構(gòu)成的微分方程組，如式(3)所示。

(3)

式中：t為時間；Nb為基礎(chǔ)狀態(tài)下的血漿胰高血糖素濃度；N(t)為血漿胰高血糖素的濃度，pg/mL；Y(t)為胰高血糖素對葡萄糖產(chǎn)生的作用，min-1；Z1(t)、Z2(t)為皮下注射胰高血糖素的兩室吸收，pg/kg；SG為部分葡萄糖效果，min-1；SN為胰高血糖素的敏感性，min-1/(pg/dL)；P4為描述胰高血糖素作用的時間常數(shù)的倒數(shù)，min-1；kN為胰高血糖素清除速率，min-1；tN為胰高血糖素吸收時間常數(shù)，min；VN為血漿胰高血糖素分布體積，mL/kg。

2.2 目標函數(shù)的設(shè)計與求解

由于零輸注子系統(tǒng)的控制輸入始終為零，因此無需設(shè)計零輸注子系統(tǒng)的目標函數(shù)。血糖調(diào)節(jié)的目的是將血糖控制在[70,180]的區(qū)間內(nèi)，為進一步提高控制效果，在區(qū)間控制的前提下進一步調(diào)節(jié)血糖使其接近設(shè)定值。由于低血糖風(fēng)險遠高于高血糖風(fēng)險，因此取正常區(qū)間中間偏上值130 mg/dL為設(shè)定值。設(shè)計包含設(shè)定值控制和區(qū)間控制2方面的目標函數(shù)，胰島素子系統(tǒng)和胰高血糖素子系統(tǒng)的目標函數(shù)包括控制輸入、血糖預(yù)測值與設(shè)定值偏差和血糖預(yù)測值與區(qū)間偏差3部分，如式(4)和式(5)所示。

(4)

(5)

s.t. Δumin≤Δu(k+i)≤Δumax

umin≤u(k+i)≤umax, (i=0,1,…，M-1)

式中：J1、J3分別是所設(shè)計的胰島素子系統(tǒng)、胰高血糖素子系統(tǒng)的目標函數(shù)；y1(t)、y3(t)分別表示胰島素子系統(tǒng)和胰高血糖素子系統(tǒng)的血糖預(yù)測值；ysp(t)=130 mg/dL為血糖設(shè)定值；ΔU(t)為最優(yōu)控制輸入增量；u(k+i)為系統(tǒng)的最優(yōu)控制輸入；R1、R3，Q1、Q3，S1、S3分別為目標函數(shù)3部分對應(yīng)的權(quán)值。

區(qū)間偏差部分是一個分段函數(shù)，當血糖預(yù)測值高于區(qū)間上限時，偏差值是血糖預(yù)測值與區(qū)間上限間的距離；當血糖預(yù)測值處于正常區(qū)間范圍內(nèi)時，偏差值始終為0；當血糖預(yù)測值低于區(qū)間下限時，偏差值是血糖預(yù)測值與區(qū)間下限間的距離。因此區(qū)間偏差值越小越好，與對目標函數(shù)進行最小優(yōu)化一致，可作為目標函數(shù)的一部分。利用粒子群優(yōu)化算法對2個目標函數(shù)進行求解得到2個子系統(tǒng)的最優(yōu)控制輸入序列，選擇第1項作為當前時刻的最優(yōu)輸入值即胰島素和胰高血糖素輸注速率。

2.3 三支決策切換規(guī)則

利用三支決策理論進行決策，包含兩個狀態(tài)集和3個行動集，并根據(jù)貝葉斯準則選擇期望損失最小的行動集執(zhí)行?；谌Q策的雙激素人工胰臟模型預(yù)測控制算法中，3種決策行動分別對應(yīng)著執(zhí)行胰島素子系統(tǒng)、執(zhí)行零輸注子系統(tǒng)和執(zhí)行胰高血糖素子系統(tǒng)，兩個狀態(tài)集為血糖預(yù)測值處于區(qū)間[180,+∞]和血糖預(yù)測值處于區(qū)間[0,180]。所設(shè)計的三支決策代價矩陣如表1所示，該矩陣由上述狀態(tài)集、行動集組成，矩陣值表示在該狀態(tài)下采用該行動集所付出的代價。

表1 雙激素人工胰臟三支決策切換的代價矩陣Tab.1 The cost matrix of the three-branch decision switching of double-hormone artificial pancreas

代價矩陣中，價格因素由激素的價格和輸注速率2部分構(gòu)成，高低血糖風(fēng)險由預(yù)測血糖與區(qū)間上下限間的偏差組成。假定在血糖預(yù)測值處于區(qū)間[180,+∞]時執(zhí)行3個行動集的代價分別為λ11、λ12、λ13，在血糖預(yù)測值處于區(qū)間[0,180]時執(zhí)行3個行動集的代價分別為λ21、λ22、λ23，有：

(6)

式中：y1、y2、y3分別為3個子系統(tǒng)的血糖預(yù)測值向量；y1i、y2i、y3i分別表示3個向量的分量；yH、yL分別表示由區(qū)間上限和區(qū)間下限組成的向量；m1、m2分別為胰島素和胰高血糖素的價格；R2、R4分別為胰島素和胰高血糖素經(jīng)濟代價的權(quán)重系數(shù)；Q2、Q4分別為高血糖風(fēng)險和低血糖風(fēng)險的權(quán)重系數(shù)；u1代表胰島素輸注速度向量；u3代表胰高血糖素輸注速度向量。

假設(shè)血糖預(yù)測值處于[180,+∞]和[0,180]區(qū)間的概率分別為p1、p2，執(zhí)行3個行動集的期望損失分別為S1、S2、S3，則有如下關(guān)系：

S1=p1λ11+p2λ21

S2=p1λ12+p2λ22

S3=p1λ13+p2λ23

(7)

由于只有2個狀態(tài)集，p1+p2=1，式(7)可轉(zhuǎn)化為：

S1=λ21+p1(λ11-λ21)

S2=λ22+p1(λ12-λ22)

S3=λ23+p1(λ13-λ23)

(8)

根據(jù)貝葉斯準則可知，當S1≤S2且S1≤S3時，執(zhí)行胰島素子系統(tǒng)；當S2≤S1且S2≤S3時，執(zhí)行零輸入子系統(tǒng)；當S3≤S1且S3≤S2時，執(zhí)行胰高血糖素子系統(tǒng)。結(jié)合式(8)可知，選擇執(zhí)行胰島素子系統(tǒng)等價于式(9)成立。

λ21+p1(λ11-λ21)≤λ22+p1(λ12-λ22)

λ21+p1(λ11-λ21)≤λ23+p1(λ13-λ23)

(9)

血糖預(yù)測值處于[180,+∞]時執(zhí)行胰島素子系統(tǒng)是正確決策，血糖預(yù)測值處于[0,180]時執(zhí)行胰島素子系統(tǒng)是錯誤決策。實際進行決策時正確決策的代價應(yīng)小于延遲決策的代價，延遲決策的代價應(yīng)小于錯誤決策的代價，因此有λ11≤λ12≤λ13和λ23≤λ21且λ22≤λ21，則式(9)可轉(zhuǎn)化為：

(10)

綜上可知，選擇執(zhí)行胰島素子系統(tǒng)等價于式(10)成立。同理可得，選擇執(zhí)行零輸注子系統(tǒng)等價于式(11)成立，選擇執(zhí)行胰高血糖素子系統(tǒng)等價于式(12)成立。

(11)

(12)

令α=(λ22-λ21)/(λ11-λ21-λ12+λ22)，β=(λ22-λ23)/(λ22-λ23-λ12+λ13)可知，三支決策理論選擇執(zhí)行期望損失最小的行動集的規(guī)則等價于：當p1≥α?xí)r執(zhí)行胰島素子系統(tǒng)；當β≤p1≤α?xí)r執(zhí)行零輸注子系統(tǒng)；當p1≤β時，執(zhí)行胰高血糖素子系統(tǒng)[14]。

設(shè)計三支決策切換規(guī)則有2種方法，一是根據(jù)設(shè)計的代價矩陣直接計算選擇執(zhí)行3個子系統(tǒng)的決策風(fēng)險值，切換到風(fēng)險最小的子系統(tǒng)；二是根據(jù)設(shè)計的代價矩陣值計算閾值α、β，并與血糖預(yù)測值處于區(qū)間[180,+∞]的概率進行比較，決策切換到相應(yīng)的子系統(tǒng)。

2.4 粒子群算法求解目標函數(shù)

本文算法中通過優(yōu)化胰島素子系統(tǒng)和胰高血糖素子系統(tǒng)的目標函數(shù)得到當前最優(yōu)的控制輸入即激素輸注速率，采用粒子群優(yōu)化算法求解目標函數(shù)[15,16]包括如下步驟：

步驟1：對如式(4)、式(5)所示的目標函數(shù)的參數(shù)、初值和粒子速度進行初始化設(shè)置。

步驟2：根據(jù)初始化值計算設(shè)計的目標函數(shù)值。

步驟3：計算個體極值和群體極值，個體極值是每個潛在解移動過程中使目標函數(shù)值最小時的輸注速率值，群體極值是所有潛在解所搜索到的使目標函數(shù)值最小時的輸注速率值。

步驟4：利用式(13)更新個體速度和位置。

(13)

步驟5：更新目標函數(shù)值，根據(jù)更新的粒子速度和位置計算相應(yīng)的目標函數(shù)值。

步驟6：更新個體極值和群體極值，根據(jù)更新的目標函數(shù)值，從中選出最小值對應(yīng)的激素輸注速率即為個體極值和群體極值。

步驟7：將迭代次數(shù)設(shè)為終止條件，若已達到所設(shè)置的最大迭代次數(shù)，輸出最優(yōu)解；否則，轉(zhuǎn)第4步，尋找下一個最優(yōu)解。

3 仿真結(jié)果分析

為驗證本文所提算法的可行性及有效性，采用美國弗吉尼亞大學(xué)和意大利帕多瓦大學(xué)聯(lián)合開發(fā)的UVA仿真平臺進行仿真實驗[17]。該平臺提供了33個虛擬病人，其中11個成人、11個青少年和11個兒童，基本涵蓋了各個年齡段的糖尿病患者。選擇33個虛擬病人，設(shè)計分別在早上7:00，中午12:00和下午18:00攝入45 g、70 g和80 g的碳水化合物的場景，取預(yù)測步數(shù)為24，控制步數(shù)為10進行為期一天的仿真實驗。實驗中設(shè)置初始血糖值為100 mg/dL，設(shè)置體現(xiàn)個體差異性的胰島素基礎(chǔ)輸注速率，控制目標是將血糖調(diào)節(jié)到[70,180]的區(qū)間內(nèi)，胰島素價格設(shè)為0.3元/U，胰高血糖素價格為 1元/U[8]。

對于仿真結(jié)果分別從個體和整體2部分進行分析，從33個實驗對象中選擇1個個體對其結(jié)果進行分析，該個體的血糖曲線、對應(yīng)的胰島素輸注速率曲線和胰高血糖素輸注速率曲線如圖2所示。其中，圖2(a)是該個體在控制過程中血糖濃度隨時間變化的曲線，圖2(b)是控制過程中胰島素輸注速率隨時間變化的曲線，圖2(c)是控制過程中胰高血糖素輸注速率隨時間變化的曲線。對3個子圖分析得到3個子系統(tǒng)的工作時間分布如圖3所示，其中采用胰島素子系統(tǒng)輸出值為1，采用零輸注子系統(tǒng)輸出值為2，采用胰高血糖素子系統(tǒng)輸出值為3。

圖2 基于三支決策切換規(guī)則的人工胰臟個體控制效果圖Fig.2 Artificial pancreas individual control effect based on three decision switching rules

圖3 3個子系統(tǒng)工作時間分布圖Fig.3 Working time distribution of the three subsystems

實驗過程中設(shè)置該個體的胰島素基礎(chǔ)輸注速率為55.475 6 pmol/min，從圖2(a)可以看出，該個體的血糖值始終處于[70,180]的正常范圍內(nèi)，在設(shè)定值130 mg/dL附近上下浮動；虛擬病人分別在420 min、720 min和1 080 min時用餐，并且血糖隨之升高。綜合分析圖2可知，在0 min時虛擬病人經(jīng)過長時間的睡眠，血糖水平下降，因此輸注胰高血糖素以維持血糖，185 min后血糖逐漸升高，停止輸注胰高血糖素。420 min時病人用餐，血糖大幅度升高但可依靠胰島素基礎(chǔ)輸注速率將血糖控制在正常范圍內(nèi)，因此無需輸注大劑量胰島素，胰島素子系統(tǒng)不工作。隨著餐后血糖不斷被消耗利用，497 min時血糖開始下降，693 min時血糖低于設(shè)定值，輸注胰高血糖素減緩血糖下降速率。720 min時病人再次用餐，血糖不斷升高，782 min時血糖高于區(qū)間上限，輸注大劑量胰島素，801 min時血糖恢復(fù)到正常區(qū)間內(nèi)，胰島素子系統(tǒng)停止工作。1 025 min時血糖低于設(shè)定值，輸注胰高血糖素減緩血糖下降。1 080 min時病人再次用餐，血糖升高，1 142 min時輸注大劑量胰島素，1 161 min時血糖恢復(fù)到正常區(qū)間內(nèi)，胰島素子系統(tǒng)停止工作。綜上所述，本文算法將該個體的血糖濃度控制在正常范圍內(nèi)，且沒有低血糖事件的發(fā)生，獲得了較好的控制效果。

從圖3可以看出，在782～801 min和1 142～ 1 161 min時，胰島素子系統(tǒng)工作，進行大劑量胰島素的輸注；在0～185 min、361～448 min、693～765 min和1 025～1 117 min時，胰高血糖素子系統(tǒng)工作，為虛擬病人輸注胰高血糖素；其余時間即為零輸注子系統(tǒng)的工作時間。

對33個虛擬病人進行仿真實驗，實驗對象血糖值隨時間變化的曲線如圖4所示，反映33個虛擬病人整體的血糖平均狀態(tài)及最值的血糖跟蹤圖如圖5所示，血糖密度圖如圖6所示，控制變異性網(wǎng)格分析(CVGA)如圖7所示。

圖4 33個虛擬病人的血糖濃度曲線Fig.4 Blood glucose concentration curves of 33 virtual patients

圖5 血糖跟蹤軌跡曲線Fig.5 Trace curve of blood glucose

圖6 血糖密度圖Fig.6 Blood glucose density

圖7 控制變異性網(wǎng)格分析圖Fig.7 Schematic diagram of Control-Variability Grid Analysis

從圖4可看出，除個別對象出現(xiàn)短時間輕微高血糖現(xiàn)象外，絕大部分虛擬病人的血糖處于正常范圍內(nèi)，得到了良好的控制效果。圖5中綠色曲線表示33位虛擬病人的平均血糖狀態(tài)，仿真過程中始終處于正常血糖范圍內(nèi)；紅色曲線分別表示33位虛擬病人血糖水平的最值情況，可看出僅有最大值存在輕微高血糖現(xiàn)象。圖6血糖密度圖給出了血糖的概率分布，可看出33個虛擬病人全部血糖值中只有4%處于輕微高血糖，96%均處于正常范圍內(nèi)，無低血糖狀況。圖7中橫、縱坐標均為血糖值，處于A區(qū)域為較好結(jié)果，可看出33個虛擬病人中，73%處于A區(qū)域，27%處于B區(qū)域，取得了較好的控制效果?，F(xiàn)有研究的控制變異性網(wǎng)格分析結(jié)果為40%處于A區(qū)域，60%處于B區(qū)域[5,6]或100%處于B區(qū)域[8]，因此本算法在群體控制質(zhì)量上有著顯著提高。此外，該算法運行時間為526 s，具有較高的運行效率。綜合分析圖4、圖5、圖6和圖7，本算法在33位虛擬病人身上取得了良好的控制效果。

最后，33位虛擬病人的血糖指數(shù)和血糖標準差數(shù)據(jù)如表2所示，主要包括血糖均值、低血糖時間百分比、正常血糖時間百分比、高血糖時間百分比、血糖指數(shù)(BGRI)和血糖標準差數(shù)據(jù)(RoC)。血糖指數(shù)值和血糖標準差越小，風(fēng)險就越小，表格中BGRI最大為4.67，RoC最大為0.93，可看出基于三支決策的雙激素人工胰臟模型預(yù)測控制算法的風(fēng)險很低，具有較高的安全性。

表2 血糖指數(shù)和血糖標準差分析Tab.2 Analysis of glucose index and glucose standard deviation

4 結(jié) 論

針對雙激素人工胰臟子系統(tǒng)切換決策機制研究較少且存在運算效率和群體控制質(zhì)量有待提高的問題，本文提出一種基于三支決策的雙激素人工胰臟模型預(yù)測控制算法。該算法綜合考慮了控制效果和經(jīng)濟代價2個方面因素設(shè)計了三支決策切換規(guī)則，在現(xiàn)有的人工胰臟預(yù)測模型基礎(chǔ)上，設(shè)計了包含設(shè)定值控制和區(qū)間控制兩部分的目標函數(shù)，并用粒子群算法求解。最后對33位虛擬病人進行仿真實驗，從血糖曲線、激素輸注速率曲線、血糖跟蹤圖、血糖密度圖、控制變異性網(wǎng)格分析、血糖指數(shù)和血糖標準差多個方面對控制效果進行分析，證明該算法能很好地將血糖控制在正常范圍內(nèi)，在跟蹤誤差、風(fēng)險指數(shù)和運行效率上取得了良好的效果，在對多個病人的群體控制質(zhì)量上有顯著提高。