CN116166921A - 一种稀疏自适应的复杂系统模型辨识方法 - Google Patents

一种稀疏自适应的复杂系统模型辨识方法 Download PDF

Info

Publication number
CN116166921A
CN116166921A CN202310268392.9A CN202310268392A CN116166921A CN 116166921 A CN116166921 A CN 116166921A CN 202310268392 A CN202310268392 A CN 202310268392A CN 116166921 A CN116166921 A CN 116166921A
Authority
CN
China
Prior art keywords
candidate
model
function
basis function
mape
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.)
Pending
Application number
CN202310268392.9A
Other languages
English (en)
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.)
Southeast University
Original Assignee
Southeast University
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 Southeast University filed Critical Southeast University
Priority to CN202310268392.9A priority Critical patent/CN116166921A/zh
Publication of CN116166921A publication Critical patent/CN116166921A/zh
Pending legal-status Critical Current

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/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
    • 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
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种稀疏自适应的复杂系统模型辨识方法，包含三个模块：基函数设计模块，设计以三角函数、幂函数等初等函数为主的候选基函数库，对时间节点向量进行变换，实现对系统的精准建模；评价指标求解模块，设计了一种针对各项候选基函数的评价指标——相关因子，衡量了各项候选基函数对模型的贡献程度，通过一种基于边缘似然函数最大化的稀疏贝叶斯学习算法进行计算；系数矩阵回归模块，在给定的可调误差阈值下利用设计的评价指标对候选基函数进行最大化的筛除，得到该阈值条件下最稀疏的模型。该方法在给定的误差阈值下能辨识出最稀疏的模型；能通过提高误差阈值提炼出主要的基函数项，有利于对系统进行特征分析。

Description

