本发明公开了一种基于凸空间滤波的工业系统故障诊断方法，属于工业系统故障诊断领域。本方法根据系统的线性模型，获取系统的测量带信息和参数真值向量；根据系统的参数真值向量，确定包裹参数真值向量的初始化凸空间；根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障；若系统没有发生故障，则求取k时刻的凸空间和正交凸空间，若系统发生故障，由凸空间测试集的状态立马隔离出故障参数，在仅扩容故障参数而不影响其他参数的滤波估计区间基础上，完成故障识别。该方法不仅降低了传统凸空间滤波方法的保守性，而且提升了系统的故障诊断精度。
1.一种基于凸空间滤波的工业系统故障诊断方法，其特征在于，所述方法包括：S1：确定工业系统的线性离散模型；
S2：根据S1确定的线性离散模型，获取工业系统的参数真值向量和k时刻的测量带信息；
S3：根据S2获取的工业系统的参数真值向量，确定包裹参数真值向量的初始化凸空间，将此初始化凸空间作为初始时刻k＝1的凸空间；
S4：k≠1时刻，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障；
若系统没有发生故障，则求取k时刻的凸空间和正交凸空间，所述正交凸空间为k‑1时刻的凸空间和k时刻的凸空间的交集；
若系统发生故障，则扩容k‑1时刻的正交凸空间，进而计算得到凸空间测试集；
对具体故障参数进行扩容得到包裹故障真值向量的k时刻的正交凸空间，将k时刻的正交凸空间作为k时刻的凸空间，重复S4，完成故障辨识。
2.根据权利要求1所述的方法，其特征在于，所述工业系统为浆距子系统时，所述S1包括：获取浆距子系统的状态空间模型为：
其中β和βa分别为桨距角和角速度大小， 和 分别为桨距角的一阶导数和角速度的一阶导数，βr为桨距的参考值，ζ和ωn分别为阻尼系数和浆距子系统的固有频率；
将式(1)的状态空间模型转化为二阶连续系统模型：其中y为输出，表示桨距角；u为输入，表示桨距的参考值，即y＝β，u＝βr，利用采样时刻Ts，将式(2)离散化为：A(z)y(k)＝B(z)u(k)+e(k) (3)‑1 ‑2 ‑1 ‑2
其中A(z)和B(z)为z变换多项式，A(z)＝1+a1z +a2z ,B(z)＝b1z +b2z ；u(k)为k时刻系统输入数据，即桨距的参考值；则由式(3)获取浆距子系统的线性离散模型为：其中k表示时刻， 为可观测的数据向量，表示浆距子系统中由桨距角历史值和桨距的参考值组成的信息向量，y(k)为k时刻系统输出数据，即当前时刻的桨距角，上标T表示向量的转置，θt为浆距子系统T参数真值向量，θt＝[a1,a2,b1,b2] ；e(k)为k时刻浆距子系统受到的噪声，|e(k)|≤δ(k)，δ(k)为浆距子系统在k时刻受到噪声绝对值的最大值。
3.根据权利要求2所述的方法，其特征在于，所述S2包括：T
由式(4)可知θt＝[a1,a2,b1,b2]为系统参数真值向量；
其中，ck为k时刻测量带信息的中心点， pk为k时刻测量带信息的方向向量，θs为测量带信息估计的包裹参数真值的范围。
4.根据权利要求3所述的方法，其特征在于，所述S3包括：根据参数真值向量，确定初始化凸空间 以包裹参数真值向量，令初始时刻k＝
其中，T为凸空间的生成矩阵，θc为凸空间的中心，θ为中间向量，‖·‖∞表示无穷范数，θ表示凸空间包裹范围内所有的参数可行集。
5.根据权利要求4所述的方法，其特征在于，所述S4中，k≠1时，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障，若交集为空，系统发生故障，即标量 或者 若交集不为空，系统没有发生故障，即标量 且
其中，p0，c0为中间变量，p0＝pk，c0＝ck，θc,k‑1为k‑1时刻凸空间的中心，tk‑1,h为k‑1时刻凸空间生成矩阵的第h列，h为整数，h＝1,...,n。
6.根据权利要求5所述的方法，其特征在于，所述S4中，当k时刻的测量带信息与k‑1时刻的凸空间的交集不为空时，系统没有发生故障，即式(7)中的标量 且 则求取k时刻包裹参数真值向量的凸空间 和正交凸空间
由式(7)中的标量 和 求得中间变量 和 即得到缩减后的测量带信息 为缩减后的测量带信息的方向向量，为缩减后的测量带信息的中心：
引入中间变量 得到缩减后的凸空间 为缩减后的凸空间的生成矩阵， 为缩减后的凸空间的生成矩阵的第g列，g为整数，为缩减后的凸空间的中心：则 可有以下公式获得
为中间凸空间生成矩阵T的第w列向量，w为整数，w＝1,...,n； 为中间凸空间的中心：由式(7)～(15)求出包裹参数真值向量的凸空间 后，利用以下式(16)～(20)可求出正交凸空间
定义n维凸空间生成矩阵的每一列为tj，其中j＝1,…,n，计算出n维凸空间的顶点Vl：n n
其中l为顶点的个数，l＝1,2,…,2，αl,j为中间变量，αl,j∈{‑1,1}；R表示n维实数集；
定义k‑1时刻n维凸空间的顶点矩阵为 凸空间有多少个顶点，顶点矩阵就有多少列，Vk‑1(p,l)表示k‑1时刻Vk‑1矩阵第l列的第p个元素，p为整数，p＝
由式(16)求出k‑1时刻n维凸空间中每个参数的边界最大值 和边界最小值将k时刻与k‑1时刻的凸空间的参数边界最大值对比取最小，参数边界最小值对比取最大，得到其中，v为整数，v＝1,…,n， 为中间变量；
其中，θc,k为k时刻正交凸空间的中心，Tk为k时刻正交凸空间的生成矩阵，‖θ‖∞为中间变量；
当系统在无故障情况下，重复式(7)～(20)可求出每个时刻包裹真值向量的凸空间和单调收敛的正交凸空间。
7.根据权利要求6所述的方法，其特征在于，所述S4中，若系统发生故障，则将k‑1时刻的正交凸空间 扩容获得n个正交凸空间测试集 其中i＝1,…,n；
其中，f为整数， 为第i个正交凸空间测试集里的第f个参数边界最大值， 为第i个正交凸空间测试集里的第f个参数边界最小值， 为正交凸空间第f个参数扩容的最大幅值；
给定时间长度L，在k时刻检测到发生故障后的L时刻内无新故障发生，获得测量带信息S(k),…,S(k+L)；利用式(7)～(15)将扩容得到的n个正交凸空间测试集分别与L个测量带信息计算获得k+L时的n个凸空间测试集 q为整数，q＝1,...,n。
8.根据权利要求7所述的方法，其特征在于，所述根据凸空间测试集的状态隔离出系统具体的故障参数，对具体故障参数进行扩容得到包裹故障真值向量的正交凸空间，完成故障识别，包括：若 为空集，说明参数真值向量的第q个分量发生故障；
若 不为空集，说明参数真值向量的第q个分量没有发生故障；
当在k时刻隔离出系统的具体故障参数后，扩容k‑1时刻的正交凸空间，作为k时刻包裹故障真值向量的正交凸空间 即仅在故障参数方向上进行扩容，假设参数真值向量的第q个分量发生故障，包裹故障真值向量的正交凸空间为：其中
其中，d为整数，d＝1,...,n； 为正交凸空间里的第d个参数边界最大值，为正交凸空间里的第d个参数边界最小值； 为正交凸空间第d个参数扩容的最大幅值；
扩容后得到的正交凸空间 是包裹故障参数真值向量的，将正交凸空间作为当前时刻的凸空间，利用上式(7)～(15)求出包裹参数真值向量的凸空间以及式(16)～(20)求出正交凸空间，不断迭代至求出的正交凸空间在预定空间阈值范围内，完成故障识别。
9.根据权利要求8所述的方法，其特征在于，所述方法应用于浆距子系统的故障检测时，取ζ＝0.6，ωn＝11.11rad/s，采样时间Ts＝0.01s。
10.根据权利要求9所述的方法，所述工业系统包括：浆距子系统、直流电动机系统、弹簧阻尼系统。
一种基于凸空间滤波的工业系统故障诊断方法
技术领域
[0001] 本发明涉及一种基于凸空间滤波的工业系统故障诊断方法，属于工业系统故障诊断领域。
背景技术
[0002] 随着科学与技术的飞速发展，现代化复杂工程系统规模越来越大，例如动力系统工程、冶金石化工业和电网系统等，这些现代工业的特点是大型化、集成化、自动化，对设备依赖程度高，因此及时快速的对工业系统进行故障诊断至关重要。
[0003] 故障诊断主要分为故障检测、故障隔离和故障识别三个环节。目前基于故障诊断的方法大体可分为三类：基于模型的方法、基于信号处理的方法和基于知识的方法，其中，基于信号处理的方法对早期潜在的故障无效，只有当故障发生到一定的程度并影响到外部特征时才有效，且多数情况下无法定位故障；基于知识的方法则需要知道大量工业对象的专业知识进行故障诊断，而基于模型的故障诊断方法研究最早，最为深入，同时最为成熟，因为有已知的数学模型，稳定性好，实时性好。
[0004] 但是传统的基于模型的故障诊断方法要求输入噪声必须满足一定的分布规律，比如高斯分布等，然而，实际工业系统领域里的噪声往往是无法测量的，仅能得到噪声的最大幅值，因此有学者提出了一种基于椭球空间滤波的故障诊断方法(参见Reppa V,Tzes A.Fault diagnosis based on set membership identification using outputerror models[J].International Journal of Adaptive Control and Signal Processing,2016,30(2):224255.)，但是鉴于一些工业系统对于故障诊断的更高要求，仍然需要对工业系统的故障诊断进行更深的研究。
发明内容
[0005] 为了进一步提高对于工业系统故障诊断的速度、获得保守型更低、故障估计区间更小的故障估计结果，本发明提供了一种基于凸空间滤波的工业系统故障诊断方法，所述方法包括：
[0006] S1：确定工业系统的线性离散模型；
[0007] S2：根据S1确定的线性离散模型，获取工业系统的参数真值向量和k时刻的测量带信息；
[0008] S3：根据S2获取的工业系统的参数真值向量，确定包裹参数真值向量的初始化凸空间，将此初始化凸空间作为初始时刻k＝1的凸空间；
[0009] S4：k≠1时刻，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障；
[0010] 若系统没有发生故障，则求取k时刻的凸空间和正交凸空间，所述正交凸空间为k‑1时刻的凸空间和k时刻的凸空间的交集；
[0011] 若系统发生故障，则扩容k‑1时刻的正交凸空间，进而计算得到凸空间测试集；
[0012] 根据凸空间测试集的状态隔离出系统的具体故障参数；
[0013] 对具体故障参数进行扩容得到包裹故障真值向量的k时刻的正交凸空间，将k时刻的正交凸空间作为k时刻的凸空间，重复S4，完成故障辨识。
[0014] 可选的，所述工业系统为浆距子系统时，所述S1包括：
[0015] 获取浆距子系统的状态空间模型为：
[0016]
[0017] 其中β和βa分别为桨距角和角速度大小， 和 分别为桨距角的一阶导数和角速度的一阶导数，βr为桨距的参考值，ζ和ωn分别为阻尼系数和浆距子系统的固有频率；
[0018] 将式(1)的状态空间模型转化为二阶连续系统模型：
[0019]
[0020] 其中y为输出，表示桨距角；u为输入，表示桨距的参考值，即y＝β，u＝βr[0021] 利用采样时刻Ts，将式(2)离散化为：
[0022] A(z)y(k)＝B(z)u(k)+e(k) (3)
[0023] 其中A(z)和B(z)为z变换多项式，A(z)＝1+a1z‑1+a2z‑2,B(z)＝b1z‑1+b2z‑2；u(k)为k时刻系统输入数据，即桨距的参考值；则由式(3)获取浆距子系统的线性离散模型为：
[0024]
[0025] 其中k表示时刻，k＝1,…,N， 为可观测的数据向量，表示浆距子系统中由桨距角历史值和桨距的参考值组成的信息向量，y(k)为k时刻系统输出数据，即当前时刻的桨距角，上标T表示向量的转置，θt为浆距子系统T
参数真值向量，θt＝[a1,a2,b1,b2] ；e(k)为k时刻浆距子系统受到的噪声，|e(k)|≤δ(k)，δ(k)为浆距子系统在k时刻受到噪声绝对值的最大值。
[0026] 可选的，所述S2包括：
[0027] 由式(4)可知θt＝[a1,a2,b1,b2]T为系统参数真值向量，
[0028] 定义k时刻的测量带信息为S(pk,ck)；
[0029]
[0030] 其中，ck为k时刻测量带信息的中心点， pk为k时刻测量带信息的方向向量， θs为测量带信息估计的包裹参数真值的范围。
[0031] 可选的，所述S3包括：
[0032] 根据参数真值向量，确定初始化凸空间 以包裹参数真值向量，令初始时刻k＝1的凸空间 其中，凸空间定义为：
[0033]
[0034] 其中，T为凸空间的生成矩阵，θc为凸空间的中心，θ为中间向量，‖·‖∞表示无穷范数，θ表示凸空间包裹范围内所有的参数可行集。
[0035] 可选的，所述S4中，
[0036] k≠1时，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障，若交集为空，系统发生故障，即标量 或者 若交集不为空，系统没有发生故障，即标量 且
[0037]
[0038] 其中，p0，c0为中间变量，p0＝pk，c0＝ck，θc,k‑1为k‑1时刻凸空间的中心，tk‑1,h为k‑1时刻凸空间生成矩阵的第h列，h为整数，h＝1,...,n。
[0039] 可选的，所述S4中，
[0040] 当k时刻的测量带信息与k‑1时刻的凸空间的交集不为空时，系统没有发生故障，即式(7)中的标量 且 则求取k时刻包裹参数真值向量的凸空间 和正交凸空间
[0041] 由式(7)中的标量 和 求得中间变量 和 即
[0042]
[0043] 得到缩减后的测量带信息 为缩减后的测量带信息的方向向量， 为缩减后的测量带信息的中心：
[0044]
[0045]
[0046] 引入中间变量 得到缩减后的凸空间为缩减后的凸空间的生成矩阵，为缩减后的凸空间的生成矩阵的第g列，g为
整数，g＝1,...,n； 为缩减后的凸空间的中心：
[0047]
[0048]
[0049] 则 可有以下公式获得
[0050]
[0051] 其中，i*， m为中间变量， 为中间*
凸空间， 为中间凸空间生成矩阵T 的第w列向量，w为整数，w＝1,...,n； 为中间凸空间的中心：
[0052]
[0053]
[0054] 由式(7)～(15)求出包裹参数真值向量的凸空间 后，利用以下式(16)～(20)可求出正交凸空间
[0055] 定义n维凸空间生成矩阵的每一列为tj，其中j＝1,...,n，计算出n维凸空间的顶点Vl：
[0056]
[0057] 其中l为顶点的个数，l＝1,2,…,2n，αl,j为中间变量，αl,j∈{‑1,1}；Rn表示n维实数集；
[0058] 定义k‑1时刻n维凸空间的顶点矩阵为 凸空间有多少个顶点，顶点矩阵就有多少列，Vk‑1(p,l)表示k‑1时刻Vk‑1矩阵第l列的第p个元素，p为整数，p＝1,…,n；
[0059] 由式(16)求出k‑1时刻n维凸空间中每个参数的边界最大值 和边界最小值[0060]
[0061] 将k时刻与k‑1时刻的凸空间的参数边界最大值对比取最小，参数边界最小值对比取最大，得到
[0062]
[0063] 其中，v为整数，v＝1,…,n， 为中间变量；
[0064] 由式(18)定义k时刻的正交凸空间为
[0065]
[0066] 其中，
[0067]
[0068] 其中，θc,k为k时刻正交凸空间的中心，Tk为k时刻正交凸空间的生成矩阵，‖θ‖∞为中间变量；
[0069] 当系统在无故障情况下，重复式(7)～(20)可求出每个时刻包裹真值向量的凸空间和单调收敛的正交凸空间。
[0070] 可选的，所述S4中，
[0071] 若系统发生故障，则将k‑1时刻的正交凸空间 扩容获得n个正交凸空间测试集 其中i＝1,…,n；
[0072]
[0073] 其中
[0074]
[0075] (1)当f≠i时，
[0076]
[0077] (2)当f＝i时，
[0078]
[0079] 其中，f为整数，f＝1,...,n； 为第i个正交凸空间测试集里的第f个参数边界最大值， 为第i个正交凸空间测试集里的第f个参数边界最小值， 为正交凸空间第f个参数扩容的最大幅值；
[0080] 给定时间长度L，在k时刻检测到发生故障后的L时刻内无新故障发生，获得测量带信息S(k),…,S(k+L)；利用式(7)～(15)将扩容得到的n个正交凸空间测试集分别与L个测量带信息计算获得k+L时的n个凸空间测试集 q为整数，q＝1,...,n。
[0081] 可选的，所述根据凸空间测试集的状态隔离出系统具体的故障参数，对具体故障参数进行扩容得到包裹故障真值向量的正交凸空间，完成故障识别，包括：
[0082] 若 为空集，说明参数真值向量的第q个分量发生故障；
[0083] 若 不为空集，说明参数真值向量的第q个分量没有发生故障；
[0084] 当在k时刻隔离出系统具体的故障参数后，扩容k‑1时刻的正交凸空间，作为k时刻包裹故障真值向量的正交凸空间 即仅在故障参数方向上进行扩容，假设参数真值向量的第q个分量发生故障，包裹故障真值向量的正交凸空间为：
[0085]
[0086] 其中
[0087]
[0088] (1)当d＝q时，
[0089]
[0090] (2)当d≠q时，
[0091]
[0092] 其中，d为整数，d＝1,...,n； 为正交凸空间里的第d个参数边界最大值，为正交凸空间里的第d个参数边界最小值； 为正交凸空间第d个参数扩容的最大幅值；
[0093] 扩容后得到的正交凸空间 是包裹故障参数真值向量的，将正交凸空间作为当前时刻的凸空间，利用上式(7)～(15)求出包裹参数真值向量的凸空间
以及式(16)～(20)求出正交凸空间，不断迭代至求出的正交凸空间在预定空间阈值范围内，完成故障识别。
[0094] 可选的，所述方法应用于浆距子系统的故障检测时，取ζ＝0.6，ωn＝11.11rad/s，采样时间Ts＝0.01s。
[0095] 可选的，所述工业系统包括浆距子系统、直流电动机系统、弹簧阻尼系统。
[0096] 本发明有益效果是：
[0097] 本发明通过获取工业系统的线性模型；根据系统的线性模型，获取系统的k时刻的测量带信息和参数真值向量；根据系统的参数真值向量，确定包裹参数真值向量的初始化凸空间，将此初始化凸空间作为初始时刻k＝1的凸空间；k≠1时根据系统的k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障；若系统没有发生故障，则求取k时刻的凸空间和正交凸空间，若系统发生故障，扩容k‑1时刻的正交凸空间与测量带信息计算得到凸空间测试集；根据凸空间测试集的状态隔离出系统具体的故障参数，对具体故障参数进行扩容得到包裹故障真值向量的正交凸空间，完成故障识别。本发明提供的故障诊断方法，在检测到故障发生时，及时进行故障隔离，立马判断出具体的故障参数，隔离出故障参数后仅仅扩容故障参数的滤波估计区间，而不扩容其他非故障参数的滤波估计区间，从而进行下一步的故障识别，降低传统凸空间滤波方法的保守性，提升了系统故障诊断精度。
附图说明
[0098] 为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0099] 图1是本申请一个实施例中公开的基于凸空间滤波的工业系统故障诊断方法的流程图。
[0100] 图2是本申请一个实施例中公开的在外加浆距子系统的故障信号后，采用本发明方法和现有基于椭球滤波分别进行故障诊断得到的浆距子系统参数真值向量中的分量θ1的对比图。
[0101] 图3是本申请一个实施例中公开的在外加浆距子系统的故障信号后，采用本发明方法和现有基于椭球滤波分别进行故障诊断得到的浆距子系统参数真值向量中的分量θ2的对比图。
[0102] 图4是本申请一个实施例中公开的在外加浆距子系统的故障信号后，采用本发明方法和现有基于椭球滤波分别进行故障诊断得到的浆距子系统参数真值向量中的分量θ3的对比图。
[0103] 图5是本申请一个实施例中公开的在外加浆距子系统的故障信号后，采用本发明方法和现有基于椭球滤波分别进行故障诊断得到的浆距子系统参数真值向量中的分量θ4的对比图。
具体实施方式
[0104] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0105] 实施例一：
[0106] 本实施例提供一种基于凸空间滤波的工业系统故障诊断方法，所述方法包括：
[0107] S1：确定工业系统的线性离散模型；
[0108] S2：根据S1确定的线性离散模型，获取工业系统的参数真值向量和k时刻的测量带信息；
[0109] S3：根据S2获取的工业系统的参数真值向量，确定包裹参数真值向量的初始化凸空间，将此初始化凸空间作为初始时刻k＝1的凸空间；
[0110] S4：k≠1时刻，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障；
[0111] 若系统没有发生故障，则求取k时刻的凸空间和正交凸空间，所述正交凸空间为k‑1时刻的凸空间和k时刻的凸空间的交集；
[0112] 若系统发生故障，则扩容k‑1时刻的正交凸空间，进而计算得到凸空间测试集；
[0113] 根据凸空间测试集的状态隔离出系统的具体故障参数；
[0114] 对具体故障参数进行扩容得到包裹故障真值向量的k时刻的正交凸空间，将k时刻的正交凸空间作为k时刻的凸空间，重复S4，完成故障辨识。
[0115] 实施例二：
[0116] 本实施例提供一种基于凸空间滤波的工业系统故障诊断方法，以浆距子系统为例进行说明如下，请参见图1，所述方法包括：
[0117] 步骤101，获取工业系统中的浆距子系统的线性离散模型。
[0118] 获取浆距子系统的状态空间模型为：
[0119]
[0120] 其中β和βa分别为桨距角和角速度大小， 和 分别为桨距角的一阶导数和角速度的一阶导数，βr为桨距的参考值，ζ和ωn分别为阻尼系数和浆距子系统的固有频率；
[0121] 将式(1)的状态空间模型转化为二阶连续系统模型：
[0122]
[0123] 其中y为输出，表示桨距角；u为输入，表示桨距的参考值，即y＝β，u＝βr[0124] 利用采样时刻Ts，将式(2)离散化为：
[0125] A(z)y(k)＝B(z)u(k)+e(k)(3)
[0126] 其中A(z)和B(z)为z变换多项式，A(z)＝1+a1z‑1+a2z‑2,B(z)＝b1z‑1+b2z‑2；u(k)为k时刻系统输入数据，即桨距的参考值；则由式(3)获取浆距子系统的线性离散模型为：
[0127]
[0128] 其中k表示时刻，k＝1,…,N， 为可观测的数据向量，表示浆距子系统中由桨距角历史值和桨距的参考值组成的信息向量，y(k)为k时刻系统输出数据，即当前时刻的桨距角，上标T表示向量的转置，θt为浆距子系统T
参数真值向量，θt＝[a1,a2,b1,b2] ；e(k)为k时刻浆距子系统受到的噪声，|e(j)|≤δ(k)，δ(k)为浆距子系统在k时刻受到噪声绝对值的最大值。
[0129] 步骤102，根据系统的线性离散模型，获取工业系统的参数真值向量和k时刻的测量带信息。
[0130] 由式(4)可知θt＝[a1,a2,b1,b2]T为系统参数真值向量。
[0131] 定义k时刻的测量带信息为S(pk,ck)：
[0132]
[0133] 其中，ck为k时刻测量带信息的中心点， pk为k时刻测量带信息的方向向量， θs为测量带信息估计的包裹参数真值的范围。
[0134] 步骤103，根据参数真值向量，确定包裹参数真值向量的初始化凸空间，将此初始化凸空间作为初始时刻k＝1的凸空间：
[0135] 根据参数真值向量，确定初始化凸空间 以包裹参数真值向量，令初始时刻k＝1的凸空间 其中，凸空间定义为：
[0136]
[0137] 其中，T为凸空间的生成矩阵，θc为凸空间的中心，θ为中间向量，‖·‖∞表示无穷范数，θ表示凸空间包裹范围内所有的参数可行集。
[0138] 步骤104，k≠1时，根据系统k时刻的测量带信息与k‑1时刻的凸空间的交集判断有无发生故障：
[0139] 若交集为空，系统发生故障，即标量 或者 若交集不为空，系统没有发生故障，即标量 且
[0140]
[0141] 其中，p0，c0为中间变量，p0＝pk，c0＝ck，θc,k‑1为k‑1时刻凸空间的中心，tk‑1,h为k‑1时刻凸空间生成矩阵的第h列，h为整数，h＝1,...,n。
[0142] 步骤105，若系统没有发生故障，则求取k时刻的凸空间和正交凸空间。
[0143] 当k时刻的测量带信息与k‑1时刻的凸空间的交集不为空时，系统没有发生故障，即式(7)中的标量 且 则求取k时刻包裹参数真值向量的凸空间 和正交凸空间
[0144] 由式(7)中的标量 和 求得中间变量 和 即
[0145]
[0146] 得到缩减后的测量带信息 为缩减后的测量带信息的方向向量， 为缩减后的测量带信息的中心：
[0147]
[0148]
[0149] 引入中间变量 得到缩减后的凸空间为缩减后的凸空间的生成矩阵，为缩减后的凸空间的生成矩阵的第g列，g为
整数，g＝1,...,n； 为缩减后的凸空间的中心：
[0150]
[0151]
[0152] 则 可有以下公式获得
[0153]
[0154] 其中，i*， m为中间变量， 为中间*
凸空间， 为中间凸空间生成矩阵T的第w列向量，w为整数，w＝1,...,n； 为中间凸空间的中心：
[0155]
[0156]
[0157] 由式(7)～(15)求出包裹参数真值向量的凸空间 后，利用以下式(16)～(20)可求出正交凸空间
[0158] 定义n维凸空间生成矩阵的每一列为tj，其中j＝1,…,n，计算出n维凸空间的顶点Vl：
[0159]
[0160] 其中l为顶点的个数，l＝1,2,…,2n，αl,j为中间变量，αl,j∈{‑1,1}；Rn表示n维实数集；
[0161] 定义k‑1时刻n维凸空间的顶点矩阵为 凸空间有多少个顶点，顶点矩阵就有多少列，Vk‑1(p,l)表示k‑1时刻Vk‑1矩阵第l列的第p个元素，p为整数，p＝1,…,n；
[0162] 由式(16)求出k‑1时刻n维凸空间中每个参数的边界最大值 和边界最小值[0163]
[0164] 将k时刻与k‑1时刻的凸空间的参数边界最大值对比取最小，参数边界最小值对比取最大，得到
[0165]
[0166] 其中，v为整数，v＝1,…,n， 为中间变量；
[0167] 由式(18)定义k时刻的正交凸空间为
[0168]
[0169] 其中，
[0170]
[0171] 其中，θc,k为k时刻正交凸空间的中心，Tk为k时刻正交凸空间的生成矩阵，‖θ‖∞为中间变量。
[0172] 当系统在无故障情况下，重复式(7)～(20)可求出每个时刻包裹真值向量的凸空间和单调收敛的正交凸空间。
[0173] 步骤106，若系统发生故障，则扩容k‑1时刻的正交凸空间，进而计算得到凸空间测试集。
[0174] 若系统发生故障，则将k‑1时刻的正交凸空间 扩容获得n个正交凸空间测试集 其中i＝1,…,n；
[0175]
[0176] 其中
[0177]
[0178] (1)当f≠i时，
[0179]
[0180] (2)当f＝i时，
[0181]
[0182] 其中，f为整数，f＝1,...,n； 为第i个正交凸空间测试集里的第f个参数边界最大值， 为第i个正交凸空间测试集里的第f个参数边界最小值， 为正交凸空间第f个参数扩容的最大幅值；
[0183] 给定时间长度L，在k时刻检测到发生故障后的L时刻内无新故障发生，获得测量带信息S(k),…,S(k+L)；利用式(7)～(15)将扩容得到的n个正交凸空间测试集分别与L个测量带信息计算获得k+L时的n个凸空间测试集 q为整数，q＝1,...,n。
[0184] 步骤107，根据凸空间测试集的状态隔离出系统具体的故障参数，对具体故障参数进行扩容得到包裹故障真值向量的正交凸空间，完成故障识别，包括：
[0185] 若 为空集，说明参数真值向量的第q个分量发生故障；
[0186] 若 不为空集，说明参数真值向量的第q个分量没有发生故障；
[0187] 当在k时刻隔离出系统系统具体的故障参数后，扩容k‑1时刻的正交凸空间，作为k时刻包裹故障真值向量的正交凸空间 即仅在故障参数方向上进行扩容，假设参数真值向量的第q个分量发生故障，包裹故障真值向量的正交凸空间为：
[0188]
[0189] 其中
[0190]
[0191] (1)当d＝q时，
[0192]
[0193] (2)当d≠q时，
[0194]
[0195] 其中，d为整数，d＝1,...,n...，n。 为正交凸空间里的第d个参数边界最大值， 为正交凸空间里的第d个参数边界最小值； 为正交凸空间第d个参数扩容的最大幅值；
[0196] 扩容后得到的正交凸空间 是包裹故障参数真值向量的，将正交凸空间作为当前时刻的凸空间，利用上式(7)～(15)求出包裹参数真值向量的凸空间
以及式(16)～(20)求出正交凸空间，不断迭代至求出的正交凸空间在预定空间阈值范围内，完成故障识别。
[0197] 为验证本申请提出的基于凸空间滤波的工业系统故障诊断方法的有效性，本申请进行仿真实验，实验中，取ζ＝0.6，ωn＝11.11rad/s，采样时间Ts＝0.01s，则式(3)中的z变换多项式为：
[0198] A(z)＝1‑1.864z‑1+0.8752z‑2,B(z)＝0.0059z‑1+0.0056z‑2
[0199] 则式(4)中的可观测的数据向量为：
[0200]
[0201] 参数真值向量为：
[0202] θt＝[‑1.864,0.8752,0.0059,0.0056]T
[0203] 待辨识的参数个数为n＝4。
[0204] 输入信号|u(k)|≤20，噪声幅值|e(k)|≤0.2，时间长度L＝5，扩容的幅值[0205] 针对于浆距子系统，其常见的三种故障如下：
[0206] 第1种故障是执行器中的压降作用导致压力降为原来的50％，由于这种故障的发生，导致参数真值向量的θ1发生变化，即由θ1＝‑1.864变为θ1＝‑1.520；
[0207] 第2种故障是泵磨损，泵磨损是一个潜在的故障，发生故障时将导致桨距闭环系统发生参数变化，这种故障的发生会导致参数真值向量的θ2,θ3发生变化，即由正常值变成θ2＝1.020,θ3＝0.4038；
[0208] 第3种故障是执行机构内的空气含量增加，在正常情况下，液压油中的空气为7％，由于这种故障的发生，这种故障的发生会导致参数真值向量的θ3,θ4发生变化，即由正常值变成θ3＝0.2038,θ4＝0.3028。
[0209] 在本仿真中分别在k＝501,1001,1501三个时刻加入以上三种故障，故障诊断要达到的空间阈值由于不同现场和技术人员的要求而不同，在本实验中，了解到发生故障后经过500次的迭代后可以得到一个满足要求很小的空间阈值，因此本实验在检测到故障发生后，设定进行迭代500次得到故障辨识的结果，在针对其他工业系统故障检测时，可由技术人员的先验知识确定具体的空间阈值。
[0210] 根据系统的参数真值向量，设初始化凸空间 包裹参数真值向量，令k＝1时刻的凸空间
[0211] 在预定的时间范围内，循环执行步骤104到107，得到故障诊断结果和滤波估计区间。
[0212] 同时与现有的采用椭球滤波进行故障诊断的方法(可参见Reppa V,Tzes A.Fault diagnosis based on set membership identification using outputerror models[J].International Journal of Adaptive Control and Signal Processing,2016,30(2):224255.)得到的滤波估计区间和故障诊断结果进行对比，对比结果如图2～5中所示。
[0213] 由图2～5中可以看出，在系统未发生故障的时候，即k<501时，与椭球空间滤波方法对比，本发明的方法的滤波估计曲线下降的更快，收敛的速度也更快，说明本发明的方法滤波估计区间更小，保守性更低，精度更高。
[0214] 当系统发生故障时，本发明的方法在检测到故障发生时，及时进行故障隔离，立马判断出具体的故障参数，隔离出故障参数后仅仅扩容故障参数的滤波估计区间，而不扩容其他非故障参数的滤波估计区间，从而进行下一步的故障辨识，例如由图2～5中可以看出在k＝501时刻，θ1扩容但θ2,θ3,θ4不扩容，说明检测出故障发生并及时隔离出故障参数θ1；同理在k＝1001时刻，检测出故障发生并及时隔离出故障参数θ2,θ3，在k＝1501时刻，检测出故障发生并及时隔离出故障参数θ3,θ4。但现有的椭球空间滤波在检测到故障发生时，并不是及时的进行故障隔离，而是扩容所有参数的滤波估计区间，通过一段时间后的故障识别曲线才能判断出具体发生故障的参数。同时由于椭球空间滤波扩容非故障参数的滤波估计区间，所以其在故障诊断时影响了非故障参数的滤波估计，例如在k＝501时刻，θ2,θ3,θ4未发生故障，从图3～5中可以看出，本发明的方法的滤波估计曲线沿着非故障时刻的曲线继续收敛，而椭球空间滤波的滤波估计曲线突然扩容，导致后面的滤波估计曲线收敛性远小于本发明的方法，同时对于故障参数θ1，从图1可以看出，在同样扩容的情况下，本发明的方法的滤波估计曲线收敛速度更快，包裹真值的上下界区间更小，保守性更低，精度更高。
[0215] 综上，针对于现有的椭球空间滤波方法，本发明的方法在快速的检测到系统故障时，能够及时的对故障参数进行隔离，故障诊断的效果更好，并且无论是否在故障状态，本发明的方法的滤波估计曲线收敛性更快，滤波估计区间都更小，保守性更低，精度更高。
[0216] 本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或硬盘等。
[0217] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
