CN111341451A - A Model-Free Adaptive Predictive Control Algorithm with Disturbance Compensation for Glycemic Control - Google Patents

Abstract

The invention discloses a model-free adaptive predictive control algorithm with interference compensation for blood sugar control, which comprises the steps of firstly establishing a very local data model with interference compensation, and deducing the structure of a prediction equation according to the model for describing the input-output relationship of a controlled system in a certain period of time in the future; carrying out online estimation on dynamic linearization parameters in a prediction equation by using historical input and output data; calculating a future optimal control input sequence by using a projection algorithm, and solving the control input at the next moment according to a rolling optimization principle; and estimating the interference of the system by using an adaptive extended state observer, and calculating a feeding interference sequence which possibly appears in the future by using a least square method for interference compensation. The control method can effectively deal with the hysteresis of the insulin action process, can inhibit the influence of feeding interference on the blood sugar system of the human body, has the safety and the rapidity of the control effect, and is very suitable for realizing the blood sugar control of the artificial pancreas.

Description

Translated fromChinese

用于血糖控制的具有干扰补偿的无模型自适应预测控制算法A Model-Free Adaptive Predictive Control Algorithm with Disturbance Compensation for Glycemic Control

技术领域technical field

本发明涉及一种人工胰腺控制算法，特别是一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法。The invention relates to an artificial pancreas control algorithm, in particular to a model-free adaptive prediction control algorithm with interference compensation for blood sugar control.

背景技术Background technique

人工胰腺在糖尿病患者人群中具有广阔的应用空间。人工胰腺主要由连续血糖监测装置、胰岛素泵、闭环控制算法三个部分组成，可根据患者的血糖水平调节胰岛素的注射剂量，使患者的血糖水平保持在正常范围内。而闭环控制算法作为人工胰腺的核心，是当今人工胰腺的发展重心，需要应对胰岛素起效的滞后性、人体血糖生理系统的非线性与不确定性等血糖控制存在的难题。Artificial pancreas has a broad application space in the diabetic patient population. The artificial pancreas is mainly composed of three parts: a continuous blood glucose monitoring device, an insulin pump, and a closed-loop control algorithm. It can adjust the insulin injection dose according to the patient's blood sugar level to keep the patient's blood sugar level within the normal range. The closed-loop control algorithm, as the core of artificial pancreas, is the focus of the development of artificial pancreas today.

目前常用的人工胰腺控制算法多为模型预测控制。如文献(Chakrabarty A,Zavitsanou S,Doyle FJ,Dassau E.Event-triggered model predictive control forembedded artificial pancreas systems.IEEE Transactions on BiomedicalEngineering.2018,65:575-586)中采用了事件触发模型预测控制算法。该方法依赖于人体血糖-胰岛素系统模型，所需调节的参数众多，需要较长的训练时间。且该方法的抗干扰性能存在局限性，在面对人体进食带来的较大血糖波动时难以保证较好的控制效果。The most commonly used artificial pancreas control algorithms are model predictive control. For example, the event-triggered model predictive control algorithm was used in the literature (Chakrabarty A, Zavitsanou S, Doyle FJ, Dassau E. Event-triggered model predictive control forembedded artificial pancreas systems. IEEE Transactions on Biomedical Engineering. 2018, 65: 575-586). This method relies on a model of the human blood glucose-insulin system, which requires many parameters to be adjusted and requires a long training time. In addition, the anti-interference performance of this method has limitations, and it is difficult to ensure a better control effect in the face of large blood sugar fluctuations caused by human eating.

发明内容SUMMARY OF THE INVENTION

本发明的目的在于提供一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，将干扰估计和补偿与无模型自适应预测控制算法相结合，从而兼顾控制效果的快速性、安全性与抗干扰性。The purpose of the present invention is to provide a model-free adaptive predictive control algorithm with interference compensation for blood sugar control, which combines interference estimation and compensation with a model-free adaptive predictive control algorithm, thereby taking into account the rapidity and safety of the control effect. performance and anti-interference.

实现本发明目的的技术解决方案为：一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，包括以下步骤：The technical solution for realizing the purpose of the present invention is: a model-free adaptive predictive control algorithm with interference compensation for blood sugar control, comprising the following steps:

步骤1、建立预测方程：首先建立用于预测单步输出的极局部数据模型，再推导出可多步预测输出的预测方程，作为算法的控制框架主体；Step 1. Establish a prediction equation: first, establish a very local data model for predicting single-step output, and then derive a prediction equation that can predict output in multiple steps, as the main body of the algorithm's control framework;

步骤2、动态线性化参数在线估计：根据所得的输入输出数据，估计预测方程中的动态线性化参数；Step 2. Online estimation of dynamic linearization parameters: according to the obtained input and output data, estimate the dynamic linearization parameters in the prediction equation;

