用于血糖控制的具有干扰补偿的无模型自适应预测控制算法的制作方法

本发明涉及一种人工胰腺控制算法，特别是一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法。

背景技术：

人工胰腺在糖尿病患者人群中具有广阔的应用空间。人工胰腺主要由连续血糖监测装置、胰岛素泵、闭环控制算法三个部分组成，可根据患者的血糖水平调节胰岛素的注射剂量，使患者的血糖水平保持在正常范围内。而闭环控制算法作为人工胰腺的核心，是当今人工胰腺的发展重心，需要应对胰岛素起效的滞后性、人体血糖生理系统的非线性与不确定性等血糖控制存在的难题。

目前常用的人工胰腺控制算法多为模型预测控制。如文献(chakrabartya,zavitsanous,doylefj,dassaue.event-triggeredmodelpredictivecontrolforembeddedartificialpancreassystems.ieeetransactionsonbiomedicalengineering.2018,65:575-586)中采用了事件触发模型预测控制算法。该方法依赖于人体血糖-胰岛素系统模型，所需调节的参数众多，需要较长的训练时间。且该方法的抗干扰性能存在局限性，在面对人体进食带来的较大血糖波动时难以保证较好的控制效果。

技术实现要素：

本发明的目的在于提供一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，将干扰估计和补偿与无模型自适应预测控制算法相结合，从而兼顾控制效果的快速性、安全性与抗干扰性。

实现本发明目的的技术解决方案为：一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，包括以下步骤：

步骤1、建立预测方程：首先建立用于预测单步输出的极局部数据模型，再推导出可多步预测输出的预测方程，作为算法的控制框架主体；

步骤2、动态线性化参数在线估计：根据所得的输入输出数据，估计预测方程中的动态线性化参数；

步骤3、计算控制输入：在已建立的预测方程的基础上，采用投影算法计算出当前时刻的控制输入；

步骤4、基于自适应扩张观测器的干扰估计：利用自适应扩张观测器，估计可能出现的进食干扰对系统带来的扰动；

步骤5、基于最小二乘法的干扰预测与补偿：使用最小二乘法，对未来可能出现的干扰序列进行预测，并将其用于补偿下一时刻的预测方程中的干扰项。

本发明与现有技术相比，其显著优点为：

(1)本发明的无模型自适应预测控制算法不基于模型，在避免模型参数辨识的前提下保证控制效果依然具有模型预测控制的安全性、有效性；(2)本发明使用的动态线性化参数估计方法可通过输入输出数据在线调整控制器参数，使控制器具有自适应性，不需要手动调节，并且能够应对被控系统性能变化等情况；(3)本发明使用的干扰补偿方法针对可能出现的进食干扰进行估计、预测与补偿，降低进食干扰对控制效果的影响，增强了算法的安全性，这是传统的无模型自适应控制没有的；(4)本发明采用了离散形式，并且采样间隔较大，在保证控制效果的前提下减少计算量，确保了算法的实用性。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明控制算法的总体流程图。

图2为具有干扰补偿的无模型自适应预测控制算法结构图。

图3为本发明控制算法的实施步骤流程图。

图4(a)为本发明控制算法的闭环血糖控制系统仿真的血糖水平曲线图。

图4(b)为本发明控制算法的闭环血糖控制系统仿真的胰岛素血糖输注速率曲线图。

具体实施方式

结合图1，本发明的一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，包括以下步骤：

步骤1、建立预测方程。预测方程可以用于描述未来的被控系统输入输出的关系，是本发明控制算法结构的主体，如图2所示。首先设计一种具有干扰补偿的极局部数据模型，用于描述输入变化量与干扰对单步输出变化量的影响；再基于极局部数据模型推导出可用于多步预测的预测方程，用于找到最佳的未来控制输入序列使被控系统输出能够尽可能接近参考轨迹。

步骤1.1、设计具有干扰补偿的单步极局部数据模型：

δy(k+1)＝f(k)+a^t(k)δu(k)(1)

其中，δu(k)＝[δu(k),...,δu(k-l+1)]^t，δu(k)＝u(k)-u(k-1)，u(k)和y(k)分别为第k个采样间隔内的输入输出；f(k)包含系统中的干扰以及其他快时变部分；a(k)＝[α1(k),...,αl(k)]^t为动态线性化参数；l为向量的长度，本发明中根据胰岛素的作用时间选为36，即描述3个小时内胰岛素输注变化量对现在血糖水平变化量的影响。

