Active Disturbance Rejection Control for Artificial Pancreas System Based on Insulin Basal Rate Estimation
-
摘要:
胰岛素基础率是人工胰腺系统实现人体血糖闭环控制的基准, 但该变量在临床治疗中难以准确确定. 针对这一问题, 本文设计了一种基于胰岛素基础率动态估计的人工胰腺自抗扰控制方法, 通过扩张状态观测器(Extended state observer, ESO)实时估计血糖代谢过程中的内部与外界干扰, 构建具备参数自适应能力的反馈控制律和胰岛素注射安全约束, 实现血糖闭环调控能力的有效改善. 在此基础上, 本文设计了基于移动设备和蓝牙模块的人工胰腺软件系统, 并通过美国食品药品监督管理局(Food and Drug Administration, FDA)接受的UVA/Padova T1DM仿真平台完成算法的比较仿真与功能测试. 本文的工作将为后续人工胰腺临床试验的开展提供方法基础和技术支持, 也为我国糖尿病患者血糖管理的改善提供精准医学治疗手段.
Abstract:Insulin basal rate provides the reference for closed-loop blood glucose regulation using artificial pancreas systems, but this quantity is usually difficult to determine accurately in clinical practice. In this regard, this paper introduces an active disturbance rejection control method for artificial pancreas systems based on dynamic estimation of the basal rate. To enable improved glucose regulation, an extended state observer (ESO) is employed to estimate the internal and external disturbances in the glucose metabolic process, and a feedback control law and insulin infusion safety constraints that both incorporate parameter adaptation are proposed. Based on the proposed method, an artificial pancreas software system is designed for mobile devices with Bluetooth modules. The proposed results are evaluated through comparative simulations and functionality tests by using the US FDA (Food and Drug Administration)-accepted UVA/Padova T1DM simulator. The obtained results provide methodological and technical support for further clinical studies of artificial pancreas systems, and introduce a precision medicine solution to enhanced glucose management for Chinese patients with diabetes mellitus.
-
数据降维方法在众多领域应用广泛, 其划分依据也不尽相同, 按照数据结构特征保持与否的准则进行划分, 则可根据数据的全局结构保持和局部结构保持分成两类[1], 前者反映了数据的外部形状, 后者反映了数据的内在属性, 可以寻找出高维观测数据中所隐藏的低维流形结构.其中, 主成分分析(Principal component analysis,PCA)[2]、独立元分析[3]和人工神经网络[4]等方法均为数据全局特征结构保持的代表方法, 核主成分分析(Kernel principal component analysis, KPCA)[5]通过非线性映射函数将线性不可分的原始样本数据输入空间通过投影变换到线性可分的高维特征空间, 然后在新的特征空间中利用线性方法完成主成分分析, 从而实现数据整体方差最大化, 但KPCA方法只能够提取数据的全局结构信息, 若数据中低维局部结构中包含较多特征信息的话, 则效果较差.另外, 流形学习能够从高维历史信息中获取数据间有效的内部联系, 从而得以保持局部结构特征, 具有良好的非线性数据内部属性的处理能力[6-7], 代表性的流行学习[8]方法主要包括等距特征映射算法(Isometric feature mapping, ISOMAP)[9], 拉普拉斯特征映射算法(Laplacian eigenmaps, LE)[10], 局部线性嵌入算法(Locally linear embedding, LLE)[11], 局部保持投影算法(Locality preserving projections, LPP)[12]等, 其中, 有学者在LPP算法中引入核方法, 提出核局部保持投影(Kernel locality preserving projection, KLPP)[13], 其在保持局部结构特征的同时实现线性计算, 反映出数据的局部结构特征, 但本质上KLPP是一种基于局部结构保持的降维方法, 它并不能有效提取出数据的全局特征信息[14].
针对以上问题, 本研究拟对KPCA与KLPP相结合的降维方法进行探讨, 提出了本文的解决办法, 并把新提出的算法命名为改进全局与局部结构保持算法(Global and local structure preserving, GLSP), 在进行原始数据的投影变换时, 既考虑全局结构得以保持, 也兼顾保持局部近邻结构.首先使用局部与全部特征提取方法, 解决数据有效降维的问题, 使用聚类分析中类内距离与类间距离等作为衡量指标, 并使用K近邻(K-nearest neighbor, KNN)方法进行故障的检测[15].本文数据使用柴油机仿真故障数据[16]和TE过程公共数据集, 用于验证方法的有效性.
1. 改进结构保持方法
作为非线性特征提取的经典方法, KPCA通过非线性映射将线性不可分的原始数据从低维空间变换到一个线性可分的高维特征空间, 运用线性方法进行数据降维与特征提取, 其目标是使得数据方差最大化, 但数据方差指标主要用来描述数据集的全局结构信息.此外, KLPP是通过建立样本点之间的近邻关系来保持数据集的局部结构, 本质是保持原始数据局部结构和内部属性.综合考虑KPCA与KLPP两种投影保持方法的思想, 本文提出GLSP, 其目标函数可以理解为由全局目标函数和局部目标函数共同组成.
1.1 局部结构保持算法描述
KLPP通过非线性投影映射, 在投影空间建立近邻图, 最大限度地保持了数据集的近邻结构, 其局部结构保持目标函数定义如下[17]:假设数据集$ {\pmb X} = {\left[ {{{\pmb x}_1}, {{\pmb x}_2}, \cdots , {{\pmb x}_n}} \right]^{\rm T}} \in {\textbf{R}^{n \times m}} $, n为样本个数, m为数据维数, 通过非线性映射$ {\pmb \Phi} $将原始数据映射到高维空间中, 记为$ {\pmb \Phi} \left( {{\pmb x}_i} \right) $, $ {{\pmb J}_{\rm{local}}}\left( {\pmb w} \right) $的目标是在特征空间中找到投影向量$ {\pmb w} $, 使得投影$ {{\pmb y}_i} = {\pmb \Phi}^{\rm T} {\left( {{{\pmb x}_i}} \right)}{\pmb w} $在高维特征空间保持数据点之间的近邻关系, 可以认为, 如果$ {\pmb \Phi} \left( {{{\pmb x}_i}} \right) $和$ {\pmb \Phi} \left( {{{\pmb x}_j}} \right) $是近邻, 那么$ {{\pmb y}_i} = {\pmb \Phi}^{\rm T}{\left( {{{\pmb x}_i}} \right)}{\pmb w} $和$ {{\pmb y}_j} = {\pmb \Phi}^{\rm T} {\left( {{{\pmb x}_j}} \right)}{\pmb w} $也是近邻的.其局部结构保持的目标函数定义为:
$$ \begin{align} &{{\pmb J}_{\rm{local}}}\left( {\pmb w} \right) = \mathop {\min }\limits_{\pmb w} \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{\left\| {{{\pmb y}_i} - {{\pmb y}_j}} \right\|}^2}{{\pmb s}_{ij}}} } = \\&\quad \mathop {\min }\limits_{\pmb w} \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{\left\| {{\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_i}} \right)}}{\pmb w} - {\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_j}} \right)}} {\pmb w}} \right\|}^2}{{\pmb s}_{ij}}} } = \\&\quad \mathop {\min }\limits_{\pmb w} \left\{ {{{\pmb w}^{\rm T}}{\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_i}} \right)}}\left( {{\pmb D} - {\pmb S}} \right){\pmb \Phi} \left( {{{\pmb x}_j}} \right){\pmb w}} \right\} \end{align} $$ (1) 式中, $ {{\pmb s}_{ij}} $为权重参数, 表示数据点之间的近邻关系, $ {\pmb S} $为权重矩阵, $ {\pmb D} $为对角矩阵, $ {{\pmb D}_{ii}} = \sum_{j = 1}^n {{{\pmb s}_{ij}}} $, $ {{\pmb s}_{ij}} $取值一般为:
$$ \begin{align*} {\pmb s}_{ij} = \begin{cases} \dfrac{{\exp \left( { - {{\left\| {{{\pmb x}_i} - {{\pmb x}_j}} \right\|}^2}} \right)}}{t}, &{{\pmb x}_i} \in \Omega _{{{\pmb x}_j}}^k, {{\pmb x}_j} \in \Omega _{{{\pmb x}_i}}^k\\ 0, &\mbox{其他} \end{cases} \end{align*} $$ (2) 式中, $ \Omega _{{{\pmb x}_i}}^k $为$ {{\pmb x}_i} $的k邻域.
引入核函数$ {{\pmb K}_{ij}} = {\pmb K}\left( {{{\pmb x}_i}, {{\pmb x}_j}} \right) = {\pmb \Phi}^{\rm T} {\left( {{{\pmb x}_i}} \right)}{\pmb \Phi} \left( {{{\pmb x}_j}} \right) $, 并存在系数$ {\pmb \alpha} = [{\alpha _1}, {\alpha _2}, \cdots , {\alpha _n}] $对特征空间中的样本线性表示为$ {\pmb w} = \sum_{i = 1}^n {{\alpha _i}{\pmb \Phi} \left( {{{\pmb x}_i}} \right)} $, 局部结构保持的目标函数转化为:
$$ \begin{align} &{{\pmb J}_{\rm{local}}}\left( {\pmb \alpha} \right) = \mathop {\min }\limits_{\pmb \alpha} \Big\{ {{\pmb \alpha} ^{\rm T}}{\pmb \Phi} \left( {{{\pmb x}_i}} \right){\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_j}} \right)}}\\&\quad \left( {{\pmb D} - {\pmb S}} \right) {\pmb \Phi} \left( {{{\pmb x}_i}} \right)\boldsymbol {\Phi} ^{\rm T}{{\left( {{{\pmb x}_j}} \right)}}\boldsymbol {\alpha} \Big\} = \\&\quad \mathop {\min }\limits_{\pmb \alpha} \left\{ {{\boldsymbol {\alpha} ^{\rm T}}{\pmb \Phi} \left( {{{\pmb x}_i}} \right){\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_j}} \right)}}{\pmb L}{\pmb \Phi} \left( {{{\pmb x}_i}} \right){\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_j}} \right)}}{\pmb \alpha} } \right\} = \\&\quad \mathop {\min }\limits_{\pmb \alpha} {{\pmb \alpha} ^{\rm T}}{\pmb{ KLK}}{\pmb \alpha} = \mathop {\min }\limits_{\pmb \alpha} {{\pmb \alpha} ^{\rm T}}{\pmb L}'{\pmb \alpha} \end{align} $$ (3) 式中, $ {\pmb L} = {\pmb D} - {\pmb S} $为Laplacian矩阵, $ {\pmb L}' = \pmb{ KLK} $.
KLPP算法的目的是使数据在高维映射空间中, 仍能保持数据之间的近邻结构, 但是算法本身忽略了对数据集的整体结构特征描述.其本质是因为局部目标函数表达式中没有显式地考虑样本点的全局特征, 只是利用局部结构来代替全局信息, 导致了数据在低维映射中全局特征的扭曲显示.
1.2 全局结构保持算法描述
PCA算法常用于数据主要成分的分析与维度约简, 其本质是将线性数据变换为各个维度线性无关表示的几组数据, 便于对数据中主要特征分量的提取. KPCA算法由核映射将数据映射到高维核空间, 然后使用PCA方法, 与KLPP算法类似, 通过非线性映射$ {\pmb \Phi} $将原始数据$ {\pmb X} = {\left[ {{{\pmb x}_1}, {{\pmb x}_2}, \cdots , {{\pmb x}_n}} \right]^{\rm T}} \in {\textbf{R}^{n \times m}} $映射到高维空间, 记为$ {\pmb \Phi}( {\pmb x} _i ) $, 经过投影向量$ {\pmb w} $投影后的映射$ {{\pmb y} _i} = {\pmb \Phi}^{\rm T}{\left( {{{\pmb x}_i}} \right)}{\pmb w} $, 在投影方向上保证数据方差最大化, 这样可以充分利用高阶统计信息和全局特征结构保持.其全局结构目标函数定义为:
$$ \begin{align} {{\pmb J} _{\rm{global}}}\left( {\pmb w} \right) = \, & \mathop {\max }\limits_{\pmb w} \sum\limits_{i = 1}^n {{{\pmb y}_i}^2} = \\& \mathop {\max }\limits_{\pmb w} \sum\limits_{i = 1}^n {{{\left( {{\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_i}} \right)}}{\pmb w}} \right)}^2}} \end{align} $$ (4) 式中, $ {{\pmb w}^{\rm T}}{\pmb w} = 1 $.
引入核函数$ {{\pmb K}_{ij}} = {\pmb K}\left( {{{\pmb x}_i}, {{\pmb x}_j}} \right) = {\pmb \Phi}^{\rm T} {\left( {{{\pmb x}_i}} \right)}{\pmb \Phi} \left( {{{\pmb x}_j}} \right) $, 可看出即使不确定核函数的具体表达形式, 但是其转化为映射后的数据的内积运算, 存在系数$ {\pmb \alpha} = [{\alpha _1}, {\alpha _2}, \cdots , {\alpha _n}] $对特征空间中的样本线性表示为$ {\pmb w} = \sum_{i = 1}^n {{\alpha _i}{\pmb \Phi} \left( {{{\pmb x}_i}} \right)} $, 全局结构保持的目标函数转化为:
$$ \begin{align} {{\pmb J}_{\rm{global}}}\left( {\pmb \alpha} \right) = \, & \mathop {\max }\limits_{\pmb \alpha} \sum\limits_{i = 1}^n {{{\left( {{\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_i}} \right)}}\sum\limits_{j = 1}^n {{\alpha _j}{\pmb \Phi} \left( {{{\pmb x}_j}} \right)} } \right)}^2}} = \\& \mathop {\max }\limits_{\pmb \alpha} \sum\limits_{i = 1}^n {{{\left( {\sum\limits_{j = 1}^n {{{\pmb \alpha}_j} {\pmb \Phi}^{\rm T} {{\left( {{{\pmb x}_i}} \right)}} {\pmb \Phi} \left( {{{\pmb x}_j}} \right)} } \right)}^2}} = \\& \mathop {\max }\limits_{\pmb \alpha} {{\pmb \alpha} ^{\rm T}}\pmb{ KK}{\pmb \alpha} = \mathop {\max }\limits_{\pmb \alpha} {{\pmb \alpha}^{\rm T} }{\pmb C}{\pmb \alpha} \end{align} $$ (5) 式中, $ {{\pmb \alpha} ^{\rm T}}{\pmb K}{\pmb \alpha} = 1 $, $ {\pmb C} = \pmb{ KK} $.
因为KPCA是一种面向全局的变换方法, 由于其保持了原始数据的大部分方差信息, 所以实现了全局结构的特征提取.然而, 保持全局结构的目标函数中没有考虑各类数据点之间的内部联系, 在低维空间里, 数据点之间的局部几何关系与内在属性可能被忽视, 甚至导致重要信息的丢失.
1.3 改进的整体目标函数
全局结构目标保持函数的思想是在最大程度保持全局信息方差不变的情况下, 提取出样本数据的非线性特征; 而局部结构保持目标函数的思想是在投影中保持样本对之间的远近亲疏关系, 在低维空间中最小化近邻样本间的距离加权平方和, 即尽量避免样本集的发散.
结合全局结构保持目标函数和局部结构保持目标函数的意义, 构造一种新的结构保持目标函数, 使得映射后的特征空间能够保留全局结构的同时, 又可以保持数据间的局部近邻结构[18-19], 构造出如式(6)的最大值选择问题, 由于局部结构保持的目标函数为求取最小值, 所以式(6)中引入其相反数:
$$ \begin{equation} {\pmb J}\left( {\pmb \alpha} \right) = \mathop {\max }\limits_{\pmb \alpha} \left\{ {{{\pmb J}_{\rm{global}}}\left( {\pmb \alpha} \right), - {{\pmb J}_{\rm{local}}}\left( {\pmb \alpha} \right)} \right\} \end{equation} $$ (6) 由于式(6)只从两个目标函数中选择其一, 是一个典型的多目标求取最值问题, 由于未考虑两个目标函数的综合效果, 通常很难求解到全局最优解, 因此, 将式(6)中两个目标函数进行求和操作, 求和后的目标函数如式(7)所示:
$$ \begin{align} {\pmb J}\left( {\pmb \alpha} \right) = \, &\mathop {\max }\limits_{\pmb \alpha} \left( {{{\pmb J}_{\rm{global}}} \left( {\pmb \alpha} \right) - {{\pmb J}_{\rm{local}}} \left( {\pmb \alpha} \right)} \right) = \\& \mathop {\max }\limits_{\pmb \alpha} \left( {{{\pmb \alpha} ^{\rm T}} {\pmb C}{\pmb \alpha} - {{\pmb \alpha} ^{\rm T}}{\pmb L}'{\pmb \alpha} } \right) = \\& \mathop {\max }\limits_{\pmb \alpha} \left( {{{\pmb \alpha} ^{\rm T}} \left( {{\pmb C} - {\pmb L}'} \right) {\pmb \alpha} } \right) \end{align} $$ (7) 同理, 这两个目标函数很难同时达到最佳效果.考虑到两个目标函数之间的差异, 所以, 引入一个权重参数$ \beta $来平衡上述两个目标函数, $ \beta $是一个介于0和1之间的值, 其值的大小对于新的目标函数$ {\pmb J}\left( {\pmb \alpha} \right) $有很大影响, 因为它决定着两个原始目标函数的重要性问题[20].实际上$ \beta $可以看作是平衡两个目标函数的能量变化. $ \beta $越小越侧重于全局特征的提取, $ \beta $越大越侧重于局部特征的提取, $ \beta $值的选取按照如下准则:
$$ \begin{equation} \beta {S_{\rm{global}}} = \left( {1 - \beta } \right){S_{\rm{local}}} \end{equation} $$ (8) 其中, $ {S_{\rm{global}}} $和$ {S_{\rm{local}}} $分别表示$ {{\pmb J}_{\rm{global}}}\left( {\pmb w} \right) $和$ {{\pmb J}_{\rm{local}}}\left( {\pmb w} \right) $的规模大小, 受到参考文献[18, 21]启发, 定义为:
$$ \begin{equation} {S_{\rm{global}}} = \rho \left( {\pmb C} \right) \end{equation} $$ (9) $$ \begin{equation} {S_{\rm{local}}} = \rho \left( {{\pmb L}'} \right) \end{equation} $$ (10) 式中, $ \rho \left( \cdot \right) $是相关矩阵谱半径.结果表明, 该平衡参数的引入策略能够很好地平衡全局和局部的行为, GLSP的降维性能也可以得到保证.事实上, 权重参数也可以根据不同背景下原始数据的特性, 赋予不同定义, 而不仅限于本文所提方法, 这使得改进的结构保持算法更加灵活.结合式(8)~(10), 可得到权重参数$ \beta $的计算公式如下:
$$ \begin{equation} {\beta _{\pmb C} } = \frac{{\rho \left( {\pmb C} \right)}}{{\rho \left( {\pmb C} \right) + \rho \left( {{\pmb L}'} \right)}}, \beta _{{\pmb L}'} = 1 - {\beta _{\pmb C}} \end{equation} $$ (11) 因此, 改进的整体目标函数表示为:
$$ \begin{align} {\pmb J}\left( {\pmb \alpha} \right) = \, &\mathop {\max } \limits_{\pmb \alpha} \left( {{{\beta _{\pmb C} }{\pmb J}_{\rm{global}}} \left( {\pmb \alpha} \right) - {{\beta _{{\pmb L}'} }{\pmb J}_{\rm{local}}} \left( {\pmb \alpha} \right)} \right) = \\& \mathop {\max }\limits_{\pmb \alpha} \left( {{{\pmb \alpha} ^{\rm T}} \left( {{\beta _{\pmb C} }{\pmb C} - \beta _{{\pmb L}'}{\pmb L}'} \right){\pmb \alpha} } \right) = \\& \mathop {\max }\limits_{\pmb \alpha} {{\pmb \alpha} ^{\rm T}}{\pmb M}{\pmb \alpha} \end{align} $$ (12) 式中$ {\pmb M} = {\beta _{\pmb C}}{\pmb C} - \beta _{{\pmb L}'}{\pmb L}', 0 \le {\beta _{\pmb C}} \le 1, 0 \le \beta _{{\pmb L}'} \le 1 $.
最后将上述目标函数的优化问题转化为求解特征向量问题.确定权重参数$ \beta $的值之后, 结合式(5)中的条件, 引入拉格朗日乘子法, 求解特征向量:
$$ \begin{equation} {\pmb L} = {{\pmb \alpha} ^{\rm T}}{\pmb M}{\pmb \alpha} - \lambda \left( {{{\pmb \alpha} ^{\rm T}}{\pmb K}{\pmb \alpha} - 1} \right) \end{equation} $$ (13) 当$ \frac{{\partial {\pmb L}}}{{\partial {\pmb \alpha} }} = 0 $, 可得:
$$ \begin{equation} {\pmb M}{\pmb \alpha} = \lambda {\pmb K}{\pmb \alpha} \end{equation} $$ (14) 非线性问题的求解过程中, 引入正则化方法, 我们用$ {\pmb K} + \eta {{\pmb I}_n} $来代替式(14)中的$ {\pmb K} $, 其中, $ \eta $是一个很小的正整数, $ {{\pmb I}_n} $是一个$ n \times n $的单位向量.
与PCA类似, 本文使用累积方差贡献率准则选取满足贡献率达到要求的主成分个数, 依据式(15)选取前$ p $个特征值确定主成分个数, 本研究中贡献率选定为$ 85\, \% $.
$$ \begin{equation} \frac{{\sum\limits_{k = 1}^p {{\lambda _k}} }}{{\sum\limits_{k = 1}^N {{\lambda _k}} }} > 85\, \% \end{equation} $$ (15) 2. 降维方法实现
原始数据在降维时, 将面临难以全面提取有用信息的困难, 为解决这一问题, 本文提供一种结合全局与局部结构保持的数据降维思想, 实现方法为:采集原始数据, 建立能够从多角度反映数据信息的高维数据集, 并加以验证; 再将数据集输入所提全局与局部结构保持算法中进行降维处理:将低维特征子集输入KNN最近邻分类器, 计算KNN的识别率, 并将聚类分析中类间距与类内距的比值$ {{\pmb S}_B}/{{\pmb S}_W} $作为衡量降维效果指标[22].
2.1 降维评价指标
Fisher判别分析是模式识别方法中的一种数据降维与分类方法.其通过投影将测试数据映射到不同方向, 使得不同类别的测试样本的投影的类间离散度最大, 类内离散度最小[23].类内距$ {{\pmb S}_W} $描述同一类样本内部分布的紧密程度, 而类间距$ {{\pmb S}_B} $用反映不同类别之间的分离程度, 定义如下[15, 24]:
$$ \begin{align} {{\pmb S}_B} = \, & \frac{1}{N}\sum\limits_{i = 1}^C {{N_i}\left( {{{\pmb m}_i} - {\pmb m}} \right){{\left( {{{\pmb m}_i} - {\pmb m}} \right)}^{\rm T}}} \end{align} $$ (16) $$ \begin{align} {{\pmb S}_W} = \, & \frac{1}{N}\sum\limits_{i = 1}^C \sum\limits_{j = 1}^{{N_i}} \left[ {{\pmb \Phi} \left( {{\pmb x}_i^j} \right) - {{\pmb m}_i}} \right]\times\\&{{\left[ {{\pmb \Phi} \left( {{\pmb x}_i^j} \right) - {{\pmb m}_i}} \right]}^{\rm T}} \end{align} $$ (17) 其中, $ {{\pmb m}_i} $表示特征空间中第$ i $类采样均值, $ {\pmb m} $表示所有样本点在特征空间中的均值.显然, $ {{\pmb S}_B}/{{\pmb S}_W} $越大说明该方法的分类与聚类效果越好, 因此将该指标作为降维效果的综合衡量指标之一.
KNN是对不同类别的数据信息根据训练样本特征进行分类的方法, 具有操作直观、效果稳定、时效性强等优点, 广泛应用到各类数据分类领域, 尤其是故障数据的诊断与分类中.原始数据进行降维操作, 其最终目的是实现不同故障类别的准确分类, 故KNN方法的识别率越高, 其反映出对数据的初始降维方法越好[25-26].
2.2 降维方法流程
总结局部与全局结构保持算法的流程图如图 1所示.
算法主要流程如下.
步骤1. 对于数据集$ {\pmb X} = {\left[ {{{\pmb x}_1}, {{\pmb x}_2}, \cdots , {{\pmb x}_n}} \right]^{\rm T}} \in $ $ {\textbf{R}^{n \times m}} $, 构造局部结构保持函数$ {{\pmb J}_{\rm{local}}}\left( {\pmb \alpha} \right) $.
步骤2. 构造全局方差最大目标函数$ {{\pmb J}_{\rm{global}}}\left( {\pmb \alpha} \right) $.
步骤3. 构造整体目标函数$ {\pmb J}\left( {\pmb \alpha} \right) $.
步骤4. 根据式(14)求解特征值$ {\lambda _1}, {\lambda _2}, \cdots , {\lambda _n} $与对应的特征向量$ {{\pmb A}} = \left[ {{\alpha _1}, {\alpha _2}, \cdots , {\alpha _n}} \right] $.
步骤5. 根据式(15)求解前$ p $个特征值$ {\lambda _1} $, $ {\lambda _2} $, $ \cdots $, $ {\lambda _p} $与对应的特征向量$ {{\pmb A'}} = \left[ {{\alpha _1}, {\alpha _2}, \cdots , {\alpha _p}} \right] $.
步骤6. 根据公式$ {\pmb T} = {{\pmb K}^{\rm T}}{{\pmb A'}} $, 获得样本集在低维正交特征子空间的投影.
步骤7. 通过映射矩阵对训练及测试样本进行维数约简, 再将得到的低维特征子集输入到KNN, 并计算低维特征子集的$ {{\pmb S}_B} $, $ {{\pmb S}_W} $及$ {{\pmb S}_B}/{{\pmb S}_W} $指标.
3. 仿真实验与分析
3.1 柴油机故障数据仿真
船舶柴油机广泛应用于实际航运工程中, 其安全稳定的运行状态对整个系统起着至关重要的影响.因此, 在船舶柴油机发生故障时, 如果能够准确将故障信号的有效特征提取并分析, 则可提供足够多有效信息, 便于故障的分类和诊断[27-28].
1) 模型设计
本文以MAN公司S35ME-B9型柴油机为主要研究对象, 利用专业模拟软件AVL Boost完成柴油机故障模型仿真模拟系统, 图 2为柴油机仿真模型.
图 2中, SB1、SB2、SB3为系统边界, 外界气体通过SB1进入系统, 系统工质通过SB3排出系统, MP1$ \sim $MP8为测点, MP1和MP2测量中冷器前后的气体压力、温度, MP3和MP4测量气体进入和流出进气管PL1的气体压力、温度, MP5和MP6测量进入和流出1号缸的气体压力、温度, MP7和MP8测量废气进入和流出涡轮增压的气体压力、温度, C1$ \sim $C6为气缸, 1$ \sim $29为管道连接, PL1为进气管, CO1为中冷器, TC1为涡轮增压器.
本文对正常工况以及三种常见的船舶柴油机故障进行仿真模拟, 包括空冷器冷却不足, 排气口堵塞以及涡轮增压器效率降低, 如表 1所示, 由于本文采用的数值仿真模型, 因此采用设置关键参数的方式对故障进行模拟.在每种工况下, 记录模型中8个测量点的15个状态参数作为原始数据, 分别为功率(kW), 最大爆发压力(100 kPa), 压力机流量(kg/c), 压力机出口温度($ ^\circ $C), 压力机出口压力(100 kPa), 中冷器后温度($ ^\circ $C), 中冷器温差($ ^\circ $C), 中冷器后压力(100 kPa), 中冷器压差(100 kPa), 扫气温度($ ^\circ $C), 扫气压力(100 kPa), 排气管温度($ ^\circ $C), 排气管压力(100 kPa), 废气进涡轮机温度($ ^\circ $C), 涡轮增压出口温度($ ^\circ $C).
表 1 正常工况与故障工况模拟Table 1 The simulation of normal and fault conditionsNo. 工况类型 样本个数 数据维数 1 正常工况 960 15 2 故障1_空冷器冷却不足 960 15 3 故障2_排气口堵塞 960 15 4 故障3_涡轮增压效率降低 960 15 将故障数据集经本文所提降维方法进行处理, 选取KPCA、KLPP、核Fisher判别分析(Kernel fisher discriminant analysis, KFDA)[29]、局部和全局主成分分析(Local and global principal component analysis, LGPCA)[30]、全局–局部结构张量分析(Tensor global-local structure analysis, TGLSA)[20]和本文共6种算法进行对比, 在本研究中选取的目标维数为3维, 研究中采用了交叉验证方法选取最优高斯核参数, 实验从降维效果可视化、降维效果综合衡量指标和特征提取速度分析三方面验证方法的有效性.
2) 模型验证
使用AVL Boost进行柴油机工作状态仿真模型的建立, 选择台架实验的关键状态参数数据和AVL Boost对应状态参数进行对比验证.以柴油机功率、排气温度对柴油机模型正确性进行验证, 以排气阀堵塞和空冷器冷却不足等故障对柴油机模型进行故障模拟验证.
表 2是使用三种不同的工作状态的特定工作参数平均值与台架实验中对应参数的比较, 可以发现台架实验数据与AVL Boost相差较小, 可以认为使用AVL Boost建立的柴油机模型与台架实验使用柴油机具有相同的工作参数.
表 2 数据与台架实验数据多工况对比Table 2 The data contrast between AVL Boost and bench test under multiple working conditions负荷 排气温度(℃) 相对误差(%) 功率(kW) 相对误差(%) 模型数据 台架实验数据 模型数据 台架实验数据 90%负荷 329.89 328.50 0.42 3 281.40 3 277.00 0.13 75%负荷 304.39 307.30 0.95 2 839.20 2 844.00 0.17 75%推进 319.23 320.90 0.37 2 866.85 2 864.00 0.10 排气口堵塞是采用逐渐减小单个气缸的排气口直径的方法进行仿真验证, 排气口堵塞会造成排气效率不佳, 气缸内废气无法及时排出, 会造成气缸内气体逐渐增加, 缸内压力逐渐增大, 引起扫气压力增大, 进气逐渐减少, 也就是压力机流量减小.排气口堵塞还会造成气缸内高速积碳, 引起后燃现象, 使得排气温度上升, 燃烧效率下降, 但由于上述故障因素只是增加在一个气缸中, 其余五个气缸的燃烧过程影响不大, 因此对柴油机的功率影响不大.使用AVL Boost仿真这种故障, 上述提到的理论上的参数变化均获得了较好的验证.
空冷器冷却不足是使用逐渐增加冷却液温度, 逐渐降低空冷器冷却效率的方法仿真验证.这种故障最直观的反映就是空冷器前后温差降低, 另外由于空冷器冷却效率下降, 增压之后的气体无法较好地得到冷却, 扫气温度上升, 进入气缸内的气体质量会下降, 随之扫气压力上升, 流经压力机的空气流量会相应减小, 由于进入气缸的新鲜空气减少, 功率和最大爆发压力都会出现下降趋势, 同时由于进入气缸的空气温度上升, 排气温度也会上升.
3) 降维效果可视化
为验证本文降维方法的有效性, 将故障数据集经KPCA、KLPP、KFDA、LGPCA、TGLSA和本文所提降维方法进行处理, 选取前480个样本作为训练样本, 后480个样本作为测试样本.根据本文所提算法选取降维后的前三个主成分即可较为直观有效的表现降维效果, 图 3为根据贡献率原则, 选择不同降维维数与对应的主元贡献率, 由图中可知, 本文中选择主成分为前三维, 其贡献率之和即可达到总贡献率的85$ \% $, 因此仿真实验中的贡献率选择为85$ \% $即可.得到降维后的测试样本三维特征量分布见图 4~9.
从图 4可以看出, KPCA可以分辨正常工况与故障1, 但对于故障2和故障3有较严重的数据重叠现象.从图 5可以看出, KLPP对于4种工况均不能较为有效地进行区分.
图 6~8分别为KFDA、LGPCA、TGLSA三种不同降维方法的对比实验效果图, KFDA效果较差, 样本数据大部分呈现混叠状态, LGPCA和TGLSA算法也综合保持了全局结构和局部结构, 因此其效果略好于KFDA、LGPCA和TGLSA三种算法, 但由于其主要应用于维数较低的数据, 因此对于船舶柴油机的高维复杂数据, 效果一般.
本文所提降维方法效果图如图 9所示, 其中, 权重参数根据计算得到$ {\beta _{\pmb C}} = 0.79 $, $ {\beta _{{\pmb L}'}} = 0.21 $, 可见对于柴油机故障数据, KPCA方法所侧重的全局特征占主导地位.由于全局和局部结构特征提取过程的综合考虑, 在完成数据约简和可视化的同时, 有效地分离了四种状态, 同时具有良好的聚类能力.因此, 本文提出的方法能够提取故障特征, 解决了数据降维可视化问题.
4) 降维效果综合衡量指标
为直观有效地可视化各类方法的降维效果, 将类间距$ {{\pmb S}_B} $, 类内距$ {{\pmb S}_W} $及二者的比值$ {{\pmb S}_B}/{{\pmb S}_W} $作为衡量指标, 类间距$ {{\pmb S}_B} $及二者的比值$ {{\pmb S}_B}/{{\pmb S}_W} $越大, 表明分类效果明显, 得到的评价结果见图 10, 从图中可以看出, 六种方法中, 本文所提GLSP方法对应的类间距及二者的比值具有最大值.分析如下:
a) KPCA降维效果及识别率一般, 虽然可以去除特征空间中的数据冗余信息, 但并未达到能够有效表达最佳情况的条件; KLPP局部结构保持的方法的类间距明显提高, 该方法能有效提取出数据集中局部结构, 但对于综合全局与结构方法所得到的类间距与类内距的比值, 还有一定差距; KFDA方法具有较高类内距, 但是对于类内距没有很好的聚合效果; LGPCA和TGLSA两种方法均具有较好类内聚合作用, 但是由于不同类别的类间距离较小, 导致综合类间距与类内距比值较小.
b) 本文所提方法的降维效果及识别率要高于其他对比方法, 该方法能够避免子空间重构, 更利于故障类别的划分, 且具有较强的全局与局部判别信息的挖掘能力.
故障诊断的实质是模式识别, 考虑船舶柴油机故障诊断的实船应用性, 选择极限学习机(Extreme learning machine, ELM)、支持向量机(Support vector machine, SVM)、相关向量机(Relevance vector machine, RVM)、KNN四种基础有效的分类方法进行检验, 将低维样本分别应用于上述四种方法进行分类效果比较, 表 3~5给出了测试样本的故障诊断结果, 可知相较传统的分类方法, 对于仿真所得故障数据的诊断率均不高, 但是所提出的GLSP降维方法所获得的低维有效数据在大部分情况下获得最高的故障识别精度, 未获得最高识别精度情况下, 其精度与最高精度相差不大.
表 3 故障1识别准确率($ \% $)Table 3 The accuracy of fault1 diagnosis ($ \% $)方法 Fault1 KPCA KLPP KFDA LGPCA TGLSA GLSP ELM 55.32 61.38 60.58 54.21 58.69 62.97 SVM 58.69 70.61 71.68 65.34 68.49 69.27 RVM 72.77 69.59 74.21 68.98 63.40 76.35 KNN 72.26 66.86 70.38 75.49 77.36 78.53 表 4 故障2识别准确率($ \% $)Table 4 The accuracy of fault2 diagnosis ($ \% $)方法 Fault2 KPCA KLPP KFDA LGPCA TGLSA GLSP ELM 80.95 76.85 79.65 77.49 70.28 82.62 SVM 78.36 77.32 77.05 74.39 72.15 80.09 RVM 79.74 74.16 78.66 85.68 81.29 83.62 KNN 82.35 82.63 75.39 78.91 86.54 88.84 表 5 故障3识别准确率($ \% $)Table 5 The accuracy of fault3 diagnosis ($ \% $)方法 Fault3 KPCA KLPP KFDA LGPCA TGLSA GLSP ELM 70.65 72.39 77.16 74.29 70.53 79.26 SVM 66.34 68.29 68.49 65.39 60.87 66.58 RVM 59.38 62.58 55.21 59.86 60.13 66.34 KNN 58.62 62.38 65.98 63.24 61.09 65.08 5) 降维效果综合衡量指标
时间性能比较是在实验室中台式机电脑上进行, 其配置为Intel Core i3 CPU 3.3GHz, RAM 4GB, Win7操作系统, 仿真软件为MATLAB 2010a, 计算结果如表 6所示.由表 6可知, 本文所提方法降维所需时间相对于其他方法有所增加, 但是都是在一个数量级, 且运行速度均在5s以内, 满足实际情况中对于实时的要求.
表 6 特征提取所需时间(s)Table 6 Feature extraction time (s)维度 特征提取方法 KPCA KLPP KFDA LGPCA TGLSA GLSP 3 0.651 1.155 1.039 2.598 2.134 1.596 5 0.795 1.159 1.118 2.019 1.495 1.632 8 0.815 1.209 0.975 1.069 1.396 1.885 10 0.867 1.344 1.185 1.563 2.098 1.962 3.2 TE数据仿真
TE过程是一个公认的对比各种控制和监控方案的平台, 为验证本文所提方法的通用性, 将本文算法应用于故障检测与诊断领域被广泛使用的TE化工数据集上. TE过程是基于真实工业过程的仿真平台, 包含了一组正常状态和21组不同故障状态, 分别涵盖了12个操纵变量和41个测量变量, 每组状态包含480组训练数据和960组测试数据, 每一组故障从第160个数据点引入, 过程的详细描述、工艺流程图以及其故障形式的具体介绍见文献[31].
从TE数据集中选取1组正常工况数据和3组故障数据(故障4、故障8和故障14), 使用本文所提GLSP算法得到降维后的前两维特征矢量, 对于TE数据, 前两维即可较好体现降维效果, 本文权重参数根据计算得到$ {\beta _{\pmb C}} = 0.64, {\beta _{{\pmb L}'}} = 0.36 $, 并将本文算法与KPCA、KLPP、KFDA、LGPCA和TGLSA五种降维方法进行比较, 如图 11至图 16.
从图 11和图 12前两维降维效果可以看出, KPCA和KLPP对于4类数据分离效果较差, 特别是正常工况与故障4和故障14存在大量重合现象; 图 13和图 14中KFDA和LGPCA方法对于正常工况和故障14难以有效区分; 图 15和图 16的特征提取效果略好于前几种方法, 但TGLSA方法虽然能够较好地区分故障4和故障14, 但是正常工况和故障8仍有较大重叠, 而本文所提算法提取了更为丰富的全局与局部结构信息, 在完成数据降维可视化的同时能够将四类数据有效分离.
4. 结论
本文通过融合全局特征提取KPCA与局部特征提取KLPP两种降维方法, 提出一种结合两种降维方法的数据维数约简方法GLSP, 增强了数据低维可视化效果, 同时提高了识别精度, 并将其应用于故障诊断中.所提方法将流形学习保持局部结构的思想融入核主成分分析的目标函数中, 使得到的特征空间不仅具有原始样本空间的整体结构, 还保持样本空间相似的局部近邻结构, 可以包含更丰富的特征信息.使用AVL Boost软件对船舶柴油机工作过程进行故障仿真, 提取正常工况与故障工况下的仿真数据, 并将AVL Boost软件仿真数据和TE化工公共故障数据应用到所提方法中, 实验结果证明, 本文所提GLSP算法具有较好的维度约简效果, 并具有较高分类精度的优势.
-
表 1 线性反馈参数设计
Table 1 Parameter design for linear feedback
参数 血糖上升阶段 血糖下降阶段 ${y \le 275\;{\rm{mg/dL} } }$ ${ {y > 275}\;{\rm{mg/dL} } }$ k1(y) 0.0001 × exp(h) 0.00001 × exp(−h) 0.00015 × exp(h) k2(y) 0.002 × exp(h) 0.001 × exp(h) 0.002 × exp(h) k3 1 0.15 1 注: $h = (y - {G_r})/y$ 表 2 正常胰岛素基础率下本控制算法与zMPC控制算法的评估结果
Table 2 Evaluation results of the proposed controller and the zMPC controller at normal insulin basal rate
血糖控制性能指标 餐前补充大剂量胰岛素 餐前未补充大剂量胰岛素 zMPC 本算法 zMPC 本算法 < 54 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) < 70 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) 70 ~ 180 mg/dL 时间比率 (%) 92.4 (10.2) 92.5 (7.8) 72.4 (13.4) 71.8 (9.4) > 250 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 1.2 (9.5) 2.5 (9.2) 血糖平均值 (mg/dL) 135.4 (6.1) 133.9 (6.4) 151.7 (19.9) 152.2 (13.9) 标准差 (mg/dL) 28.2 (8.2) 28.0 (4.9) 44.9 (14.1) 45.8 (14.6) 7:00 血糖平均值 (mg/dL) 121.0 (15.2) 116.0 (15.5) 121 (15.5) 118.7 (15.0) 注: 统计数据以中位数 (四分位距) 的形式表示. 表 4 2倍正常胰岛素基础率下本控制算法与zMPC控制算法的评估结果
Table 4 Evaluation results of the proposed controller and the zMPC controller at 2 times of normal insulin basal rate
血糖控制性能指标 餐前补充大剂量胰岛素 餐前未补充大剂量胰岛素 zMPC 本算法 zMPC 本算法 < 54 mg/dL 时间比率 (%) 7.9 (10.0) 0 (0.0) 7.8 (10.9) 0 (0.0) < 70 mg/dL 时间比率 (%) 19.9 (9.3) 0 (0.8) 18.4 (9.2) 0 (0.0) 70 ~ 180 mg/dL 时间比率 (%) 77.4 (10.2) 93.4 (6.2) 67.8 (12.8) 74.4 (12.6) > 250 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 0 (3.9) 2.2 (9.0) 血糖平均值 (mg/dL) 101.2 (6.2) 130.2 (6.8) 113.5 (14.4) 149.5 (11.9) 标准差 (mg/dL) 35.7 (6.8) 28.3 (5.8) 50.8 (11.0) 45.2 (14.1) 7:00 血糖平均值 (mg/dL) 83.0 (38.0) 111.0 (26.0) 78.5 (35.0) 113.0 (27.0) 注: 统计数据以中位数 (四分位距) 的形式表示. 表 3 0.5倍正常胰岛素基础率下本控制算法与zMPC控制算法的评估结果
Table 3 Evaluation results of the proposed controller and the zMPC controller at 0.5 times of normal insulin basal rate
血糖控制性能指标 餐前补充大剂量胰岛素 餐前未补充大剂量胰岛素 zMPC 本算法 zMPC 本算法 < 54 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) < 70 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) 70 ~ 180 mg/dL 时间比率 (%) 66.6 (9.7) 92.0 (7.4) 38.5 (42.1) 69.6 (9.3) > 250 mg/dL 时间比率 (%) 0 (0.0) 0 (0.0) 18.5 (27.9) 5.0 (11.6) 血糖平均值 (mg/dL) 171.8 (5.9) 135.5 (6.9) 204.9 (49.0) 159.8 (11.7) 标准差 (mg/dL) 24.1 (5.9) 27.7 (5.5) 43.2 (13.5) 45.8 (15.6) 7:00 血糖平均值 (mg/dL) 164.0 (26.5) 117.0 (26.0) 167.5 (28.0) 124.0 (27.0) 注: 统计数据以中位数 (四分位距) 的形式表示. -
[1] 纪立农, 孙子林, 李启富, 秦贵军, 徐海岩. 中国四城市糖尿病患者胰岛素注射相关皮下脂肪增生的横断面研究. 中国糖尿病杂志, 2019, 27(10): 721−727 doi: 10.3969/j.issn.1006-6187.2019.10.001Ji Li-Nong, Sun Zi-Lin, Li Qi-Fu, Qin Gui-Jun, Xu Hai-Yan. Cross-sectional study of insulin injection-related lipohypertrophy in diabetic patients in four cities in China. Chinese Journal of Diabetes, 2019, 27(10): 721−727 doi: 10.3969/j.issn.1006-6187.2019.10.001 [2] 纪立农. 活的更长、活得更好——糖尿病与“健康中国2030”. 中国糖尿病杂志, 2019, 27(1): 1−2 doi: 10.3969/j.issn.1006-6187.2019.01.001Ji Li-Nong. Live longer and live better—diabetes and “Healthy China 2030”. Chinese Journal of Diabetes, 2019, 27(1): 1−2 doi: 10.3969/j.issn.1006-6187.2019.01.001 [3] Liu W, McGuire H C, Kissimova-Skarbek K, Zhou X H, Han X Y, Wang Y N, et al. Factors associated with acute complications among individuals with type 1 diabetes in China: The 3C study. Endocrine Research, 2020, 45(1): 1−8 doi: 10.1080/07435800.2019.1624567 [4] Liu W, Wang Y N, Han X Y, Cai X L, Zhu Y, Zhang M X, et al. Factors associated with resistance to complications in long-standing type 1 diabetes in China. Endocrine Connections, 2020, 9(2): 187−193 doi: 10.1530/EC-19-0521 [5] Beck R W, Bergenstal R M, Laffel P L M, Pickup P J C. Advances in technology for management of type 1 diabetes. The Lancet, 2019, 394(10205): 1265−1273 doi: 10.1016/S0140-6736(19)31142-0 [6] Benhamou P Y, Franc S, Reznik Y, Thivolet C, Schaepelynck P, Renard E, et al. Closed-loop insulin delivery in adults with type 1 diabetes in real-life conditions: A 12-week multicentre, open-label randomised controlled crossover trial. The Lancet Digital Health, 2019, 1(1): e17−e25 doi: 10.1016/S2589-7500(19)30003-2 [7] Weisman A, Bai J W, Cardinez M, Kramer C K, Perkins B A. Effect of artificial pancreas systems on glycaemic control in patients with type 1 diabetes: A systematic review and meta-analysis of outpatient randomised controlled trials. The Lancet Diabetes & Endocrinology, 2017, 5(7): 501−512 [8] Davis G M, Galindo R J, Migdal A L, Umpierrez G E. Diabetes technology in the inpatient setting for management of hyperglycemia. Endocrinology and Metabolism Clinics of North America, 2020, 49(1): 79−93 doi: 10.1016/j.ecl.2019.11.002 [9] Quiroz G. The evolution of control algorithms in artificial pancreas: A historical perspective. Annual Reviews in Control, 2019, 48: 222−232 doi: 10.1016/j.arcontrol.2019.07.004 [10] Kovatchev B. A century of diabetes technology: Signals, models, and artificial pancreas control. Trends in Endocrinology and Metabolism, 2019, 30(7): 432−444 [11] Patra K, Nanda A, Panigrahi S, Mishra A K. The fractional order PID controller design for BG control in type-I diabetes patient. Advances in Intelligent Computing and Communication. Singapore: Springer, 2020. 321−329 [12] Ly T T, Keenan D B, Roy A, Han J, Grosman B, Cantwell M, et al. Automated overnight closed-loop control using a proportional-integral-derivative algorithm with insulin feedback in children and adolescents with type 1 diabetes at diabetes camp. Diabetes Technology & Therapeutics, 2016, 18(6): 377−384 [13] Forlenza G P, Buckingham B A, Christiansen M P, Wadwa R P, Peyser T A, Lee J B, et al. Performance of omnipod personalized model predictive control algorithm with moderate intensity exercise in adults with type 1 diabetes. Diabetes Technology & Therapeutics, 2019, 21(5): 265−272 [14] Villa-Tamayo M F, Rivadeneira P S. Adaptive impulsive offset-free MPC to handle parameter variations for type 1 diabetes treatment. Industrial & Engineering Chemistry Research, 2020, 59(13): 5865−5876 [15] Garcia-Tirado J, Colmegna P, Corbett J P, Ozaslan B, Breton M D. In silico analysis of an exercise-safe artificial pancreas with multistage model predictive control and insulin safety system. Journal of Diabetes Science and Technology, 2019, 13(6): 1054−1064 doi: 10.1177/1932296819879084 [16] Bahremand S, Ko H S, Balouchzadeh R, Lee H F, Park S, Kwon G. Neural network-based model predictive control for type 1 diabetic rats on artificial pancreas system. Medical & Biological Engineering & Computing, 2019, 57(1): 177−191 [17] Huyett L M, Dassau E, Zisser H C, Doyle III F J. Design and evaluation of a robust PID controller for a fully implantable artificial pancreas. Industrial & Engineering Chemistry Research, 2015, 54(42): 10311−10321 [18] 韩月起, 张凯, 宾洋, 秦闯, 徐云霄, 李小川, 等. 基于凸近似的避障原理及无人驾驶车辆路径规划模型预测算法. 自动化学报, 2020, 46(1): 153−167Han Yue-Qi, Zhang Kai, Bin Yang, Qin Chuang, Xu Yun-Xiao, Li Xiao-Chuan, et al. Convex approximation based avoidance theory and path planning MPC for driver-less vehicles. Acta Automatica Sinica, 2020, 46(1): 153−167 [19] 魏永涛, 高原, 孙文义, 王秀蒙. 交通流动态扰动下的区域交通信号协调控制. 自动化学报, 2019, 45(10): 1983−1994Wei Yong-Tao, Gao Yuan, Sun Wen-Yi, Wang Xiu-Meng. Regional traffic signal control considering the dynamic characteristics of traffic flow. Acta Automatica Sinica, 2019, 45(10): 1983−1994 [20] 姜頔, 刘向杰, 孔小兵. 核电站蒸汽发生器水位的软约束预测控制. 自动化学报, 2019, 45(6): 1111−1121Jiang Di, Liu Xiang-Jie, Kong Xiao-Bing. Soft constrained MPC on water level control in steam generator of a nuclear power plant. Acta Automatica Sinica, 2019, 45(6): 1111−1121 [21] Boiroux D, Duun-Henriksen A K, Schmidt S, NΦrgaard K, Poulsen N K, Madsen H, et al. Adaptive control in an artificial pancreas for people with type 1 diabetes. Control Engineering Practice, 2017, 58: 332−342 doi: 10.1016/j.conengprac.2016.01.003 [22] Incremona G P, Messori M, Toffanin C, Cobelli C, Magni L. Artificial pancreas: From control-to-range to control-to-target. IFAC-PapersOnLine, 2017, 50(1): 7737−7742 doi: 10.1016/j.ifacol.2017.08.1152 [23] Shi D W, Dassau E, Doyle III F J. Adaptive zone model predictive control of artificial pancreas based on glucose-and velocity-dependent control penalties. IEEE Transactions on Biomedical Engineering, 2019, 66(4): 1045−1054 doi: 10.1109/TBME.2018.2866392 [24] Palerm C C, Zisser H, Jovanovič L, Doyle III F J. A run-to-run control strategy to adjust basal insulin infusion rates in type 1 diabetes. Journal of Process Control, 2008, 18(3−4): 258−265 doi: 10.1016/j.jprocont.2007.07.010 [25] Toffanin C, Visentin R, Messori M, Di Palma F, Magni L, Cobelli C. Toward a run-to-run adaptive artificial pancreas: In silico results. IEEE Transactions on Biomedical Engineering, 2018, 65(3): 479−488 doi: 10.1109/TBME.2017.2652062 [26] Shi D W, Dassau E, Doyle III F J. Multivariate learning framework for long-term adaptation in the artificial pancreas. Bioengineering & Translational Medicine, 2019, 4(1): 61−74 [27] Wang Y Q, Dassau E, Doyle III F J. Closed-loop control of artificial pancreatic β-cell in type 1 diabetes mellitus using model predictive iterative learning control. IEEE Transactions on Biomedical Engineering, 2010, 57(2): 211−219 doi: 10.1109/TBME.2009.2024409 [28] Wang Y Q, Zhang J P, Zeng F M, Wang N, Chen X P, Zhang B, et al. “Learning” can improve the blood glucose control performance for type 1 diabetes mellitus. Diabetes Technology and Therapeutics, 2017, 19(1): 41−48 [29] Han J P. From PID to active disturbance rejection control. IEEE Transactions on Industrial Electronics, 2009, 56(3): 900−906 doi: 10.1109/TIE.2008.2011621 [30] 杨飞, 谈树萍, 薛文超, 郭金, 赵延龙. 饱和约束测量扩张状态滤波与无拖曳卫星位姿自抗扰控制. 自动化学报, 2020, 46(11): 2337−2349Yang Fei, Tan Shu-Ping, Xue Wen-Chao, Guo Jin, Zhao Yan-Long. Extended state filtering with saturation-constrainted observations and active disturbance rejection control of position and attitude for drag-free satellites. Acta Automatica Sinica, 2020, 46(11): 2337−2349 [31] 刘昊, 魏承, 谭春林, 刘永健, 赵阳. 空间充气展开绳网系统捕获目标自抗扰控制研究. 自动化学报, 2019, 45(9): 1691−1700Liu Hao, Wei Cheng, Tan Chun-Lin, Liu Yong-Jian, Zhao Yang. Research on capturing target of space inflatable net capture system based on active disturbance rejection control. Acta Automatica Sinica, 2019, 45(9): 1691−1700 [32] Huang Y, Xue W C. Active disturbance rejection control: Methodology and theoretical analysis. ISA Transactions, 2014, 53(4): 963−976 doi: 10.1016/j.isatra.2014.03.003 [33] Feng H, Guo Z U. Active disturbance rejection control: Old and new results. Annual Reviews in Control, 2017, 44: 238−248 doi: 10.1016/j.arcontrol.2017.05.003 [34] Guo B Z, Zhao Z L. Active disturbance Rejection Control for Nonlinear Systems: An Introduction. Hoboken, NJ: John Wiley & Sons, 2016. [35] Zhao Z L, Guo B Z. On convergence of nonlinear active disturbance rejection control for SISO nonlinear systems. Journal of Dynamical and Control Systems, 2016, 22(2): 385−412 doi: 10.1007/s10883-015-9304-5 [36] Gao Z Q. Scaling and bandwidth-parameterization based controller tuning. In: Proceedings of the 2003 American Control Conference. Denver, CO, USA: IEEE, 2003. [37] Pu Z Q, Yuan R Y, Yi J Q, Tan X M. A class of adaptive extended state observers for nonlinear disturbed systems. IEEE Transactions on Industrial Electronics, 2015, 62(9): 5858−5869 doi: 10.1109/TIE.2015.2448060 [38] Guo B Z, Zhao Z L. On convergence of the nonlinear active disturbance rejection control for MIMO systems. SIAM Journal on Control and Optimization, 2013, 51(2): 1727−1757 doi: 10.1137/110856824 [39] Keith-Hynes P, Guerlain S, Mize B, Hughes-Karvetski C, Khan M, McElwee-Malloy M, et al. DiAs user interface: A patient-centric interface for mobile artificial pancreas systems. Journal of Diabetes Science and Technology, 2013, 7(6): 1416−1426 doi: 10.1177/193229681300700602 [40] Lazaro C, Oruklu E, Sevil M, Turksoy K, Cinar A. Implementation of an artificial pancreas system on a mobile device. In: Proceedings of the 2016 IEEE International Conference on Electro Information Technology (EIT). Grand Forks, USA: IEEE, 2016. 642−647 [41] Deshpande S, Pinsker J E, Zavitsanou S, Shi D W, Tompot R, Church M M, et al. Design and clinical evaluation of the interoperable artificial pancreas system (iAPS) smartphone app: Interoperable components with modular design for progressive artificial pancreas research and development. Diabetes Technology & Therapeutics, 2019, 21(1): 35−43 [42] Man C D, Micheletto F, Lv D, Breton M, Kovatchev B, Cobelli C. The UVA/PADOVA type 1 diabetes simulator: New features. Journal of Diabetes Science and Technology, 2014, 8(1): 26−34 doi: 10.1177/1932296813514502 [43] Gondhalekar R, Dassau E, Doyle III F J. Periodic zone-MPC with asymmetric costs for outpatient-ready safety of an artificial pancreas to treat type 1 diabetes. Automatica, 2016, 71: 237−246 doi: 10.1016/j.automatica.2016.04.015 [44] Gondhalekar R, Dassau E, Doyle F J. Velocity-weighting and velocity-penalty MPC of an artificial pancreas: Improved safety and performance. Automatica, 2018, 91: 105−117 doi: 10.1016/j.automatica.2018.01.025 [45] Cai D H, Song J L, Wang J Z, Shi D W. Glucose regulation for subjects with type 1 diabetes using active disturbance rejection control. In: Proceedings of the 2019 Chinese Control Conference (CCC). Guangzhou, China: IEEE, 2019. 6970−6975 [46] 韩京清. 自抗扰控制技术——估计补偿不确定因素的控制技术. 北京: 国防工业出版社, 2008.Han Jing-Qing. Active Disturbance Rejection Control Technique: The Technique for Estimating and Compensating the Uncertainties. Beijing: National Defense Industry Press, 2008. [47] Levine W S. The Control Handbook (2nd ed.) Boca Raton, USA: CRC Press, 2011. [48] Huang Y, Wang J Z, Shi D W, Wu J F, Shi L. Event-triggered sampled-data control: An active disturbance rejection approach. IEEE/ASME Transactions on Mechatronics, 2019, 24(5): 2052−2063 doi: 10.1109/TMECH.2019.2933411 [49] Huang Y, Wang J F, Shi D W, Shi L. Toward event-triggered extended state observer. IEEE Transactions on Automatic Control, 2018, 63(6): 1842−1849 doi: 10.1109/TAC.2017.2754340 -