步骤3、计算控制输入：在已建立的预测方程的基础上，采用投影算法计算出当前时刻的控制输入；Step 3. Calculate the control input: on the basis of the established prediction equation, use the projection algorithm to calculate the control input at the current moment;

步骤4、基于自适应扩张观测器的干扰估计：利用自适应扩张观测器，估计可能出现的进食干扰对系统带来的扰动；Step 4. Disturbance estimation based on the adaptive expansion observer: use the adaptive expansion observer to estimate the disturbance to the system caused by possible eating disturbances;

步骤5、基于最小二乘法的干扰预测与补偿：使用最小二乘法，对未来可能出现的干扰序列进行预测，并将其用于补偿下一时刻的预测方程中的干扰项。Step 5. Interference prediction and compensation based on the least square method: use the least square method to predict the interference sequence that may appear in the future, and use it to compensate the interference term in the prediction equation at the next moment.

本发明与现有技术相比，其显著优点为：Compared with the prior art, the present invention has the following significant advantages:

(1)本发明的无模型自适应预测控制算法不基于模型，在避免模型参数辨识的前提下保证控制效果依然具有模型预测控制的安全性、有效性；(2)本发明使用的动态线性化参数估计方法可通过输入输出数据在线调整控制器参数，使控制器具有自适应性，不需要手动调节，并且能够应对被控系统性能变化等情况；(3)本发明使用的干扰补偿方法针对可能出现的进食干扰进行估计、预测与补偿，降低进食干扰对控制效果的影响，增强了算法的安全性，这是传统的无模型自适应控制没有的；(4)本发明采用了离散形式，并且采样间隔较大，在保证控制效果的前提下减少计算量，确保了算法的实用性。(1) The model-free adaptive predictive control algorithm of the present invention is not based on a model, and the control effect is guaranteed to still have the safety and effectiveness of the model predictive control under the premise of avoiding model parameter identification; (2) The dynamic linearization used in the present invention The parameter estimation method can adjust the parameters of the controller online through the input and output data, so that the controller has self-adaptation without manual adjustment, and can cope with the performance changes of the controlled system; (3) The interference compensation method used in the present invention aims at possible Estimate, predict and compensate for the eating disturbance that occurs, reduce the influence of eating disturbance on the control effect, and enhance the security of the algorithm, which is not available in the traditional model-free adaptive control; (4) The present invention adopts a discrete form, and The sampling interval is large, which reduces the amount of calculation on the premise of ensuring the control effect, and ensures the practicability of the algorithm.

下面结合附图对本发明作进一步详细描述。The present invention will be described in further detail below with reference to the accompanying drawings.

附图说明Description of drawings

图1为本发明控制算法的总体流程图。FIG. 1 is an overall flow chart of the control algorithm of the present invention.

图2为具有干扰补偿的无模型自适应预测控制算法结构图。FIG. 2 is a structural diagram of a model-free adaptive predictive control algorithm with disturbance compensation.

图3为本发明控制算法的实施步骤流程图。FIG. 3 is a flow chart of the implementation steps of the control algorithm of the present invention.

图4(a)为本发明控制算法的闭环血糖控制系统仿真的血糖水平曲线图。FIG. 4( a ) is a blood sugar level graph simulated by the closed-loop blood sugar control system of the control algorithm of the present invention.

图4(b)为本发明控制算法的闭环血糖控制系统仿真的胰岛素血糖输注速率曲线图。Fig. 4(b) is a graph of the insulin blood glucose infusion rate simulated by the closed-loop blood glucose control system of the control algorithm of the present invention.

具体实施方式Detailed ways

结合图1，本发明的一种用于血糖控制的具有干扰补偿的无模型自适应预测控制算法，包括以下步骤：1, a model-free adaptive predictive control algorithm with interference compensation for blood sugar control of the present invention includes the following steps:

步骤1、建立预测方程。预测方程可以用于描述未来的被控系统输入输出的关系，是本发明控制算法结构的主体，如图2所示。首先设计一种具有干扰补偿的极局部数据模型，用于描述输入变化量与干扰对单步输出变化量的影响；再基于极局部数据模型推导出可用于多步预测的预测方程，用于找到最佳的未来控制输入序列使被控系统输出能够尽可能接近参考轨迹。Step 1. Establish a prediction equation. The prediction equation can be used to describe the relationship between the input and output of the controlled system in the future, and is the main body of the control algorithm structure of the present invention, as shown in FIG. 2 . Firstly, an extremely local data model with disturbance compensation is designed to describe the influence of input change and disturbance on the single-step output change; then, based on the extremely local data model, a prediction equation that can be used for multi-step prediction is derived, which is used to find The optimal future control input sequence enables the controlled system output to be as close as possible to the reference trajectory.

步骤1.1、设计具有干扰补偿的单步极局部数据模型：Step 1.1. Design a single-step pole local data model with interference compensation:

Δy(k+1)＝f(k)+A^T(k)ΔU(k) (1)Δy(k+1)=f(k)+^AT (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)＝[α₁(k),...,α_L(k)]^T为动态线性化参数；L为向量的长度，本发明中根据胰岛素的作用时间选为36，即描述3个小时内胰岛素输注变化量对现在血糖水平变化量的影响。where ΔU(k)=[Δu(k),...,Δu(k-L+1)]^T , Δu(k)=u(k)-u(k-1), u(k) and y(k) are the input and output in the kth sampling interval respectively; f(k) includes the interference and other fast time-varying parts in the system; A(k)=[α₁ (k),...,α_L (k)]^T is the dynamic linearization parameter; L is the length of the vector. In the present invention, 36 is selected according to the action time of insulin, that is, the influence of the change of insulin infusion within 3 hours on the change of the current blood sugar level is described.

步骤1.2、推导多步的极局部数据模型：Step 1.2. Derive a multi-step extremely local data model:

从单步极局部数据模型进行逐步推导，可推导出多步极局部数据模型的形式。为了便于推导，可设：Step-by-step derivation from the single-step pole local data model leads to the form of the multi-step pole local data model. In order to facilitate the derivation, we can set:

利用P与Q，可将过去的输入与未来的控制时域内的输入分离，多步极局部数据模型推导过程如下：Using P and Q, the past input can be separated from the future input in the control time domain. The derivation process of the multi-step pole local data model is as follows:

根据式(3)可总结出多步极局部数据模型的规律：According to formula (3), the rules of the multi-step pole local data model can be summarized:

步骤1.3、推导多步的预测方程：Step 1.3. Derive the multi-step prediction equation:

设计为了便于整理出预测方程，可定义如下向量：Design In order to facilitate sorting out the prediction equation, the following vector can be defined:

Y_P(k+1)＝[y(k+1),...,y(k+H_p)]^T (5)Y_P (k+1)=[y(k+1),...,y(k+H_p )]^T (5)

其中，Y_P(k+1)为预测时域中各采样时刻的系统输出序列，ΔU_P(k)为控制时域中各采样时刻内的控制输入差分序列，F(k)为干扰在预测时域中各采样时刻上对系统输出的影响序列；H_u为ΔU_P(k)的维数，H_p为Y_P(k+1)的维数。由式(4)-(7)可整理出预测方程：Among them, Y_P (k+1) is the system output sequence at each sampling time in the prediction time domain, ΔU_P (k) is the control input difference sequence at each sampling time in the control time domain, and F(k) is the disturbance in the prediction time domain. The sequence of influences on the system output at each sampling time in the time domain; H_u is the dimension of_ΔUP (k), and H_p is the dimension of Y_P (k+1). The prediction equation can be sorted out from equations (4)-(7):

Y_P(k+1)＝Ey(k)+F(k)+A(k)ΔU(k-1)+B(k)ΔU_P(k) (8)Y_P (k+1)=Ey(k)+F(k)+A(k)ΔU(k-1)+B(k)ΔU_P (k) (8)

其中E＝[1,1,...,1]^T；A与B为时变系数矩阵，如下两式所示：Where E=[1,1,...,1]^T ; A and B are time-varying coefficient matrices, as shown in the following two formulas:

步骤2、动态线性化参数在线估计。式(9)(10)中的动态线性化参数A(k+i)，i＝0,...,H_P-1是体现系统动态特性的重要参数，描述了现在与将来的系统输入输出的关系，但是并不能直接获取其数值。因此，这些参数将由其估计值代替，由

i＝0,...,H_P-1表示，并且将根据药代动力学，通过建立代价函数的方式估计。Step 2: Online estimation of dynamic linearization parameters. The dynamic linearization parameter A(k+i) in equations (9) and (10), i=0,...,H_P -1 is an important parameter reflecting the dynamic characteristics of the system, describing the current and future system input and output relationship, but its value cannot be directly obtained. Therefore, these parameters will be replaced by their estimated values, given by

i₌ 0,...,HP-1 represents and will be estimated by means of establishing a cost function according to the pharmacokinetics.

步骤2.1、基于胰岛素作用特性设计动态线性化参数向量Step 2.1. Design dynamic linearization parameter vector based on insulin action characteristics

根据药代动力学，被注射胰岛素在人体内的扩散过程及降血糖作用可以通过一个高阶惯性环节近似，这也表明系统输入变化量与系统输出变化量之间也可用高阶惯性环节来近似。因此动态线性化参数向量A(k)的估计值可以用一组序列进行表示：According to pharmacokinetics, the diffusion process and hypoglycemic effect of injected insulin in the human body can be approximated by a high-order inertial link, which also indicates that the system input change and system output change can also be approximated by a high-order inertial link . Therefore, the estimated value of the dynamic linearization parameter vector A(k) can be represented by a set of sequences:

其中，v₁(k)和v₂(k)为两个待估计的参数。M被定义为:Among them, v₁ (k) and v₂ (k) are two parameters to be estimated. M is defined as:

M＝-[g(T),g(2T),...,g(H_pT)]^T (12)M=-[g(T),g(2T),...,g(H_p T)]^T (12)

其中，g(t)为一个2阶惯性环节脉冲响应函数，T为采样间隔。由于胰岛素对血糖具有促进下降的作用，所以M值取负。式(11)中，H(v₂)是用来调整响应特性的参数矩阵，可表示为：Among them, g(t) is a second-order inertial link impulse response function, and T is the sampling interval. Since insulin has the effect of promoting the decline of blood sugar, the M value is negative. In formula (11), H(v₂ ) is the parameter matrix used to adjust the response characteristics, which can be expressed as:

这种

的表示形式与(Finite Impulse Response)滤波器类似。但与FIR滤波器不同的是，这一参数序列描述的是系统输入差分与输出差分之间的动态特性。使用式(11)的形式，可以减少需要调整的参数个数，并且可以将参数序列

对动态特性的描述能够贴近实际的药代动力学。this

The representation is similar to the (Finite Impulse Response) filter. But unlike an FIR filter, this parameter sequence describes the dynamics between the system's input differential and output differential. Using the form of formula (11), the number of parameters to be adjusted can be reduced, and the parameter sequence can be

The description of the dynamic properties can approximate the actual pharmacokinetics.

步骤2.2、为了求取v₁(k)、v₂(k)的值，建立一个代价函数：Step 2.2. In order to obtain the values of v₁ (k) and v₂ (k), establish a cost function:

等式(14)右边第二项为所设的软约束，防止参数估计对噪声过于敏感，μ为其权重。因此v₁(k)、v₂(k)可表示为：The second term on the right side of equation (14) is the soft constraint set to prevent parameter estimation from being too sensitive to noise, and μ is its weight. Therefore v₁ (k), v₂ (k) can be expressed as:

对于2型糖尿病患者，其血糖生理系统具有慢时变的特性，因此，可以假设在预测时域内A(k)的数值保持不变，在实际应用中对预测方程准确度的影响可忽略不计。因此，

i＝1,...,H_P-1均可被设为与

相等。For patients withtype 2 diabetes, their blood glucose physiological system has slow time-varying characteristics. Therefore, it can be assumed that the value of A(k) remains unchanged in the prediction time domain, and the impact on the accuracy of the prediction equation in practical applications is negligible. therefore,

i=1,...,H_P -1 can all be set to

equal.

步骤3、在式(8)已建立的预测方程的基础上，采用投影算法计算出当前时刻的控制输入。Step 3: On the basis of the established prediction equation in formula (8), the projection algorithm is used to calculate the control input at the current moment.

步骤3.1、为找出未来的最佳控制输入差分序列，基于投影算法理论建立代价函数：Step 3.1. In order to find the optimal control input differential sequence in the future, a cost function is established based on the projection algorithm theory:

其中y^*(k)为k时刻的为参考输出，等式右边的第二项为限制控制输入变化量的软约束，λ_i+1为各时刻控制输入变化量的权重。Where y^* (k) is the reference output at time k, the second term on the right side of the equation is the soft constraint limiting the amount of control input change, and λ_i+1 is the weight of the control input change at each time.

步骤3.2、将式(16)中的代价函数转化为二项式形式：Step 3.2. Convert the cost function in equation (16) into a binomial form:

J＝(Y_P^*(k+1)-Y_P(k+1))^T(Y_P^*(k+1)-Y_P(k+1))+ΔU_P(k)TλΔU_P(k) (17)J=(Y_P^* (k+1)-Y_P (k+1))^T (Y_P^* (k+1)-Y_P (k+1))+ΔU_P (k)TλΔU_P (k) (17)

其中Y_P^*(k+1)＝[y_P^*(k+1),...,y_P^*(k+H_P)]^T，λ＝diag(λ₁,...,λ_Hu)^T。where Y_P^* (k+1)=[y_P^* (k+1),...,y_P^* (k+H_P )]^T , λ=diag(λ₁ ,...,λ_Hu )^T.

步骤3.3、使

并将式(8)代入式(17)，可算出未来最佳控制输入差分序列：Step 3.3, make

Substituting equation (8) into equation (17), the optimal control input differential sequence in the future can be calculated:

步骤3.4、根据滚动优化策略，取序列中的第一个值为控制输入，由下式表示：Step 3.4. According to the rolling optimization strategy, take the first value in the sequence as the control input, which is represented by the following formula:

u(k)＝u(k-1)+d^TΔU_P(k) (19)u(k)=u(k-1)+d^T_ΔUP (k) (19)