一种稀疏自适应的复杂系统模型辨识方法
发明领域：
本发明为一种稀疏自适应的复杂系统模型辨识方法，涉及系统辨识、机器学习等领域。
背景技术：
随着信息技术的飞速发展，自动化技术发生了巨大的变化，现代控制理论的应用越来越广泛。对于一个不知道动态的陌生系统，人们需要对其建立有效的数学或物理模型以用于识别、分析和控制。系统辨识便因适应这一需要而诞生，其最重要的三个要素是数据、模型类和准则，本质就是按照一个准则在一组模型类中选择一个与数据拟合得最好的模型。系统辨识是以数据为基础,信息为手段,模型为媒体,以减少系统、信号、环境不确定性为目标的学科。它自发展以来在自动驾驶与机器人，复杂电子技术，工业过程监测，智能电网识别，电网故障检测，心房动作电位等领域都有广泛的应用。随着智能控制理论研究的不断深入，科学家们也在不断地进行各种尝试。系统辨识中的关键步骤是对建立的模型进行回归，经典的模型回归方法有很多，常见的有：利用神经网络进行深度学习、最小二乘法、Lasso回归——通过在最小二乘法的目标函数中加入l1范数的正则化项以达到对系数的稀疏回归等，然而这些方法存在以下不足：
1.得到的模型太复杂或者无模型，可扩展性弱，可解释性差。例如通过对LSTM、Transformer等神经网络进行深度学习的方法得出的模型十分复杂，带有众多的参数，缺乏可扩展性和可解释性，大多不方便用于对系统的调参控制。而小波网络等神经网络则无显式的模型，无法对模型进行数学或物理上的解释。
2.不能利用稀疏性简化模型复杂度。例如大多数神经网络的深度学习、最小二乘法等方法对于具有稀疏性的大多数物理系统，无法筛选出少数几个定义动力学的相关项，从而大大提高了模型复杂度。
3.难以实现在给定误差阈值条件下实现模型的进一步稀疏化。例如直接使用Lasso回归，由于Lasso正则化参数的取值具有连续性，较难通过调整正则化参数的取值在误差阈值下辨识出最稀疏的模型，往往易辨识出误差低于阈值条件很多，但复杂度仍较高的模型，而有时则需要更低的复杂度。
4.当候选基函数数量较多时，很多稀疏回归方法会导致模型的高复杂度。例如直接使用Lasso回归，当候选基函数增多时，一些与原始数据略有相关度但是次要的候选基函数项会被筛选出来，并且对于与其他基函数相似度较高的基函数项，Lasso回归不能对其进行有效的筛除，从而会导致模型的稀疏性降低，复杂度升高。
为此，本发明提供了一种稀疏自适应的复杂系统模型辨识方法，能有效缓解上述不足，在给定误差阈值下实现对系统的最优稀疏辨识，在候选基函数数目较多时仍能回归出较稀疏的模型，并且还能通过适当提高误差阈值来对基函数进行进一步提炼，有利于对系统进行特征分析。
发明内容
技术问题：
本发明旨在针对具有稀疏性的大规模复杂系统进行稀疏辨识，在人为设定的误差阈值下辨识出最稀疏的模型，在候选基函数数量较多时辨识出较稀疏的模型，能通过提高误差阈值来对基函数进行进一步提炼，实现对系统的有效识别、分析和控制。
技术方案：
本发明为一种稀疏自适应的复杂系统模型辨识方法，包含：基函数设计模块、评价指标求解模块、系数矩阵回归模块，共三个模块。
1.基函数设计模块
对于存在N个时间节点的一维时间序列数据x，其预测向量
Figure BDA0004133794740000026
考虑如下形式：
Figure BDA0004133794740000021
其中t是时间节点向量，共有N维，即t＝(t1,...,tN)T，ti(1≤i≤N)为第i个时间节点；预测向量
Figure BDA0004133794740000022
共有N维，即
Figure BDA0004133794740000023
为与第i个时间节点相对应的预测数据；目标向量(原始数据)x共有N维，即x＝(x1,...,xN)T，xi(1≤i≤N)为与第i个时间节点相对应的原始数据；Φ(t)为候选基函数库，它是一个N×M的矩阵，M是候选基函数的总个数，即
Figure BDA0004133794740000024
Figure BDA0004133794740000025
为设计的第j个候选基函数；w为系数矩阵，它是一个M×1的矩阵，即w＝(w1,...,wM)T，wj(1≤j≤M)为与第j个候选基函数
Figure BDA0004133794740000031
相对应的系数。
关于候选基函数
Figure BDA0004133794740000032
的选取，可以选用利于表现周期性的全频段三角基函数sin(2πkt/Δt),cos(2πkt/Δt)(1≤k≤K)(Δt为时间跨度,k为可能频率，K为最大可能频率，因为k一般为单位时间上的频率，所以需在三角基函数中除以Δt将实际时间转化为单位长度上的时间)、利于表现趋势性的幂函数tl(l≥0)等初等函数，以达到对系统的精准辨识。
2.评价指标求解模块
该模块设计了一种衡量各项候选基函数对模型贡献程度的重要评价指标——相关因子，并通过对系数矩阵进行基于边缘似然函数最大化的稀疏贝叶斯学习算法进行求解，该指标是模块3中进行候选基函数筛除的主要判据。
首先，假设数据的噪声ε服从独立的零均值、方差为σ2(β-1)的高斯分布:
Figure BDA0004133794740000033
其中N为时间节点的个数，β＝σ-2是噪声精度参数(噪声方差的倒数)。通过该条件，得到目标向量x的多元似然函数p(x|t,w,β)服从高斯分布
Figure BDA0004133794740000034
其中预测向量
Figure BDA0004133794740000035
由模块1中公式(1)给出。通过贝叶斯公式可知，后验概率正比于先验概率与似然函数的积。因此，需假设系数矩阵w的先验分布也服从高斯分布，才能得到系数矩阵w的后验分布也服从高斯分布，进而能将后验分布作为下一次的先验分布进行迭代学习。所以考虑如下形式的高斯先验：
Figure BDA0004133794740000036
其中M为候选基函数的总个数，先验精度参数向量α＝(α1,...,αM)T，即为每一个候选基函数
Figure BDA0004133794740000037
的系数wi都引入了一个单独的先验精度参数αi(先验分布方差的倒数)。
通过贝叶斯公式可得，系数w的后验分布p(w|x,t,α,β)服从高斯分布N(w|μ,Σ)，其中均值μ和方差矩阵Σ由如下公式给出：
Figure BDA0004133794740000038
其中A＝diag(α1,...,αM)，候选基函数库Φ(t)通过模块1设计并给出。
经过每一次迭代学习，通过后验分布重新得到了系数w的概率分布，进而将该后验分布作为新的先验分布继续迭代。对于与新先验分布和新似然函数相对应的先验精度参数α和噪声精度参数β，需通过最大化边缘似然函数p(x|α,β)＝∫p(x|w,β)p(w|α)dw来进行求解。由积分运算及取对数，得到边缘似然函数p(x|α,β)的对数值：
Figure BDA0004133794740000041
其中矩阵C＝β-1I+ΦA-1ΦT。
首先，令对数边缘似然函数L(α，β)关于噪声精度参数β的导数等于零，得到噪声精度参数β的最优值如下：
Figure BDA0004133794740000042
其中Σii表示方差矩阵Σ的第i个对角元素。
其次，为了方便对先验精度参数α的最优值进行快速求解，并引入评价指标——相关因子，需先为每一个候选基函数
Figure BDA0004133794740000043
引入稀疏度si和质量qi两个量：
Figure BDA0004133794740000044
通过如上公式可看出，稀疏度si衡量了基函数
Figure BDA0004133794740000045
与模型中其他基函数的重叠程度，质量qi衡量了基函数
Figure BDA0004133794740000046
与误差向量之间的对齐程度，这里的误差向量指目标向量x与基函数
Figure BDA0004133794740000047
从模型中被删除掉的预测向量
Figure BDA0004133794740000048
之间的差值。在对各项基函数的重要性进行评估时，往往考虑与其他基函数之间的相似度，以及对误差造成的影响大小，这正好与以上定义的稀疏度si和质量qi相对应。当si相对于qi较大时，候选基函数
Figure BDA0004133794740000049
与其他基函数的重叠程度较高，与误差的对齐程度较低，从而与其他基函数之间的相似度高，对误差造成的影响小，对模型的贡献程度也较小；当si相对于qi较小时，候选基函数
Figure BDA00041337947400000410
与其他基函数的重叠程度较低，与误差的对齐程度较高，从而与其他基函数之间的相似度低，对误差造成的影响大，对模型的贡献程度也较大。由此，通过稀疏度si和质量qi定义了一种衡量候选基函数
Figure BDA0004133794740000051
对模型贡献程度的评价指标——相关因子ri：
Figure BDA0004133794740000052
当相关因子ri较大时，与之对应的候选基函数
Figure BDA0004133794740000053
对模型贡献程度也较大；当相关因子ri较小时，与之对应的候选基函数
Figure BDA0004133794740000054
对模型贡献程度也较小。因此可用该指标来评价各项候选基函数的优劣性，便于后续进行候选基函数的筛除。
再来看对先验精度参数α的最优值的快速求解，可将对先验精度参数α的对数边缘似然函数L(α)表示为L(α)＝L(α-i)+λ(αi)，其中L(α-i)是省略了基函数
Figure BDA0004133794740000055
求得的对数边缘似然函数，中间参量λ(αi)由如下公式定义：
Figure BDA0004133794740000056
通过中间参量λ(αi)，可将对数边缘似然函数L(α)关于αi的导数转化成λ(αi)关于αi的导数，令导数
Figure BDA0004133794740000057
等于零，对数边缘似然函数L(α)取得最大值，先验精度参数αi取得最优解。
由公式(8)可看出，如果
Figure BDA0004133794740000058
此时αi→+∞，这是关于αi的解的第一种情况；如果
Figure BDA0004133794740000059
可解得αi的值如下：
Figure BDA00041337947400000510
这是关于αi的解的第二种情况。
由此得到了与新的先验分布和似然函数相对应的先验精度参数α和噪声精度参数β，同时还算得每个候选基函数
Figure BDA00041337947400000511
的相关因子ri。通过不断地迭代学习，每个候选基函数
Figure BDA00041337947400000512
的相关因子ri将会收敛到定值，可得到最终的相关因子，并以此作为模块3中候选基函数筛除的重要判据。
在实际操作中，往往根据上次迭代的先验精度参数α'，快速确定每个候选基函数
Figure BDA0004133794740000061
的稀疏度si和质量qi，进而求得所有候选基函数的相关因子r＝(r1,...,rM)，
若αi'＝∞，采用如下公式计算与候选基函数
Figure BDA0004133794740000062
相对应的稀疏度qi和质量si
Figure BDA0004133794740000063
若αi'≠∞，采用如下公式计算与候选基函数
Figure BDA0004133794740000064
相对应的稀疏度qi和质量si
Figure BDA0004133794740000065
紧接着利用公式(6)，可求得所有候选基函数的相关因子r。
该模块算法的具体步骤如下:
步骤1：使用无偏标准差S初始化噪声精度参数β的值，公式如下：
Figure BDA0004133794740000066
步骤2：选择与标签投影长度(余弦相似度)最大的候选基函数作为初始基函数
Figure BDA0004133794740000067
将其余候选基函数按原始下标的顺序排在其后。通过最大化边缘似然函数确定初始基函数
Figure BDA0004133794740000068
对应的先验精度参数α1，将其余j≠1的候选基函数的先验精度参数αj初始化为无穷大，公式如下：
Figure BDA0004133794740000071
步骤3：使用公式(2)计算所有候选基函数对应的均值μ和方差矩阵Σ。接着对于每一个候选基函数
Figure BDA0004133794740000072
使用公式(10)、公式(11)计算对应的稀疏度qi和质量si，使用公式(6)计算对应的相关因子ri，由此得到所有候选基函数的稀疏度、质量与相关因子；
步骤4：按下标的顺序(或随机)选择一个候选基函数
Figure BDA0004133794740000073
可重复多轮；
步骤5：对于步骤4中选择的候选基函数
Figure BDA0004133794740000074
若
Figure BDA0004133794740000075
且αi＜∞，或者
Figure BDA0004133794740000076
且αi＝∞，则使用公式(9)更新αi；若
Figure BDA0004133794740000077
且αi＜∞，则更新αi＝∞；
步骤6：使用公式(4)更新噪声精度参数β的值；
步骤7：若所有的先验精度参数αi都收敛或无αi可更新，则算法终止，所求即为此次迭代在步骤3中得到的相关因子r＝(r1,...,rM)；否则返回步骤3；
3.系数矩阵回归模块
该模块通过模块2中得到的相关因子r，在给定的误差阈值下最大化地筛除候选基函数，进而通过线性回归计算出最稀疏的系数矩阵，得到该误差阈值下最稀疏的模型，该模块算法具体步骤如下：
步骤1：给定一个误差阈值δ，含义是可容许的误差范围，即模型的误差需不超过该阈值。该阈值可以是平均绝对百分比误差MAPE、均方根误差RMSE等，这里默认使用平均绝对百分比误差MAPE。平均绝对百分比误差MAPE的定义如下：
Figure BDA0004133794740000078
其中xi为原始数据，
Figure BDA0004133794740000079
为预测值。该误差指标需越小越好，为0表示完美模型，大于1则表示劣质模型；
步骤2：将相关因子r＝(r1,...,rM)按照从小到大的顺序排列。首先，在不筛除候选基函数的前提下直接对系数矩阵进行线性回归，得到第1个模型，当前候选基函数库为Φ1＝Φ，对应的系数矩阵为w1，预测向量为
Figure BDA0004133794740000081
通过计算求得其平均绝对百分比误差为MAPE1；
其次，在排列后的相关因子r中找出最小的相关因子ri对应的候选基函数
Figure BDA0004133794740000082
并将其从当前的候选基函数库Φ1中筛除，得到更加稀疏的第2个模型，当前候选基函数库为Φ2，对应的系数矩阵为w2，预测向量为
Figure BDA0004133794740000083
平均绝对百分比误差为MAPE2；
接着，找出当前还未筛除的候选基函数所对应的相关因子中最小的相关因子rj所对应的候选基函数
Figure BDA0004133794740000084
并将其从当前的候选基函数库Φ2中筛除，得到更加稀疏的第3个模型，当前候选基函数库为Φ3，对应的系数矩阵为w3，预测向量为
Figure BDA0004133794740000085
平均绝对百分比误差为MAPE3；
以此类推，每次通过还未筛除的候选基函数所对应的相关因子中的最小相关因子选出对应的一个候选基函数，并将其从模型中筛除，因而每次回归得到的模型都比上一次回归得到的模型更加稀疏；
至于每次都要找最小的相关因子，是因为相关因子越小，其对应的候选基函数对模型的贡献程度就越小。从贡献程度小的候选基函数开始筛除，才能在确保精度的情况下逐步回归出愈加稀疏的模型；
最后，得到只剩下一个基函数
Figure BDA0004133794740000086
的第M个模型，对应的系数矩阵为wM，预测向量为
Figure BDA0004133794740000087
平均绝对百分比误差为MAPEM，再筛除该基函数后，所有候选基函数都被筛除，由此进入步骤3；
步骤3：在步骤2中所得模型的MAPE＝(MAPE1,...,MAPEM)中，找到小于等于误差阈值δ且最接近于误差阈值δ的平均绝对百分比误差MAPEi(1≤i≤M)，其对应的第i个模型，即为本发明所求解的最稀疏模型；若所得模型的MAPE全都大于误差阈值δ，则说明误差阈值δ设置的太小或者候选基函数库设计的太简单，需要适当提高误差阈值δ至大于第1个模型的MAPE1，或者再设计并增添候选基函数。该方法在给定的误差阈值下能辨识出最稀疏的模型；在候选基函数数量较多的情况下，仍能辨识出较稀疏的模型；能通过提高误差阈值提炼出主要的基函数项，有利于对系统进行特征分析。
有益效果：
1.本发明设计了一种针对各项候选基函数的评价指标，能够用于对候选基函数的筛除，从而能在给定的误差阈值下辨识出最稀疏的模型，极大地简化了模型的复杂度，使得模型能够更好地被用于系统的识别、分析和控制。
2.本发明与实际生产中系统的辨识、控制的需求相吻合，即在给定的误差范围内，对系统实行尽可能方便简单的识别、控制。
3.当候选基函数数目较多时，本发明中通过设定误差阈值来进行模型回归的方法能够保证该误差阈值下模型的最高稀疏性，对模型的复杂度具有良好的控制效果。
4.本发明可通过提高误差阈值继续降低模型的复杂度，从而进一步简化模型，提炼出主要的基函数项，便于系统的特征分析与控制的简化。
附图说明：
构成本申请的一部分的附图用来提供对本发明技术方案与有益效果的进一步理解，本发明的附图仅用于解释本发明，不用于限定本发明。在附图中：
图1为整体流程图；
图2为评价指标求解模块流程图；
图3为系数矩阵回归模块流程图；
图4为有益效果展示图。
具体实施方法：
为了使本发明的技术方案及有益效果更加清楚，下面结合具体实施例与说明书附图，对本发明进行进一步的详细说明。此处所描述的具体实施例仅用于解释本发明，不用于限定本发明。
以某地区一周的电力负荷量数据为实施例，选取每小时为采样点，使用本发明进行模型辨识的具体实施步骤如下：
步骤1：该步骤为基函数设计模块。选取的时间采样点构成了时间节点向量t，对应的采样数据构成了目标向量x，均为N维。由于该数据具有明显的周期性与趋势性特征，可设计以三角基函数、幂函数为底的候选基函数，具体流程如下：
步骤1.1：为了表现数据的周期性特征，先往候选基函数库中添加三角基函数，由采样点个数为N，将一周看作单位时间可得采样频率为N，考虑三角基函数的频率k范围为1～K(K＜N)，可设计如下形式的三角基函数：
Figure BDA0004133794740000101
其中Δt为时间跨度，因为频率k为单位时间上的频率，所以需在三角基函数中除以Δt将实际时间转化为单位长度上的时间，K为三角基函数的最大可能频率；
步骤1.2：为了表现数据的趋势性特征，再往候选基函数库中添加幂函数，形式如下：
Figure BDA0004133794740000102
其中l为幂函数的次数；
由此得到候选基函数库
Figure BDA0004133794740000103
步骤2：该步骤为评价指标求解模块，通过稀疏贝叶斯学习得到各个候选基函数的相关因子，反映了各个候选基函数对模型的贡献程度，可以此作为筛除候选基函数的重要评价指标，具体流程如下：
步骤2.1：使用无偏标准差初始化噪声精度参数β的值；
步骤2.2：选择与标签投影长度(余弦相似度)最大的候选基函数作为初始基函数
Figure BDA0004133794740000104
并确定对应的先验精度参数α1，其余候选基函数的先验精度参数αj被初始化为无穷大；
步骤2.3：计算所有候选基函数对应的均值μ和方差矩阵Σ。接着对于每一个候选基函数
Figure BDA0004133794740000105
计算对应的稀疏度qi、质量si与相关因子ri；
步骤2.4：按下标的顺序选择一个候选基函数
Figure BDA0004133794740000106
可重复多轮；
步骤2.5：对于步骤2.4中选择的候选基函数
Figure BDA0004133794740000111
若
Figure BDA0004133794740000112
且αi＜∞，或者
Figure BDA0004133794740000113
且αi＝∞，则使用更新αi的值；若
Figure BDA0004133794740000114
且αi＜∞，则更新αi＝∞；
步骤2.6：更新噪声精度参数β的值；
步骤2.7：若所有的先验精度参数αi都收敛或无αi可更新，则终止步骤2，所求即为此次迭代在步骤2.3中得到的相关因子r＝(r1,...,rM)；否则返回步骤2.3继续迭代；
步骤3：该步骤为系数矩阵回归模块，通过相关因子并利用线性回归，在给定误差阈值下对候选基函数进行最大化的筛除，得到给定误差阈值下的最稀疏模型，具体流程如下：
步骤3.1：给定一个误差阈值δ(MAPE值)，依次找出当前还未筛除的候选基函数所对应的相关因子中最小的相关因子ri所对应的候选基函数
Figure BDA0004133794740000115
并将其从候选基函数库Φ中筛除，通过线性回归得到越来越稀疏的模型，并计算每个模型的MAPE值(未筛除任何候选基函数时也要计算)，直到候选基函数全部被筛除为止；
步骤3.2：在所得的全部模型的MAPE值中，找到小于等于误差阈值δ且最接近于误差阈值δ的MAPE值，其对应的第i个模型，即为所求解的最稀疏模型；若所得模型的MAPE全都大于误差阈值δ，则适当提高误差阈值δ至大于第1个模型的MAPE值，或者再增添候选基函数。
此处结合实施例，对本发明提出的有益效果作详细展示：
1.给定误差阈值δ＝0.03，通过步骤1设计候选基函数至数目M＝50(三角基函数的可能频率范围取1～23，设计的三角基函数共46项，幂次项基函数取4项，次数分别为0,1,2,3)。通过步骤2与步骤3，结合相关因子指标对基函数作最大化的筛除，得到该阈值下的最稀疏模型。回归出的模型对应的MAPE值为0.029，基函数的总数目为8。可以看出，得到的模型在精度与复杂度上取得了良好的平衡。
2.扩大步骤1.1中三角基函数的最大可能频率K，增加候选基函数至数目M＝160并在误差阈值δ＝0.03下再次进行回归，得到的模型对应的MAPE值为0.029，基函数的总数目为10，通过基函数的总数目可看出模型复杂度与M＝50时得到的模型大体一致，这说明本发明在候选基函数较多时对模型复杂度具有良好的控制效果，不会导致基函数的总数目急剧增加。
3.候选基函数数目M＝50时，增大误差阈值至δ＝0.05并再次进行回归，得到的模型的MAPE值为0.041，基函数的总数目为4，较之前更加稀疏，同时提取出了主要的基函数项：x0、sin(2π·7t/Δt)、cos(2π·7t/Δt)、sin(2π·14t/Δt)共四项，与真实电力负荷量数据的周期性和趋势性基本吻合，说明本发明通过增大误差阈值提取出的主要基函数项可用于系统的特征分析。
可视化的有益效果展示见图4。
需要说明的是上述实施例，并非用来限定本发明的保护范围，在上述技术方案的基础上所作出的等同变换或替代均落入本发明权利要求所保护的范围。

