CN110826021B - 一种非线性工业过程鲁棒辨识和输出估计方法 - Google Patents

一种非线性工业过程鲁棒辨识和输出估计方法 Download PDF

Info

Publication number
CN110826021B
CN110826021B CN201911053823.XA CN201911053823A CN110826021B CN 110826021 B CN110826021 B CN 110826021B CN 201911053823 A CN201911053823 A CN 201911053823A CN 110826021 B CN110826021 B CN 110826021B
Authority
CN
China
Prior art keywords
model
distribution
parameter
output
local
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911053823.XA
Other languages
English (en)
Other versions
CN110826021A (zh
Inventor
杨宪强
刘新鹏
高会军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201911053823.XA priority Critical patent/CN110826021B/zh
Publication of CN110826021A publication Critical patent/CN110826021A/zh
Application granted granted Critical
Publication of CN110826021B publication Critical patent/CN110826021B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/04—Manufacturing
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Business, Economics & Management (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Economics (AREA)
  • Tourism & Hospitality (AREA)
  • Strategic Management (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Primary Health Care (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Marketing (AREA)
  • Human Resources & Organizations (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • General Business, Economics & Management (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Health & Medical Sciences (AREA)
  • Algebra (AREA)
  • Manufacturing & Machinery (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Abstract

一种非线性工业过程鲁棒辨识和输出估计方法，涉及工业过程建模及模型参数辨识领域，针对现有技术当输出数据中存在异常值时，会导致系统辨识精度降低的问题，包括以下步骤：步骤一：选取系统局部模型，并基于拉普拉斯分布，建立多模型非线性系统的鲁棒概率模型；步骤二：根据变分贝叶斯框架，建立隐变量后验分布和待估计参数的迭代更新公式；步骤三：设定步骤二中所建立的隐变量后验分布和待估计参数迭代更新公式的终止条件，当迭代终止时，记录最终的迭代结果作为估计的最优参数，进而通过局部模型插值得到模型输出值。

Description

一种非线性工业过程鲁棒辨识和输出估计方法
技术领域
本发明涉及工业过程建模及模型参数辨识领域，具体为一种非线性工业过程鲁棒辨识和输出估计方法。
背景技术
实际工业过程中，出于资源、成本等方面因素的考虑，往往会在多个工况进行转移，进而会导致系统存在非线性特性。对这些动态特性进行精确建模，是实现后续状态估计和控制的必要前提。单一的线性模型往往难以在大的工作范围内描述系统的非线性特性。而采取多个线性模型加权组合的方式，凭借局部线性结构和全局非线性特性，能较好地反映工业过程中存在的非线性动态，因此得到广泛应用。
由于工业过程不可避免会出现传感器故障及外界干扰等因素，往往会导致收集的过程输出数据中存在异常值。若在辨识算法中为考虑这些异常值的作用，往往会导致系统辨识精度的降低。
发明内容
本发明的目的是：针对现有技术由于工业过程不可避免会出现传感器故障及外界干扰等因素，往往会导致收集的过程输出数据中存在异常值。若在辨识算法中为考虑这些异常值的作用，往往会导致系统辨识精度的降低的问题，提出一种非线性工业过程鲁棒辨识和输出估计方法。
本发明为了解决上述技术问题采取的技术方案是：一种非线性工业过程鲁棒辨识和输出估计方法，包括以下步骤：
步骤一：选取系统局部模型，并基于拉普拉斯分布，建立多模型非线性系统的鲁棒概率模型；
步骤二：根据变分贝叶斯框架，建立隐变量后验分布和待估计参数的迭代更新公式；
步骤三：设定步骤二中所建立的隐变量后验分布和待估计参数迭代更新公式的终止条件，当迭代终止时，记录最终的迭代结果作为估计的最优参数，进而通过局部模型插值得到模型输出值。
进一步的，所述步骤一的具体步骤为：
步骤一一：选取自回归各态历经模型，即ARX模型，作为多模型系统的局部模型，该局部模型的回归量中包含收集到的工业生产过程输入和输出数据；
步骤一二：基于拉普拉斯分布，建立系统的鲁棒概率模型
Figure GDA0002798188100000021
进而得到系统输出的条件概率分布为拉普拉斯分布
Figure GDA0002798188100000022
其中，zn＝[zn1,...,znK]为二值隐变量，当znk＝1时，表示第n个采样点处的局部模型身份为k；λ＝[λ1,...,λK]为拉普拉斯分布的尺度参数，拉普拉斯分布满足
Figure GDA0002798188100000023
即
Figure GDA0002798188100000024
分解为
Figure GDA0002798188100000025
和
Figure GDA0002798188100000026
其中，vn为拉普拉斯分布的隐变量；
步骤一三：当辨识得到K个局部ARX模型的参数值，非线性系统的全局模型由K个局部ARX模型插值得到，系统输出的估计值表示为
Figure GDA0002798188100000027
其中，
Figure GDA0002798188100000028
分别为第k个子模型的参数估计值和输出估计值，加权系数αnk表示为
Figure GDA0002798188100000029
ωnk选取为高斯核函数的形式，即
Figure GDA00027981881000000210
其中Tk为预先设定的工作点，Hn为可测量的调度变量的值，ok为第k个子模型的有效宽度，表示局部模型身份的二值变量zn表示为
Figure GDA00027981881000000211
步骤一四：选取系统模型参数的先验分布，系统模型参数假定服从高斯分布，
Figure GDA0002798188100000031
其中超参数δk服从Gamma分布
Figure GDA0002798188100000032
进一步的，所述步骤一一中ARX模型可表示为：
Figure GDA0002798188100000033
其中，
Figure GDA0002798188100000034
为模型回归量，yn为收集到的非线性过程的输出值，un为输入值，
Figure GDA0002798188100000035
为第In个局部模型的模型参数，In∈{1,2,...,K}表示局部模型的身份，en为随机噪声，n＝1,2,...,N表示过程的采样点。
进一步的，所述步骤二的详细步骤为：
步骤二一：根据全概率公式，得到系统鲁棒概率模型的联合概率分布；
步骤二二：引入概率密度函数
Figure GDA0002798188100000036
其中hi分别对应Z,Θ,δ,v，得到
输出变量的对数似然函数
Figure GDA0002798188100000037
其中，
Figure GDA0002798188100000038
为变分下界，KL(·)为变分分布q(h)与
Figure GDA0002798188100000039
之间的KL散度，当
Figure GDA00027981881000000310
时，KL散度为零；所述变分贝叶斯框架包含两大步骤：
一、VB E步:
Figure GDA00027981881000000311
步中每个隐变量的后验分布通过如下公式更新：
Figure GDA00027981881000000312
二、VB M步:
Figure GDA00027981881000000313
该步中未知参数通过优化算法得到；
步骤二三：VB E步,更新隐变量的后验概率如下：
(1)q(Z)服从多项分布:
Figure GDA00027981881000000314
其中
Figure GDA0002798188100000041
Figure GDA0002798188100000042
<a>b表示a关于b的数学期望；
(2)q(θk)服从正态分布
Figure GDA0002798188100000043
其中
Figure GDA0002798188100000044
Figure GDA0002798188100000045
(3)q(δk)服从Gamma分布
Figure GDA0002798188100000046
其中
Figure GDA0002798188100000047
Figure GDA0002798188100000048
p为局部模型的阶次，tr(·)表示矩阵的迹；
(4)q(vn)服从广义逆高斯分布
q(vn)＝GIG(vn|p,a,b)
其中
Figure GDA0002798188100000049
Figure GDA00027981881000000410
Figure GDA00027981881000000411
根据以上分布，得各变量的期望值如下
Figure GDA0002798188100000051
步骤二四：VB M步,根据
Figure GDA0002798188100000052
更新未知参数如下：
Figure GDA0002798188100000053
Figure GDA0002798188100000054
其中最优的参数可以通过Matlab的“fmincon”优化函数求解得到。
进一步的，所述步骤二一中系统鲁棒概率模型的联合概率分布为：
Figure GDA0002798188100000055
其中，h＝{Z,Θ,δ,v}视为隐变量，
Figure GDA0002798188100000056
视为待估计参数，
Z＝{zn}n＝1,...,N,Θ＝{θk}k＝1,...,K,δ＝{δk}k＝1,...,K,v＝{vn}n＝1,...,N,λ＝{λk}k＝1,...,K,o＝{ok}k＝1,...,K分别为二值模型身份矩阵、局部模型参数矩阵、Gamma分布的超参数、拉普拉斯分布的尺度参数和局部模型的有效宽度。
进一步的，所述步骤三中所建立的隐变量后验分布和待估计参数的迭代更新的终止条件如下：
Figure GDA0002798188100000057
其中，ε为迭代终止阈值，当迭代终止时，记录最终的迭代结果，作为估计的最优参数
Figure GDA0002798188100000061
等，并根据
Figure GDA0002798188100000062
估计得到模型输出值。
进一步的，所述优化算法为拉格朗日乘子法。
本发明的有益效果是：针对工业过程中存在的异常值问题，基于拉普拉斯分布，对非线性系统进行了鲁棒建模；并基于变分贝叶斯框架，将系统参数估计过程中后验概率密度函数涉及的高维积分(求和)难以求解问题，转化为近似变分分布的估计问题，进而提高了参数估计效率，同时保证了参数估计和输出估计的精度，对于非线性过程的鲁棒辨识理论和实际工业过程应用具有重要的意义。
附图说明
图1为本发明的流程图。
图2为信噪比SNR＝25dB，输出异常值比例5％的情况下50次蒙特卡洛仿真参数估计误差棒图。
图3为信噪比SNR＝25dB，输出异常值比例10％的情况下50次蒙特卡洛仿真参数估计误差棒图
图4为不同信噪比和输出异常值比例下的参数偏差范数(BN)图。
图5为进行鲁棒建模与未进行鲁棒建模VB算法的输出估计对比图。
图6为本发明算法的参数相对估计误差(RPEE)随迭代次数变化曲线图。
具体实施方式
具体实施方式一：参照图1具体说明本实施方式，本实施方式所述的一种非线性工业过程鲁棒辨识和输出估计方法，包括以下步骤：
步骤一：选取系统局部模型，并基于拉普拉斯分布，建立多模型非线性系统的鲁棒概率模型；
步骤二：根据变分贝叶斯框架，建立隐变量后验分布和待估计参数的迭代更新公式；
步骤三：设定步骤二中所建立的隐变量后验分布和待估计参数迭代更新公式的终止条件，当迭代终止时，记录最终的迭代结果作为估计的最优参数，进而通过局部模型插值得到模型输出值。
本发明的详细步骤为：
第一步：基于拉普拉斯分布，建立多模型非线性系统的鲁棒概率模型，具体过程为：
步骤1.1：由于ARX模型结构能较好地刻画系统的输入输出特性，因此选取作为多模型系统的局部模型，可表示为：
Figure GDA0002798188100000071
其中yn为系统输出，
Figure GDA0002798188100000072
为模型回归量，
Figure GDA0002798188100000073
为第In个局部模型的模型参数，In∈{1,2,...,K}表示局部模型的身份，en为随机噪声，n＝1,2,...,N表示过程的采样点。
步骤1.2：传统方法假设随机噪声en服从零均值的高斯分布，但高斯分布模型对实际工业过程中的异常点比较敏感，进而会导致辨识算法的估计效果较差。因此，本发明基于更重尾的拉普拉斯分布，建立系统的鲁棒概率模型，即
Figure GDA0002798188100000074
进而可得到系统输出的条件概率分布仍为拉普拉斯分布
Figure GDA0002798188100000075
其中zn＝[zn1,...,znK]为二值隐变量，当znk＝1时，表示第n个采样点处的局部模型身份为k；λ＝[λ1,...,λK]为拉普拉斯分布的尺度参数。拉普拉斯分布满足
Figure GDA0002798188100000076
故(3)式可分解为
Figure GDA0002798188100000077
Figure GDA0002798188100000078
其中vn为拉普拉斯分布的隐变量。
步骤1.3：当辨识得到K个局部ARX模型的参数值后，非线性系统的全局模型由K个局部ARX模型插值得到，进而系统输出的估计值可表示为
Figure GDA0002798188100000079
其中
Figure GDA00027981881000000710
分别为第k个子模型的参数估计值和输出估计值，加权系数αnk可表示为
Figure GDA0002798188100000081
ωnk选取为高斯核函数的形式，即
Figure GDA0002798188100000082
其中Tk为预先设定的工作点，Hn为可测量的调度变量的值，ok为第k个子模型的有效宽度(待定参数)。
因此，表示局部模型身份的二值变量zn可表示为如下多项分布的形式
Figure GDA0002798188100000083
步骤1.4：选取系统模型参数的先验分布。为形成共轭先验，系统模型参数假定服从高斯分布
Figure GDA0002798188100000084
其中超参数δk服从Gamma分布
Figure GDA0002798188100000085
至此，基于拉普拉斯分布，建立了非线性系统的鲁棒概率模型。本发明的主要目的是，根据观测数据集W＝{u,y,H,T}，基于变分贝叶斯框架，鲁棒辨识得到非线性系统的局部模型参数，并估计得到系统的真实输出值。
第二步：根据变分贝叶斯框架，建立隐变量后验分布和待估计参数的迭代更新公式，具体过程为：
步骤2.1：根据全概率公式，得到系统鲁棒概率模型的联合概率分布为
Figure GDA0002798188100000086
其中h＝{Z,Θ,δ,v}视为隐变量，
Figure GDA0002798188100000087
视为待估计参数。Z＝{zn}n＝1,...,N,Θ＝{θk}k＝1,...,K,δ＝{δk}k＝1,...,K,v＝{vn}n＝1,...,N,λ＝{λk}k＝1,...,K,o＝{ok}k＝1,...,K分别为二值模型身份矩阵、局部模型参数矩阵、Gamma分布的超参数、拉普拉斯分布的尺度参数和局部模型的有效宽度。
步骤2.2：引入概率密度函数
Figure GDA0002798188100000091
其中hi分别对应Z,Θ,δ,v。则输出变量的对数似然函数可表示为
Figure GDA0002798188100000092
其中
Figure GDA0002798188100000093
为变分下界，KL(·)为变分分布q(h)与
Figure GDA0002798188100000094
之间的KL散度，当
Figure GDA0002798188100000095
时，KL散度为零。故求解隐变量的后验概率分布问题等价于极大化变分下界问题，从而可避免直接求解后验概率时所需计算的高维积分(求和)。因此，变分贝叶斯框架包含两大步骤：
(1)VB E步:
Figure GDA0002798188100000096
该步中每个隐变量的后验分布可通过如下公式更新：
Figure GDA0002798188100000097
(2)VB M步:
Figure GDA0002798188100000098
该步中未知参数可通过已有的优化算法，如拉格朗日乘子法计算得到。
步骤2.3：VB E步,更新隐变量的后验概率如下：
(1)q(Z)服从多项分布:
Figure GDA0002798188100000099
其中
Figure GDA00027981881000000910
Figure GDA00027981881000000911
<a>b表示a关于b的数学期望。
(2)q(θk)服从正态分布
Figure GDA0002798188100000101
其中
Figure GDA0002798188100000102
(3)q(δk)服从Gamma分布
Figure GDA0002798188100000103
其中
Figure GDA0002798188100000104
p为局部模型的阶次，tr(·)表示矩阵的迹。
(4)q(vn)服从广义逆高斯分布
q(vn)＝GIG(vn|p,a,b) (23)
其中
Figure GDA0002798188100000105
根据以上分布，可得各变量的期望值如下
Figure GDA0002798188100000111
步骤2.4：VB M步,根据(15),可更新未知参数如下：
Figure GDA0002798188100000112
其中最优的参数可以通过Matlab的“fmincon”优化函数求解得到。
第三步：第二步所建立的隐变量后验分布和待估计参数的迭代更新的终止条件如下：
Figure GDA0002798188100000113
其中，ε为迭代终止阈值。当迭代终止时，记录最终的迭代结果，作为估计的最优参数
Figure GDA0002798188100000114
等，并根据(6)式估计得到模型输出值。
实施例：
(1)选取一个一阶过程，其传递函数如下：
Figure GDA0002798188100000115
其中K(w)＝w2+0.6为系统增益，τ(w)＝0.5w3+3为系统时间常数，调度变量w的取值范围为w∈[1,4]。由于在该工作范围内，系统增益和时间常数可变化10倍以上，因此单一的线性模型难以描述该系统的动态特性。现采取多个局部ARX模型加权组合的方式，进行该非线性过程的辨识实验。选取三个工作点：w＝1,w＝2.25和w＝4。调度变量按如下方式变化：
1～100s：位于工作点w＝1；
101～400s：由工作点w＝1线性变化至工作点w＝2.25；
401～550s：位于工作点w＝2.25；
551～750s：由工作点w＝1线性变化至工作点w＝4；
751～900s：位于工作点w＝4；
为辨识该非线性过程，输入信号设计为随机二值序列，并在收集到的输出数据上，添加不同程度的噪声和异常值(均匀分布在[-5,5]范围内)，以验证算法的有效性。
(2)仿真结果：
固定随机噪声的峰值信噪比为25dB，异常值比例为5％和10％时，进行50次蒙特卡洛仿真，并将所估计的50个模型参数求取均值和标准差，得到参数估计误差棒图如图2和图3所示，其中圆圈为局部模型参数的真实值，菱形为蒙特卡洛仿真估计参数的平均值，竖直棒条为蒙特卡洛仿真估计参数的标准差(棒条长度越短，估计效果越好)。由图可知，所估计参数与真实参数值较为接近，但由于噪声、异常值等的影响，所估计参数与真实参数值难免会存在偏差。
当选取不同的峰值信噪比(15dB，20dB和25dB)和不同比例的异常值(5％，10％和20％)，并分别进行50次蒙特卡洛仿真时，计算得到的偏差范数(bias norm,BN)‖θ-E(θ*)‖结果图如图4所示。由图4可知，随着输出数据质量的提升(即较高的峰值信噪比和较低的异常值比例)，参数偏差范数变小。
为进一步说明算法的有效性，当峰值信噪比为20dB，输出数据异常值比例为10％时，交叉验证的输出估计图如图5所示，其中实线为真实输出值，虚线为高斯噪声假设下VB算法的输出估计值，点画线为拉普拉斯噪声假设下鲁棒VB算法的输出估计值。可见鲁棒算法能较好地刻画该非线性过程的动态特性。参数相对估计误差
Figure GDA0002798188100000121
随迭代次数的变化曲线如图6所示，可见所估计的参数能在较少的迭代次数后收敛至真值。
从以上仿真结果可知，本发明所公开的一种非线性工业过程鲁棒辨识和输出估计方法，能够在避免后验概率密度高维积分(求和)、提高算法效率的同时，保证辨识算法的精度，具有一定的理论和实际工程价值。
需要注意的是，具体实施方式仅仅是对本发明技术方案的解释和说明，不能以此限定权利保护范围。凡根据本发明权利要求书和说明书所做的仅仅是局部改变的，仍应落入本发明的保护范围内。

Claims (6)

1.一种非线性工业过程鲁棒辨识和输出估计方法，包括以下步骤：
步骤一：选取系统局部模型，并基于拉普拉斯分布，建立多模型非线性系统的鲁棒概率模型；
步骤二：根据变分贝叶斯框架，建立隐变量后验分布和待估计参数的迭代更新公式；
步骤三：设定步骤二中所建立的隐变量后验分布和待估计参数迭代更新公式的终止条件，当迭代终止时，记录最终的迭代结果作为估计的最优参数，进而通过局部模型插值得到模型输出值；其特征在于所述步骤一的具体步骤为：
步骤一一：选取自回归各态历经模型，即ARX模型，作为多模型系统的局部模型，该局部模型的回归量中包含收集到的工业生产过程输入和输出数据；
步骤一二：基于拉普拉斯分布，建立系统的鲁棒概率模型
Figure FDA0002847797750000011
进而得到系统输出的条件概率分布为拉普拉斯分布
Figure FDA0002847797750000012
其中，zn＝[zn1,...,znk,...,znK]为二值隐变量，当znk＝1时，表示第n个采样点处的第k个子模型；λ＝[λ1,...,λK]为拉普拉斯分布的尺度参数，拉普拉斯分布满足
Figure FDA0002847797750000013
即
Figure FDA0002847797750000014
分解为
Figure FDA0002847797750000015
和
Figure FDA0002847797750000016
其中，vn为拉普拉斯分布的隐变量，
Figure FDA0002847797750000017
为模型回归量，yn为收集到的非线性过程的输出值；
步骤一三：当辨识得到K个局部ARX模型的参数值，非线性系统的全局模型由K个局部ARX模型插值得到，系统输出的估计值表示为
Figure FDA0002847797750000018
其中，
Figure FDA0002847797750000019
分别为第k个子模型的参数估计值和输出估计值，加权系数αnk表示为
Figure FDA0002847797750000021
ωnk选取为高斯核函数的形式，即
Figure FDA0002847797750000022
其中Tk为预先设定的工作点，Hn为可测量的调度变量的值，ok为第k个子模型的有效宽度，表示局部模型身份的二值变量zn表示为
Figure FDA0002847797750000023
步骤一四：选取系统模型参数的先验分布，系统模型参数假定服从高斯分布，
Figure FDA0002847797750000024
其中超参数δk服从Gamma分布
Figure FDA0002847797750000025
2.根据权利要求1所述的一种非线性工业过程鲁棒辨识和输出估计方法，其特征在于所述步骤一一中ARX模型可表示为：
Figure FDA0002847797750000026
un为输入值，
Figure FDA0002847797750000027
为第In个局部模型的模型参数，In∈{1,2,...,K}表示局部模型的身份，en为随机噪声，n＝1,2,...,N表示过程的采样点。
3.根据权利要求2所述的一种非线性工业过程鲁棒辨识和输出估计方法，其特征在于所述步骤二的详细步骤为：
步骤二一：根据全概率公式，得到系统鲁棒概率模型的联合概率分布；
步骤二二：引入概率密度函数
Figure FDA0002847797750000028
其中hi分别对应Z,Θ,δ,v，得到
输出变量的对数似然函数
Figure FDA0002847797750000029
其中，W为观测数据，
Figure FDA00028477977500000210
视为待估计参数：o为模型有效宽度、λ为拉普拉斯分布的尺度参数；
Figure FDA00028477977500000211
为变分下界，KL(·)为变分分布q(h)与后验分布
Figure FDA00028477977500000212
之间的KL散度，当
Figure FDA0002847797750000031
时，KL散度为零，Z,Θ,δ,v分别表示二值模型身份矩阵、局部模型参数矩阵、Gamma分布的超参数、拉普拉斯分布的隐变量；所述变分贝叶斯框架包含两大步骤：
一、VBE步:
Figure FDA0002847797750000032
该步中每个隐变量的后验分布通过如下公式更新：
Figure FDA0002847797750000033
二、VBM步:
Figure FDA0002847797750000034
该步中未知参数通过优化算法得到；
步骤二三：VB E步,更新隐变量的后验概率如下：
(1)q(Z)服从多项分布：
Figure FDA0002847797750000035
其中Z表示二值模型身份，
Figure FDA0002847797750000036
Figure FDA0002847797750000037
<a>b表示a关于b的数学期望；
(2)q(θk)服从正态分布
Figure FDA0002847797750000038
其中
Figure FDA0002847797750000039
Figure FDA00028477977500000310
(3)q(δk)服从Gamma分布
Figure FDA00028477977500000311
其中，
Figure FDA0002847797750000041
Figure FDA0002847797750000042
d为局部模型的阶次，tr(·)表示矩阵的迹；
(4)q(vn)服从广义逆高斯分布
q(vn)＝GIG(vn|p,a,b)
其中
Figure FDA0002847797750000043
Figure FDA0002847797750000044
Figure FDA0002847797750000045
根据以上分布，得各变量的期望值如下
Figure FDA0002847797750000046
步骤二四：VBM步,根据
Figure FDA0002847797750000047
更新未知参数如下：
Figure FDA0002847797750000048
Figure FDA0002847797750000049
其中最优的参数可以通过Matlab的fmincon优化函数求解得到。
4.根据权利要求3所述的一种非线性工业过程鲁棒辨识和输出估计方法，其特征在于所述步骤二一中系统鲁棒概率模型的联合概率分布为：
Figure FDA0002847797750000051
其中，h＝{Z,Θ,δ,v}视为隐变量，
Figure FDA0002847797750000052
视为待估计参数，
Z＝{zn}n＝1,...,N,Θ＝{θk}k＝1,...,K,δ＝{δk}k＝1,...,K,v＝{vn}n＝1,...,N,λ＝{λk}k＝1,...,K，o＝{ok}k＝1,...,K分别为二值模型身份矩阵、局部模型参数矩阵、Gamma分布的超参数、拉普拉斯分布的隐变量、尺度参数和局部模型的有效宽度。
5.根据权利要求4所述的一种非线性工业过程鲁棒辨识和输出估计方法，其特征在于所述步骤三中所建立的隐变量后验分布和待估计参数的迭代更新的终止条件如下：
Figure FDA0002847797750000053
其中，ε为迭代终止阈值，当迭代终止时，记录最终的迭代结果，作为估计的最优参数
Figure FDA0002847797750000054
并根据
Figure FDA0002847797750000055
估计得到模型输出值。
6.根据权利要求3所述的一种非线性工业过程鲁棒辨识和输出估计方法，其特征在于所述优化算法为拉格朗日乘子法。
CN201911053823.XA 2019-10-31 2019-10-31 一种非线性工业过程鲁棒辨识和输出估计方法 Active CN110826021B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911053823.XA CN110826021B (zh) 2019-10-31 2019-10-31 一种非线性工业过程鲁棒辨识和输出估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911053823.XA CN110826021B (zh) 2019-10-31 2019-10-31 一种非线性工业过程鲁棒辨识和输出估计方法

Publications (2)

Publication Number Publication Date
CN110826021A CN110826021A (zh) 2020-02-21
CN110826021B true CN110826021B (zh) 2021-03-12

Family

ID=69551831

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911053823.XA Active CN110826021B (zh) 2019-10-31 2019-10-31 一种非线性工业过程鲁棒辨识和输出估计方法

Country Status (1)

Country Link
CN (1) CN110826021B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111665725B (zh) * 2020-06-24 2021-06-08 哈尔滨工业大学 量测数据随机缺失的机电定位系统鲁棒参数辨识方法
CN111737883B (zh) * 2020-07-30 2021-08-13 哈尔滨工业大学 一种具有输出时滞的非线性双率电路系统鲁棒辨识方法
CN112234612B (zh) * 2020-09-30 2023-08-18 云南电网有限责任公司 一种计及随机扰动幅度的电力系统概率稳定分析方法
CN113111456A (zh) * 2021-04-07 2021-07-13 东南大学溧阳研究院 一种燃气轮机关键运行参数在线区间辨识方法
CN114598611B (zh) * 2022-02-16 2023-04-11 北京科技大学 面向二集值fir系统事件驱动辨识的输入设计方法及系统
CN115186715B (zh) * 2022-07-20 2023-07-28 哈尔滨工业大学 一种基于状态空间模型的机电定位系统贝叶斯辨识方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106021829A (zh) * 2016-07-19 2016-10-12 中南大学 一种基于rbf-arx模型稳定参数估计的非线性系统建模方法
CN108628275A (zh) * 2018-07-04 2018-10-09 杭州电子科技大学 一种化工工业过程模糊约束控制方法
CN108681317A (zh) * 2018-07-11 2018-10-19 杭州电子科技大学 一种化工工业过程鲁棒学习控制方法
CN109033524A (zh) * 2018-06-27 2018-12-18 浙江大学 一种基于鲁棒混合模型的化工过程浓度变量在线估计方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5988419B2 (ja) * 2012-01-11 2016-09-07 インターナショナル・ビジネス・マシーンズ・コーポレーションＩｎｔｅｒｎａｔｉｏｎａｌ Ｂｕｓｉｎｅｓｓ Ｍａｃｈｉｎｅｓ Ｃｏｒｐｏｒａｔｉｏｎ 予測方法、予測システムおよびプログラム

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106021829A (zh) * 2016-07-19 2016-10-12 中南大学 一种基于rbf-arx模型稳定参数估计的非线性系统建模方法
CN109033524A (zh) * 2018-06-27 2018-12-18 浙江大学 一种基于鲁棒混合模型的化工过程浓度变量在线估计方法
CN108628275A (zh) * 2018-07-04 2018-10-09 杭州电子科技大学 一种化工工业过程模糊约束控制方法
CN108681317A (zh) * 2018-07-11 2018-10-19 杭州电子科技大学 一种化工工业过程鲁棒学习控制方法

Also Published As

Publication number Publication date
CN110826021A (zh) 2020-02-21

Similar Documents

Publication Publication Date Title
CN110826021B (zh) 一种非线性工业过程鲁棒辨识和输出估计方法
CN111178385A (zh) 一种鲁棒在线多传感器融合的目标跟踪方法
CN112418051B (zh) 一种用于非线性动态系统非高斯噪声下的状态估计方法
CN109827579B (zh) 一种组合定位中滤波模型实时校正的方法和系统
CN113553755B (zh) 电力系统状态估计方法、装置及设备
Zhang et al. Decoupling control in statistical sense: minimised mutual information algorithm
CN111798494B (zh) 广义相关熵准则下的机动目标鲁棒跟踪方法
CN114567288B (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
CN112737463A (zh) 一种永磁直线同步电机的多目标优化方法及装置
CN110398942B (zh) 一种用于工业生产过程控制的参数辨识方法
Liang et al. Power system sensitivity matrix estimation by multivariable least squares considering mitigating data saturation
Ao et al. Entropy estimation via normalizing flow
CN112287605A (zh) 一种基于图卷积网络加速的潮流校核方法
Telen et al. Uncertainty in optimal experiment design: comparing an online versus offline approaches
CN111259340A (zh) 一种基于logistic回归的饱和负荷预测方法
Mohd Lip et al. Comparative study of smoothing methods and box-jenkins model in forecasting unemployment rate in Malaysia
CN111160464B (zh) 基于多隐层加权动态模型的工业高阶动态过程软测量方法
Kim et al. Derivation and extensions of the linear feedback particle filter based on duality formalisms
Lai et al. Adaptive multinoulli-based Kalman filter with randomly unknown delayed and lost measurements
CN106936628B (zh) 一种计及传感器故障的分数阶网络系统状态估计方法
Singh et al. Cubature and quadrature based continuous-discrete filters for maneuvering target tracking
Zhou et al. High-order moment UKF for Markov jump nonlinear systems
Pillonetto Identification of time-varying systems in reproducing kernel hilbert spaces
Lai et al. Data driven optimal control for stochastic systems with non-Gaussian disturbances
CN110826184B (zh) 一种在时变时滞下narx模型结构和参数的变分贝叶斯辨识方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Yang Xianqiang

Inventor after: Liu Xinpeng

Inventor after: Gao Huijun

Inventor before: Liu Xinpeng

Inventor before: Yang Xianqiang

GR01 Patent grant
GR01 Patent grant