其中d＝[1,0,...,0]^T，即u(k)只取ΔU_P(k)中的第一个元素作为控制输入增量。where d=[1,0,...,0]^T , that is, u(k) only takes the first element in_ΔUP (k) as the control input increment.

步骤4、基于自适应扩张观测器的干扰估计。利用自适应扩张观测器，估计下一采样间隔内可能出现的进食干扰对系统带来的扰动，即式(1)中的f(k)。Step 4. Interference estimation based on the adaptive extended observer. Using the adaptive expansion observer, the disturbance to the system caused by the possible feeding disturbance in the next sampling interval is estimated, that is, f(k) in equation (1).

步骤4.1、为估计扰动值f(k)，设计一种离散型的自适应扩张状态观测器，以下为离散型扩张状态观测器的结构：Step 4.1. In order to estimate the disturbance value f(k), a discrete adaptive extended state observer is designed. The following is the structure of the discrete extended state observer:

其中，

为离散型扩张状态观测器的状态量，即为对系统输出y(k)的估计值。

与

为离散型扩张状态观测器的两个扩张状态，其中，取

作为对干扰项f(k)的估计，即使干扰估计值

为提高观测器的跟踪性能，经过实践需取两个扩张状态对

进行观测，形成三阶离散型扩张状态观测器。A_O与B_O为系数矩阵。L_D(k)为自适应增益。in,

is the state quantity of the discrete extended state observer, that is, the estimated value of the system output y(k).

and

are the two extended states of the discrete extended state observer, where, take

As an estimate of the interference term f(k), even if the interference estimate

In order to improve the tracking performance of the observer, it is necessary to take two extended state pairs through practice.

Make observations to form a third-order discrete extended state observer. A_O and B_O are coefficient matrices. L_D (k) is the adaptive gain.

步骤4.2、基于卡尔曼滤波理论，通过下式在每次采样时刻更新自适应增益L_D(k)：Step 4.2. Based on the Kalman filter theory, update the adaptive gain_LD (k) at each sampling moment by the following formula:

其中C_O＝[1,0,0]^T。P_O(k)即对应卡尔曼滤波中的先验估计协方差，其初值可根据初始条件误差来选取。R_O为测量噪声协方差，为已知量。θ为一个加快收敛的系数，通常取：where_CO = [1,0,0]^T . P_O (k) corresponds to the a priori estimated covariance in the Kalman filter, and its initial value can be selected according to the initial condition error. R_O is the measurement noise covariance, which is a known quantity. θ is a coefficient to speed up the convergence, usually taken as:

Q_O为描述由

引入的误差的协方差，可归结为以下形式：Q_O is described by

The covariance of the introduced error can be reduced to the following form:

其中

可根据

在一个采样间隔T内可能出现的数值变动范围选取。in

according to

The range of possible value changes in a sampling interval T is selected.

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

进行观测。Therefore, through the selection of P_O (0), Q_O and R_O at the beginning, the gain of the extended state observer can be adaptively adjusted, and the

make observations.

步骤5、基于最小二乘法的干扰预测与补偿。首先对进食是否出现进行检测，若出现则使用最小二乘法，对未来可能出现的干扰序列进行预测，并将其用于补偿下一时刻的预测方程中的干扰项。Step 5. Interference prediction and compensation based on the least squares method. First, the occurrence of eating is detected. If it occurs, the least squares method is used to predict the interference sequence that may appear in the future, and it is used to compensate the interference term in the prediction equation at the next moment.

步骤5.1、检测进食是否出现，可采用渐消记忆法，建立进食干扰可能性指标：Step 5.1. To detect whether eating occurs, the gradual elimination memory method can be used to establish the possibility of eating interference:

其中I_m(k)＞0为进食干扰可能性指标，为0＜s_m＜1遗忘因子。一旦I_m(k)超过了设定的特定阈值I_th，即认为此刻检测出了进食干扰，并这一时刻设为k_d。假设进食干扰需在一段作用时间内才能被检测出来，由此可根据式(26)中的遗忘因子s_m和阈值的I_th取值设立一个偏移量k_o，并将k_d-k_o时刻设为进食干扰开始的时刻。另外，还需设立一个阈值I_ter＜I_th，用来判断进食干扰持续时段是否结束。若上一时刻I_m大于而这一时刻小于等于I_ter，则可将这一时刻设为进食干扰结束的时刻k_ter。并且可假设，在检测出具有进食干扰的时段内，只有一次进食，并且在该具有进食干扰的阶段中的进食都被视为该阶段开始时的这一次进食。Among them, I_m (k)>0 is an index of the possibility of eating disturbance, and it is a forgetting factor of 0 < s_m <1. Once_Im (k) exceeds a certain set threshold I_th , it is considered that eating disturbance is detected at this moment, and this moment is set as k_d . Assuming that the feeding disturbance can be detected within a period of time, an offset k_o can be established according to the forgetting factor s_m in equation (26) and the value of the threshold I_th , and k_d - k_o The moment is set to the moment when the eating disturbance begins. In addition, a threshold value_Iter <_Ith needs to be established to judge whether the duration of eating disturbance is over. If the last time_Im is greater than and this time is less than or equal to_{Iter , this time can be set as the time k terwhen} the eating disturbance ends. And it can be assumed that there is only one feeding during the period in which the feeding disturbance is detected, and that the feedings in the period with the feeding disturbance are all regarded as the one feeding at the beginning of the period.