步骤1.2、推导多步的极局部数据模型：

从单步极局部数据模型进行逐步推导，可推导出多步极局部数据模型的形式。为了便于推导，可设：

利用p与q，可将过去的输入与未来的控制时域内的输入分离，多步极局部数据模型推导过程如下：

根据式(3)可总结出多步极局部数据模型的规律：

步骤1.3、推导多步的预测方程：

设计为了便于整理出预测方程，可定义如下向量：

yp(k+1)＝[y(k+1),...,y(k+hp)]^t(5)

其中，yp(k+1)为预测时域中各采样时刻的系统输出序列，δup(k)为控制时域中各采样时刻内的控制输入差分序列，f(k)为干扰在预测时域中各采样时刻上对系统输出的影响序列；hu为δup(k)的维数，hp为yp(k+1)的维数。由式(4)-(7)可整理出预测方程：

yp(k+1)＝ey(k)+f(k)+a(k)δu(k-1)+b(k)δup(k)(8)

其中e＝[1,1,...,1]^t；a与b为时变系数矩阵，如下两式所示：

步骤2、动态线性化参数在线估计。式(9)(10)中的动态线性化参数a(k+i)，i＝0,...,hp-1是体现系统动态特性的重要参数，描述了现在与将来的系统输入输出的关系，但是并不能直接获取其数值。因此，这些参数将由其估计值代替，由i＝0,...,hp-1表示，并且将根据药代动力学，通过建立代价函数的方式估计。

步骤2.1、基于胰岛素作用特性设计动态线性化参数向量

根据药代动力学，被注射胰岛素在人体内的扩散过程及降血糖作用可以通过一个高阶惯性环节近似，这也表明系统输入变化量与系统输出变化量之间也可用高阶惯性环节来近似。因此动态线性化参数向量a(k)的估计值可以用一组序列进行表示：

其中，v1(k)和v2(k)为两个待估计的参数。m被定义为:

m＝-[g(t),g(2t),...,g(hpt)]^t(12)

其中，g(t)为一个2阶惯性环节脉冲响应函数，t为采样间隔。由于胰岛素对血糖具有促进下降的作用，所以m值取负。式(11)中，h(v2)是用来调整响应特性的参数矩阵，可表示为：

这种的表示形式与(finiteimpulseresponse)滤波器类似。但与fir滤波器不同的是，这一参数序列描述的是系统输入差分与输出差分之间的动态特性。使用式(11)的形式，可以减少需要调整的参数个数，并且可以将参数序列对动态特性的描述能够贴近实际的药代动力学。

步骤2.2、为了求取v1(k)、v2(k)的值，建立一个代价函数：

等式(14)右边第二项为所设的软约束，防止参数估计对噪声过于敏感，μ为其权重。因此v1(k)、v2(k)可表示为：

对于2型糖尿病患者，其血糖生理系统具有慢时变的特性，因此，可以假设在预测时域内a(k)的数值保持不变，在实际应用中对预测方程准确度的影响可忽略不计。因此，i＝1,...,hp-1均可被设为与相等。

步骤3、在式(8)已建立的预测方程的基础上，采用投影算法计算出当前时刻的控制输入。

步骤3.1、为找出未来的最佳控制输入差分序列，基于投影算法理论建立代价函数：

其中y^*(k)为k时刻的为参考输出，等式右边的第二项为限制控制输入变化量的软约束，λi+1为各时刻控制输入变化量的权重。

步骤3.2、将式(16)中的代价函数转化为二项式形式：

j＝(yp^*(k+1)-yp(k+1))^t(yp^*(k+1)-yp(k+1))+δup(k)tλδup(k)(17)

其中yp^*(k+1)＝[yp^*(k+1),...,yp^*(k+hp)]^t，λ＝diag(λ1,...,λhu)^t。

步骤3.3、使并将式(8)代入式(17)，可算出未来最佳控制输入差分序列：

步骤3.4、根据滚动优化策略，取序列中的第一个值为控制输入，由下式表示：

u(k)＝u(k-1)+d^tδup(k)(19)

其中d＝[1,0,...,0]^t，即u(k)只取δup(k)中的第一个元素作为控制输入增量。

步骤4、基于自适应扩张观测器的干扰估计。利用自适应扩张观测器，估计下一采样间隔内可能出现的进食干扰对系统带来的扰动，即式(1)中的f(k)。