Claims (7)

1.一种稀疏自适应的复杂系统模型辨识方法，其特征在于，所述方法中包含以下步骤：
步骤一：基函数设计模块，设计以初等函数为主的候选基函数库；
步骤二：评价指标求解模块，计算衡量各项候选基函数对模型贡献程度的评价指标——相关因子；
步骤三：系数矩阵回归模块，在给定的可调误差阈值下通过相关因子求解出最稀疏的模型。
2.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤一：基函数设计模块，设计以初等函数为主的候选基函数库，具体包括：在候选基函数库中设计并增添三角函数、幂函数等有利于表现数据周期性和趋势性的初等函数。
3.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤二：利用评价指标求解模块，计算衡量各项候选基函数对模型贡献程度的评价指标——相关因子，具体包括：定义了相关因子指标，并通过一种基于边缘似然函数最大化的稀疏贝叶斯学习算法进行求解。
4.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤三：系数矩阵回归模块，在给定的可调误差阈值下通过相关因子求解出最稀疏的模型，具体包括：根据步骤2中定义并求得的相关因子，在给定的误差阈值下最大化地筛除候选基函数，进而通过线性回归计算出最稀疏的系数矩阵。
5.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤一中基函数设计模块，具体如下：
对于存在N个时间节点的一维时间序列数据x，其预测向量
Figure FDA0004133794730000011
考虑如下形式：
Figure FDA0004133794730000012
其中t是时间节点向量，共有N维，即t＝(t1,...,tN)T，ti(1≤i≤N)为第i个时间节点；预测向量
Figure FDA0004133794730000013
共有N维，即
Figure FDA0004133794730000014
为与第i个时间节点相对应的预测数据；目标向量(原始数据)x共有N维，即x＝(x1,...,xN)T，xi(1≤i≤N)为与第i个时间节点相对应的原始数据；Φ(t)为候选基函数库，它是一个N×M的矩阵，M是候选基函数的总个数，即
Figure FDA0004133794730000021
Figure FDA0004133794730000022
为设计的第j个候选基函数；w为系数矩阵，它是一个M×1的矩阵，即w＝(w1,...,wM)T，wj(1≤j≤M)为与第j个候选基函数
Figure FDA0004133794730000023
相对应的系数，
关于候选基函数
Figure FDA0004133794730000024
的选取，选用利于表现周期性的全频段三角基函数sin(2πkt/Δt),cos(2πkt/Δt)(1≤k≤K),其中Δt为时间跨度,k为可能频率，K为最大可能频率，因为k一般为单位时间上的频率，所以需在三角基函数中除以Δt将实际时间转化为单位长度上的时间，利于表现趋势性的幂函数tl(l≥0)等初等函数，以达到对系统的精准辨识。
6.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤二：评价指标求解模块中，该模块算法的具体步骤如下:
步骤1：使用无偏标准差S初始化噪声精度参数β的值，公式如下：
Figure FDA0004133794730000025
步骤2：选择与标签投影长度(余弦相似度)最大的候选基函数作为初始基函数
Figure FDA0004133794730000026
将其余候选基函数按原始下标的顺序排在其后，通过最大化边缘似然函数确定初始基函数
Figure FDA0004133794730000027
对应的先验精度参数α1，将其余j≠1的候选基函数的先验精度参数αj初始化为无穷大，公式如下：
Figure FDA0004133794730000028
步骤3：使用公式(2)计算所有候选基函数对应的均值μ和方差矩阵Σ，接着对于每一个候选基函数
Figure FDA0004133794730000029
使用公式(10)、公式(11)计算对应的稀疏度qi和质量si，使用公式(6)计算对应的相关因子ri，由此得到所有候选基函数的稀疏度、质量与相关因子；
步骤4：按下标的顺序(或随机)选择一个候选基函数
Figure FDA0004133794730000031
可重复多轮；
步骤5：对于步骤4中选择的候选基函数
Figure FDA0004133794730000032
若
Figure FDA0004133794730000033
且αi＜∞，或者
Figure FDA0004133794730000034
且αi＝∞，则使用公式(9)更新αi；若
Figure FDA0004133794730000035
且αi＜∞，则更新αi＝∞；
步骤6：使用公式(4)更新噪声精度参数β的值；
步骤7：若所有的先验精度参数αi都收敛或无αi可更新，则算法终止，所求即为此次迭代在步骤3中得到的相关因子r＝(r1,...,rM)；否则返回步骤3。
7.根据权利要求1所述的稀疏自适应的复杂系统模型辨识方法，其特征在于，步骤三：系数矩阵回归模块，具体如下：
该模块通过模块2中得到的相关因子r，在给定的误差阈值下最大化地筛除候选基函数，进而通过线性回归计算出最稀疏的系数矩阵，得到该误差阈值下最稀疏的模型，该模块算法具体步骤如下：
步骤1：给定一个误差阈值δ，含义是可容许的误差范围，即模型的误差需不超过该阈值，该阈值可以是平均绝对百分比误差MAPE、均方根误差RMSE等，这里默认使用平均绝对百分比误差MAPE，平均绝对百分比误差MAPE的定义如下：
Figure FDA0004133794730000036
其中xi为原始数据，
Figure FDA0004133794730000037
为预测值，该误差指标需越小越好，为0表示完美模型，大于1则表示劣质模型；
步骤2：将相关因子r＝(r1,...,rM)按照从小到大的顺序排列，首先，在不筛除候选基函数的前提下直接对系数矩阵进行线性回归，得到第1个模型，当前候选基函数库为Φ1＝Φ，对应的系数矩阵为w1，预测向量为
Figure FDA0004133794730000038
通过计算求得其平均绝对百分比误差为MAPE1；
其次，在排列后的相关因子r中找出最小的相关因子ri对应的候选基函数
Figure FDA0004133794730000039
并将其从当前的候选基函数库Φ1中筛除，得到更加稀疏的第2个模型，当前候选基函数库为Φ2，对应的系数矩阵为w2，预测向量为
Figure FDA0004133794730000041
平均绝对百分比误差为MAPE2；
接着，找出当前还未筛除的候选基函数所对应的相关因子中最小的相关因子rj所对应的候选基函数
Figure FDA0004133794730000042
并将其从当前的候选基函数库Φ2中筛除，得到更加稀疏的第3个模型，当前候选基函数库为Φ3，对应的系数矩阵为w3，预测向量为
Figure FDA0004133794730000045
平均绝对百分比误差为MAPE3；
以此类推，每次通过还未筛除的候选基函数所对应的相关因子中的最小相关因子选出对应的一个候选基函数，并将其从模型中筛除，因而每次回归得到的模型都比上一次回归得到的模型更加稀疏；
至于每次都要找最小的相关因子，是因为相关因子越小，其对应的候选基函数对模型的贡献程度就越小，从贡献程度小的候选基函数开始筛除，才能在确保精度的情况下逐步回归出愈加稀疏的模型；
最后，得到只剩下一个基函数
Figure FDA0004133794730000043
的第M个模型，对应的系数矩阵为wM，预测向量为
Figure FDA0004133794730000044
平均绝对百分比误差为MAPEM，再筛除该基函数后，所有候选基函数都被筛除，由此进入步骤3；
步骤3：在步骤2中所得模型的MAPE＝(MAPE1,...,MAPEM)中，找到小于等于误差阈值δ且最接近于误差阈值δ的平均绝对百分比误差MAPEi(1≤i≤M)，其对应的第i个模型，即为本发明所求解的最稀疏模型；若所得模型的MAPE全都大于误差阈值δ，则说明误差阈值δ设置的太小或者候选基函数库设计的太简单，需要适当提高误差阈值δ至大于第1个模型的MAPE1，或者再设计并增添候选基函数。
CN202310268392.9A 2023-03-20 2023-03-20 一种稀疏自适应的复杂系统模型辨识方法 Pending CN116166921A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310268392.9A CN116166921A (zh) 2023-03-20 2023-03-20 一种稀疏自适应的复杂系统模型辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310268392.9A CN116166921A (zh) 2023-03-20 2023-03-20 一种稀疏自适应的复杂系统模型辨识方法