步骤5.2、在进食干扰被检测出后，需要对预测方程中的f(k+i)，i＝0,...,H_P-1，即预测时域内的进食干扰进行预测并补偿。设立与血糖生理系统时间常数相近的高阶惯性环节，并对其单位脉冲响应进行采样，由下式表示：Step 5.2. After the eating disturbance is detected, f(k+i),i₌ 0, . Set up a high-order inertial link similar to the time constant of the blood glucose physiological system, and sample its unit impulse response, which is expressed by the following formula:

G_a＝[g_a(1·T),g_a(2·T),...,g_a(jT)]^T_j×1 (27)Ga₌ [_ga (1·T),ga (2·T),_... ,ga₍ jT)]^T_j×1 (27)

其中g_a(t)为脉冲响应函数，j＞H_P。where_ga (t) is the impulse response function, j>H_P .

步骤5.3、先定义一个含进食量大小信息的参数ω^*(k)，并假设ω^*(k)在没有被检测出有进食干扰时为0，由下式表示：Step 5.3. First define a parameter ω^* (k) containing food intake size information, and assume that ω^* (k) is 0 when no eating interference is detected, which is represented by the following formula:

其中

为在检测出进食干扰后被估计的参数。当进食干扰在k_d时刻被检测出，可利用最小二乘法来对

进行在线估计，由式(29)表示：in

are parameters estimated after detection of feeding disturbances. When eating disturbance is detected at time k_d , the least squares method can be used to

Perform online estimation, which is represented by equation (29):

其中

为自适应扩张状态观测器估计出的历史干扰值序列，其维数为k-k_d+k_o+1，并且将随着采样时刻的增加不断增大直至进食干扰时段结束。且

为k-k_d+k_o+1维的单位脉冲响应差分序列，其中Δg_a(kT)＝g_a(kT)-g_a((k-1)T)。in

The sequence of historical disturbance values estimated for the adaptive expansion state observer has a dimension of kk_d +k_o +1, and will continue to increase as the sampling time increases until the end of the feeding disturbance period. and

is a unit impulse response difference sequence of dimension kk_d +_ko +1, where_Δga (kT)₌ ga (kT)−ga ((_k −1)T).

步骤5.4、预测方程中的中的向量[f(k),...,,f(k+H_p-1)]^T能够被其预测值所代替，由下式表示：Step 5.4. The vector [f(k),...,,f(k+H_p -1)]^T in the prediction equation can be replaced by its predicted value, which is represented by the following formula:

这种干扰补偿方法能够在进食干扰出现时，随着采样时刻的增加获取更多的数据，使估计结果更加精确，并且对进食干扰的在线估计与预测能够不断更新预测方程中的干扰补偿，从而能够对控制输入的控制结果滚动优化。This disturbance compensation method can obtain more data with the increase of sampling time when eating disturbance occurs, so as to make the estimation result more accurate, and the online estimation and prediction of eating disturbance can continuously update the disturbance compensation in the prediction equation, so that Able to scroll optimization of control result of control input.

图3是本发明算法的实施步骤流程图。为了验证本发明控制算法的有效性，利用MATLAB/Simulink对闭环血糖控制系统进行数值仿真。给定的血糖水平参考值为90mg/dL，并且在第300分钟给定一次含75g碳水化合物的进食。安全血糖水平区间为70mg/dL至180mg/dL。Fig. 3 is a flow chart of the implementation steps of the algorithm of the present invention. In order to verify the effectiveness of the control algorithm of the present invention, the closed-loop blood glucose control system is numerically simulated by using MATLAB/Simulink. A blood glucose level reference value of 90 mg/dL was given, and a 75 g carbohydrate meal was given at the 300th minute. Safe blood glucose levels range from 70 mg/dL to 180 mg/dL.

图4(a)为本发明控制算法的闭环血糖控制系统仿真的血糖水平曲线；图4(b)为本发明控制算法的闭环血糖控制系统仿真的胰岛素血糖输注速率曲线。由仿真结果可知本发明的控制算法能够较快将血糖水平控制至预定参考值，且不会超出安全血糖水平区间，兼具快速性与安全性。Figure 4(a) is the blood glucose level curve simulated by the closed-loop blood glucose control system of the control algorithm of the present invention; Figure 4(b) is the insulin blood glucose infusion rate curve simulated by the closed-loop blood glucose control system of the control algorithm of the present invention. It can be seen from the simulation results that the control algorithm of the present invention can control the blood sugar level to a predetermined reference value quickly, and will not exceed the safe blood sugar level range, and has both rapidity and safety.