步骤4.1、为估计扰动值f(k)，设计一种离散型的自适应扩张状态观测器，以下为离散型扩张状态观测器的结构：

其中，为离散型扩张状态观测器的状态量，即为对系统输出y(k)的估计值。与为离散型扩张状态观测器的两个扩张状态，其中，取作为对干扰项f(k)的估计，即使干扰估计值为提高观测器的跟踪性能，经过实践需取两个扩张状态对进行观测，形成三阶离散型扩张状态观测器。ao与bo为系数矩阵。ld(k)为自适应增益。

步骤4.2、基于卡尔曼滤波理论，通过下式在每次采样时刻更新自适应增益ld(k)：

其中co＝[1,0,0]^t。po(k)即对应卡尔曼滤波中的先验估计协方差，其初值可根据初始条件误差来选取。ro为测量噪声协方差，为已知量。θ为一个加快收敛的系数，通常取：

qo为描述由引入的误差的协方差，可归结为以下形式：

其中可根据在一个采样间隔t内可能出现的数值变动范围选取。

因此在初始时通过对po(0)、qo与ro的选取，即能实现自适应整定扩张状态观测器的增益，对进行观测。

步骤5、基于最小二乘法的干扰预测与补偿。首先对进食是否出现进行检测，若出现则使用最小二乘法，对未来可能出现的干扰序列进行预测，并将其用于补偿下一时刻的预测方程中的干扰项。

步骤5.1、检测进食是否出现，可采用渐消记忆法，建立进食干扰可能性指标：

其中im(k)＞0为进食干扰可能性指标，为0＜sm＜1遗忘因子。一旦im(k)超过了设定的特定阈值ith，即认为此刻检测出了进食干扰，并这一时刻设为kd。假设进食干扰需在一段作用时间内才能被检测出来，由此可根据式(26)中的遗忘因子sm和阈值的ith取值设立一个偏移量ko，并将kd-ko时刻设为进食干扰开始的时刻。另外，还需设立一个阈值iter＜ith，用来判断进食干扰持续时段是否结束。若上一时刻im大于而这一时刻小于等于iter，则可将这一时刻设为进食干扰结束的时刻kter。并且可假设，在检测出具有进食干扰的时段内，只有一次进食，并且在该具有进食干扰的阶段中的进食都被视为该阶段开始时的这一次进食。

步骤5.2、在进食干扰被检测出后，需要对预测方程中的f(k+i)，i＝0,...,hp-1，即预测时域内的进食干扰进行预测并补偿。设立与血糖生理系统时间常数相近的高阶惯性环节，并对其单位脉冲响应进行采样，由下式表示：

ga＝[ga(1·t),ga(2·t),...,ga(jt)]^tj×1(27)

其中ga(t)为脉冲响应函数，j＞hp。

步骤5.3、先定义一个含进食量大小信息的参数ω^*(k)，并假设ω^*(k)在没有被检测出有进食干扰时为0，由下式表示：

其中为在检测出进食干扰后被估计的参数。当进食干扰在kd时刻被检测出，可利用最小二乘法来对进行在线估计，由式(29)表示：

其中为自适应扩张状态观测器估计出的历史干扰值序列，其维数为k-kd+ko+1，并且将随着采样时刻的增加不断增大直至进食干扰时段结束。且为k-kd+ko+1维的单位脉冲响应差分序列，其中δga(kt)＝ga(kt)-ga((k-1)t)。

步骤5.4、预测方程中的中的向量[f(k),...,,f(k+hp-1)]^t能够被其预测值所代替，由下式表示：

这种干扰补偿方法能够在进食干扰出现时，随着采样时刻的增加获取更多的数据，使估计结果更加精确，并且对进食干扰的在线估计与预测能够不断更新预测方程中的干扰补偿，从而能够对控制输入的控制结果滚动优化。

图3是本发明算法的实施步骤流程图。为了验证本发明控制算法的有效性，利用matlab/simulink对闭环血糖控制系统进行数值仿真。给定的血糖水平参考值为90mg/dl，并且在第300分钟给定一次含75g碳水化合物的进食。安全血糖水平区间为70mg/dl至180mg/dl。

图4(a)为本发明控制算法的闭环血糖控制系统仿真的血糖水平曲线；图4(b)为本发明控制算法的闭环血糖控制系统仿真的胰岛素血糖输注速率曲线。由仿真结果可知本发明的控制算法能够较快将血糖水平控制至预定参考值，且不会超出安全血糖水平区间，兼具快速性与安全性。