Publications (1)

Publication Number Publication Date
CN116166921A true CN116166921A (zh) 2023-05-26

Family

ID=86416452

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310268392.9A Pending CN116166921A (zh) 2023-03-20 2023-03-20 一种稀疏自适应的复杂系统模型辨识方法

Country Status (1)

Country Link
CN (1) CN116166921A (zh)

Similar Documents

Publication Publication Date Title
CN110414788B (zh) 一种基于相似日和改进lstm的电能质量预测方法
CN109871977A (zh) 基于离散小波变换和优化最小支持向量机的负荷预测方法
CN110571792A (zh) 一种电网调控系统运行状态的分析评估方法及系统
CN111768000A (zh) 在线自适应微调深度学习的工业过程数据建模方法
CN109146186A (zh) 一种基于双重分解的短期风电功率预测方法
CN114117852B (zh) 一种基于有限差分工作域划分的区域热负荷滚动预测方法
CN111652271A (zh) 一种基于神经网络的非线性特征选择方法
CN110633846A (zh) 燃气负荷预测方法和装置
CN111985825A (zh) 一种用于滚磨机定向仪的晶面质量评估方法
CN113780420B (zh) 基于gru-gcn的变压器油中溶解气体浓度预测方法
CN113743504A (zh) 一种电压暂降源识别模型构建方法、装置、终端及介质
CN112508286A (zh) 一种基于Kmeans-BiLSTM-DMD模型的短期负荷预测方法
CN118151020B (zh) 一种电池安全性能的检测方法及系统
Wang et al. Time-weighted kernel-sparse-representation-based real-time nonlinear multimode process monitoring
CN113052373A (zh) 一种基于改进elm模型的月径流变化趋势预测方法
CN113420815B (zh) 半监督rsdae的非线性pls间歇过程监测方法
CN114500335A (zh) 基于模糊c均值和混合核最小二乘支持向量机的sdn网络流程控制方法
CN117592593A (zh) 基于改进二次模态分解和WOA优化BILSTM-attention的短期电力负荷预测方法
CN117458480A (zh) 基于改进lof的光伏发电功率短期预测方法及系统
Rocha Filho et al. Recursive fuzzy instrumental variable based evolving neuro-fuzzy identification for non-stationary dynamic system in a noisy environment
CN110738363A (zh) 一种光伏发电功率预测模型及其构建方法和应用
JPH06337852A (ja) ニューラルネットワークによる時系列予測方法
CN111369072A (zh) 一种基于稀疏化方法的核最小均方时间序列在线预测模型
CN116596396A (zh) 一种基于k近邻插值和slstm的工业聚乙烯过程质量预测方法
CN116166921A (zh) 一种稀疏自适应的复杂系统模型辨识方法

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