Claims (6)

1. A model-free adaptive predictive control algorithm with interference compensation for glycemic control, comprising the steps of:

step 1, establishing a prediction equation: firstly, establishing a very local data model for predicting single-step output, and then deriving a prediction equation capable of realizing multi-step prediction output as a control frame main body of an algorithm;

step 2, dynamic linearization parameter online estimation: estimating dynamic linearization parameters in a prediction equation according to the obtained input and output data;

step 3, calculating control input: on the basis of the established prediction equation, calculating the control input of the current moment by adopting a projection algorithm;

and 4, estimating the interference based on the adaptive extended observer: estimating the disturbance of feeding interference possibly occurring to the system by using an adaptive extended observer;

and 5, interference prediction and compensation based on a least square method: the possible future occurrence of the interference sequence is predicted using the least squares method and used to compensate for the interference term in the prediction equation at the next time instant.

2. The model-free adaptive predictive control algorithm with interference compensation for glycemic control of claim 1, wherein step 1 builds a predictive equation, specifically:

step 1.1, designing a single-step local data model with interference compensation:

Δy(k+1)＝f(k)+A^T(k)ΔU(k) (1)

wherein Δ u (k) ═ Δ u (k),.. DELTA.u (k-L +1)]^TU (k) -u (k-1), u (k) and y (k) are input and output in the k-th sampling interval, respectively, f (k) includes interference and other fast time-varying parts in the system, a (k) α₁(k),...,α_L(k)]^TIs a dynamic linearization parameter; l is the length of the vector;

step 1.2, deducing a multistep extremely local data model:

gradually deducing from the single-step local data model, deducing the form of the multi-step local data model, and setting:

and separating the past input from the input in the future control time domain by utilizing P and Q, and summarizing the rule of the multi-step extremely local data model:

step 1.3, deriving a prediction equation for multiple steps

The following vectors are defined:

Y_P(k+1)＝[y(k+1),...,y(k+H_p)]^T(5)

wherein, Y_P(k +1) is the system output sequence, Δ U, at each sampling instant in the prediction time domain_P(k) Inputting a differential sequence for controlling the control in each sampling time in the control time domain, wherein F (k) is an influence sequence of interference on system output at each sampling time in the prediction time domain; h_uIs Delta U_P(k) Dimension of (C), H_pIs Y_PThe dimension of (k + 1); the prediction equation can be compiled from equations (4) to (7):

Y_P(k+1)＝E_y(k)+F(k)+A(k)ΔU(k-1)+B(k)ΔU_P(k) (8)

wherein E ═ 1, 1.., 1]^T(ii) a A and B are time-varying coefficient matrixes, and are shown as the following two formulas:

3. the model-free adaptive predictive control algorithm with interference compensation for glycemic control of claim 2, wherein the step 2 estimates the dynamic linearization parameters in the predictive equation based on the established predictive equation and the obtained input and output data, and comprises the following steps:

step 2.1, designing dynamic linearization parameter vectors based on insulin action characteristics

The estimated value of the dynamic linearization parameter vector a (k) can be represented by a set of sequences:

wherein v is₁(k) And v₂(k) Two parameters to be estimated; m is defined as:

M＝-[g(T),g(2T),...,g(H_pT)]^T(12)

wherein g (T) is a 2-order inertia element impulse response function, and T is a sampling interval; since insulin has a function of promoting the decrease of blood sugar, the value of M is negative; in formula (11), H (v)₂) Is a parameter matrix used to adjust the response characteristics, which can be expressed as:

step 2.2, in order to find v₁(k)、v₂(k) Establishes a cost function:

the second term on the right of equation (14) is the set soft constraint to prevent the parameter estimation from being too sensitive to noise, with μ being its weight;

thus v₁(k)、v₂(k) Can be expressed as:

for the type 2 diabetes patients, the blood sugar physiological system has the characteristic of slow time variation, so that the influence on the accuracy of the prediction equation in practical application is negligible on the assumption that the value of A (k) in the prediction time domain is kept unchanged; therefore, the temperature of the molten metal is controlled,

can be set as

Are equal.

4. The model-free adaptive predictive control algorithm with interference compensation for glycemic control of claim 3, wherein the step 3 calculates the control input at the current time by using a projection algorithm based on the established prediction equation, specifically:

3.1, in order to find out the future optimal control input difference sequence, establishing a cost function based on a projection algorithm theory:

wherein y is^*(k) For the reference output at time k, the second term on the right of the equation is a soft constraint that limits the amount of change in the control input, λ_i+1Controlling the weight of the input variable quantity for each moment;

and 3.2, converting the cost function in the formula (16) into a binomial form:

J＝(Y_P^*(k+1)-Y_P(k+1))^T(Y_P^*(k+1)-Y_P(k+1))+ΔU_P(k)^TλΔU_P(k) (17)

wherein Y is_P^*(k+1)＝[y_P^*(k+1),...,y_P^*(k+H_P)]^T，

Step 3.3, make

And substituting equation (8) for equation (17) to calculate the future optimal control input differential sequence:

and 3.4, taking the first value in the sequence as a control input according to a rolling optimization strategy, wherein the first value is represented by the following formula:

u(k)＝u(k-1)+d^TΔU_P(k) (19)

wherein d ═ 1,0,. 0, 0]^TI.e. U (k) takes only Δ U_P(k) As the control input increment.

5. The model-free adaptive predictive control algorithm with disturbance compensation for glycemic control of claim 4, wherein the step 4 estimates the disturbance of the system caused by the feeding disturbance that may occur by using an adaptive extended observer, and specifically comprises:

step 4.1, in order to estimate the disturbance value f (k), a discrete type self-adaptive extended state observer is designed, and the structure of the discrete type extended state observer is as follows:

wherein,

the state quantity of the discrete extended state observer is an estimated value of the output y (k) of the system;

and

is two extended states of a discrete extended state observer, wherein

As an estimate of the interference term f (k), i.e. the interference estimate

To improve the tracking performance of the observer, two extended state pairs are taken

Observing to form a three-order discrete extended state observer; a. the_OAnd B_OIs a coefficient matrix; l is_D(k) Is an adaptive gain;

step 4.2, based on the Kalman filtering theory, updating the adaptive gain L at each sampling moment by the following formula_D(k)：

Wherein C is_O＝[1,0,0]^T；P_O(k) Corresponding to the prior estimation covariance in Kalman filtering, the initial value can be selected according to the initial condition error; r_OTo measure the noise covariance, as a known quantity; θ is a coefficient for accelerating convergence, and is taken as:

Q_Oto describe by

Covariance of introduced error:

wherein

Can be based on

The range of values that may occur within a sampling interval T is selected.

6. The model-free adaptive predictive control algorithm with interference compensation for glycemic control of claim 5, wherein the step 5 uses a least squares method to predict the interference sequences that may occur in the future and to use them to compensate the interference terms in the prediction equation at the next time, specifically:

step 5.1, detecting whether food intake occurs, and establishing a food intake interference possibility index by adopting an evanescent memory method:

wherein I_m(k) 0 is the index of the possibility of food disturbance and is 0 < s_m< 1 forgetting factor; once I has been_m(k) Exceeds a set specific threshold value I_thThat is, it is considered that the eating disturbance is detected at this time, and this time is k_d(ii) a Assuming that eating disturbances are detected only during an action time, the forgetting factor s in equation (26) can be used_mAnd I of threshold value_thSetting an offset k_oAnd will k_d-k_oSetting the time as the time of starting the eating disturbance; in addition, a threshold I is also established_ter＜I_thFor judging whether the eating disturbance duration period is over; if it was at the previous moment I_mGreater than and equal to or less than I at the moment_terThis time can be set as the end of the eating disturbance at time k_ter(ii) a And assuming that there is only one meal during the period when eating disturbances are detected, and that meals in the period with eating disturbances are all considered this meal at the beginning of the period;

step 5.2, after eating disturbance is detected, f (k + i), i ═ 0, H in the prediction equation is needed_P-1, predicting and compensating for eating disturbances in the prediction time domain; establishing a high-order inertia link, and sampling a unit impulse response of the high-order inertia link, wherein the unit impulse response is represented by the following formula:

G_a＝[g_a(1·T),g_a(2·T),...,g_a(jT)]^T_j×1(27)

wherein g is_a(t) is the impulse response function, j > H_P；

Step 5.3, defining a parameter omega containing the food intake size information^*(k) And assume ω^*(k) 0 when no eating disturbance is detected, represented by the following formula:

wherein

Is a parameter that is estimated after the detection of eating disturbances; when the food intake is disturbed at k_dThe time is detected and is paired by the least square method

An online estimation is performed, represented by equation (29):

wherein

A sequence of historical disturbance values estimated for an adaptive extended state observer having dimensions k-k_d+k_o+1 and will increase with increasing sampling instant until the end of the eating disturbance period; and is

Is k-k_d+k_oA unit impulse response differential sequence of +1 dimension, where Δ g_a(kT)＝g_a(kT)-g_a((k-1)T)；

Step 5.4, predict the vector [ f (k), ] in the equation, f (k + H)_p-1)]^TCan be replaced by its predicted value, represented by the following equation:

。

Images