CN108053077A - 一种基于区间二型t-s模糊模型的短期风速预测方法与系统 - Google Patents

一种基于区间二型t-s模糊模型的短期风速预测方法与系统 Download PDF

Info

Publication number
CN108053077A
CN108053077A CN201711463836.5A CN201711463836A CN108053077A CN 108053077 A CN108053077 A CN 108053077A CN 201711463836 A CN201711463836 A CN 201711463836A CN 108053077 A CN108053077 A CN 108053077A
Authority
CN
China
Prior art keywords
model
fuzzy
output
data
input
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
CN201711463836.5A
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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and 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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201711463836.5A priority Critical patent/CN108053077A/zh
Publication of CN108053077A publication Critical patent/CN108053077A/zh
Pending legal-status Critical Current

Links

Classifications

    • 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
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/02—Computing arrangements based on specific mathematical models using fuzzy logic
    • 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/06—Energy or water supply

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • Economics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Health & Medical Sciences (AREA)
  • Strategic Management (AREA)
  • Software Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • General Business, Economics & Management (AREA)
  • Tourism & Hospitality (AREA)
  • Marketing (AREA)
  • Fuzzy Systems (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Biomedical Technology (AREA)
  • Automation & Control Theory (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Public Health (AREA)
  • Water Supply & Treatment (AREA)
  • Molecular Biology (AREA)
  • Data Mining & Analysis (AREA)
  • Primary Health Care (AREA)
  • Computational Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Algebra (AREA)
  • Development Economics (AREA)
  • Game Theory and Decision Science (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开一种基于区间二型T‑S模糊模型的短期风速预测方法与系统。方法中对历史实际风速观测数据进行变分模态分解VMD为K个模态，对每个模态进行属性选择和归一化处理，建立预测模糊模型；利用区间二型模糊C回归聚类IT2‑FCR对模型进行结构划分，并以实际观测风速数据与模糊模型预测结果的加权均方根误差为目标函数利用引力搜索算法GSA进行模型前件参数优化；利用最小二乘法对模型参数进行辨识。从而得到以历史实际观测数据为输入的短期风速预测区间二型T‑S模糊模型。本发明采用一种新型的二型超平面隶属函数，可以提高风速模糊预测模型的辨识精度，能得到更精确的辨识参数，使得对应的风速预测结果与实际观测风速吻合更好。

Description

一种基于区间二型T-S模糊模型的短期风速预测方法与系统
技术领域
本发明属于风速短期预测技术领域，更具体地，涉及一种基于区间二型T-S模糊模型的短期风速预测方法与系统。
背景技术
传统的化石能源日益枯竭、全球变暖等问题使得新型清洁能源的利用日益迫切。风电作为一种蕴量巨大、转换快速且无公害的能源系统，中国新能源战略开始把大力发展风力发电设为重点。按照国家规划，未来15年，全国风力发电装机容量将达到2000万至3000万千瓦。风电能源的不稳定性和不可控性给电能稳定利用带来严峻挑战。有效的风速预测机制能为风电场检修计划以及电力系统调度计划提供可靠依据，保证电力系统安全稳定运行。另一方面准确的风速预测能提高风电场的能源利用效率、减轻并网运行时对电网的不利影响、降低风电场运行成本。因此，对风速的精确预测显得尤为重要。
发明内容
本发明的目的是针对风能的不稳定性导致对电能利用效率的不良影响，设计了一种基于区间二型T-S模糊模型的短期风速预测方法与系统，结合GSA优化算法能有效提高T-S模糊预测模型的辨识精度，减小风速预测的预测误差。
为了实现上述目的，本发明提供了一种基于区间二型T-S模糊模型的短期风速预测方法，如图2所示，包括：
步骤1、风速观测时间序列数据data预处理，建立T-S模糊预测模型的输入输出矩阵(input,output)：
步骤1.1、利用变分模态分解VMD将风速观测时间序列数据data分解为K个本征模态函数IMF，对每个IMF都独立执行下列步骤；
步骤1.2、取IMF前L间隔数据分量划分候选输入属性矩阵In，其中Output为风速实际观测数据矩阵，LL为风速观测时间序列长度，利用克莱姆-施密特正交化方法GSO进行候选输入属性对实际观测值的相关性排序，从中挑选前M(M≤L)个属性，建立风速预测T-S模糊模型的数据对(Input,Output)；
步骤1.3、将数据对(Input,Output)进行数据归一化处理，归一化到区间[-1,1]内，得到归一化后的数据对(input,output)；
步骤2、建立区间二型T-S模糊模型，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个模糊系统即为各个子系统的线性组合。所述IF-THEN模糊规则如下：
规则i:
i＝1…C,k＝1…N,j＝1…M
其中是线性隶属度函数的模糊集合，缩写为前件参数、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,…xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，指代超平面。再通过线性组合得到T-S模糊模型。其中记 其中ωi为高斯模糊权重；
步骤3、将数据对(input,output)分为三部分，第一部分用作模型初始化(xini,yini)，第二部分用作模型优化(xopt,yopt)，第三部分用作模型测试(xtest,ytest)。利用第一部分数据，通过区间二型模糊C回归聚类(Inter Type-2Fuzzy C-regression，IT2-FCR)进行T-S模糊模型前件参数辨识，得到上下超平面参数和α：
步骤3.1、算法初始化：设置算法参数，包括两个模糊因子m1和m2、聚类总数C、距离调节参数η、最大迭代次数Tmax、停止迭代阈值ε。特别地，利用模糊C均值进行粗聚类获取一型隶属度矩阵U＝[uik]，令初始 和u ik是xk分别隶属于第i类上超平面和下超平面的隶属度，初始划分矩阵令当前迭代次数t＝0；
步骤3.2、利用下式计算每个聚类的上下超平面参数：
α i＝[XT P iX]-1XT P iy
其中，X＝[x 1]∈RN×(M+1)
步骤3.3、利用类型还原并计算实际观测值与模糊模型预测输出的误差Eik(αi)＝yk-fi(xk,αi)，其中yk是实际观测矩阵y＝(y1,y2…yN)T的第k分量；
步骤3.4、按照下式计算上下隶属度函数和u ik：
步骤3.5、利用上下隶属度函数更新计算划分矩阵和P；
步骤3.6、令t＝t+1，转至步骤3.2，直到相邻两次迭代的差值 是当前迭代的第i类超平面参数，是上一次迭代的第i类超平面参数，ε是停止迭代阈值，或者t＞Tmax时，停止迭代；
步骤3.7、按照上述步骤计算每个聚类的上下超平面参数；
步骤4、利用引力搜索算法GSA进行模型上下超平面参数优化，得到优化后上下超平面参数：
步骤4.1、算法初始化，设置GSA算法参数，包括最大迭代次数Nitmax、群体规模N0、初始引力常数G0，衰减系数β等。选取模型上下超平面参数和α为待优化变量确定待优化变量的上下界[BL,BU]， 和α max分别为上下超平面参数的最小值和最大值。在此区间随机初始化群体中所有个体的位置向量，个体位置向量代表一组控制参数；令当前迭代次数t＝0；
步骤4.2、计算群体目标函数最小值，具有最小目标函数值的个体确定为当前最优个体Xbest；
步骤4.2.1、个体i位置向量Xi(t)解码得到控制参数，分别为模型上超平面参数和下超平面参数，基于超平面型隶属函数计算高斯模糊权重ωi：
步骤4.2.1.1、计算输入数据点距离上下超平面的距离：
步骤4.2.1.2、依据下式计算数据点分别隶属于上下超平面型的高斯隶属度函数和
步骤4.2.1.3、计算模糊权重ωi(xk)表示第k个输入分量属于第i类的程度。
步骤4.2.2、最小二乘法辨识模糊后件参数：利用第一部分数据得到Α＝(ΦTΦ)- 1XTy；
步骤4.2.3、分别以第一和第二部分数据计算模型预测输出和Φopt是按照第二部分数据计算得到。反归一化后按照下式分别计算两部分数据的模型预测值与实际观测值的均方根误差RMSEtrain和RMSEopt：
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值。
以两部分数据所计算得到的RMSE的加权平均值作为个体i的目标函数值：
fiti＝t1*RMSEtrain+t2*RMSEopt
其中，t1+t2＝1。
步骤4.2.4、重复步骤4.2.1-4.2.3，获取种群中N0个个体的目标函数值；
步骤4.3、更新所有个体的引力常数Gi；
其中，G0为引力常数初始值，β为衰减系数，t为当前迭代次数，Nitmax为最大迭代次数。
步骤4.4、计算所有个体的万有引力Fi和加速度ai；
依据牛顿引力定理，第i个粒子受到第j个粒子的作用力为：
其中，Maj是第j个粒子主动引力质量，Mpi是第i个粒子的被动引力质量，G(t)是引力时间常数，此时认为其为时变量。
对第i个粒子，受到来自其他粒子引力合力用引力的随机加权和表示为：
基于牛顿第二定理，粒子i产生的加速度为：
其中，Mii为粒子i的惯性质量；
引力及惯性质量则依据目标函数值计算。质量重的个体较质量轻的个体更为优秀。假设引力质量与惯性质量相等，根据目标函数给出粒子质量，其定义为：
Mai＝Mpi＝Mii＝Mi
对于极小化问题best＝min fitj，worst＝max fitj；
步骤4.5、更新所有个体的速度vi和位置Xi；
randi为(0，1)之间随机数，d表示位置向量的第d维；
步骤4.6、t＝t+1，如果t>Nitmax，当前最优个体位置作为终解；否则，转入步骤4.2。
步骤4.7、输出最优个体位置即为模型上下超平面的最优参数、优化后的最优模型后件参数Α和模糊权重ωi。
步骤5、得到完整的模型前件参数及后件参数，利用第三部分测试数据的当前IMF计算模型预测输出，Φtest是按照第三部分测试数据计算得到。
步骤6、每一个模态进行上述步骤1.2-步骤5，反归一化后合并得到第三部分数据的模型预测输出；
步骤7、最后以模糊模型预测输出与实际观测数据均方根误差RMSE，平均绝对值误差MAE以及平均绝对百分误差MAPE为指标，来评估模糊模型预测的准确度。其中，均方根误差RMSE，平均绝对误差MAE和平均绝对百分误差MAPE的表达式如下：
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值。
与现有技术相比，利用本发明所述方法进行短期风速预测时，IT2-FCRM可以辨识得到T-S模糊模型精确参数，更精确的预测模型辨识参数使得对应的模型预测结果精度更高，且过程简单，计算量小。
附图说明
图1为本发明实施例提供的一种基于区间二型T-S模糊模型的短期风速预测方法流程图；
图2为本发明实施例提供的一种具体的基于区间二型T-S模糊模型的短期风速预测方法流程图；
图3为本发明实施例提供的另一种完整的基于区间二型T-S模糊模型的短期风速预测方法流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
实施例1：
本发明实施例提供了一种基于区间二型T-S模糊模型的短期风速预测方法，如图1所示，包括以下执行步骤：
在步骤1中，风速观测数据data预处理，通过变分模态分解(Variational ModeDecomposition，简写为：VMD)将原始风速观测数据分解为K个本征模态函数(IntrinsicMode Function，简写为：IMF)，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵(input,output)；
在步骤2中，对应于每一个本征模态函数IMF，建立区间二型T-S模糊模型，其中，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个非线性系统为各个子系统的线性组合；所述IF-THEN模糊规则如下：
规则i：
i＝1…C,k＝1…N,j＝1…M
其中是线性隶属度函数的模糊集合、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,…xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，用于指代超平面；
再通过线性组合得到T-S模糊模型；
其中， ωi为高斯模糊权重；
在步骤3中，根据(input,output)数据对，通过区间二型模糊C回归聚类(InterType-2Fuzzy C-regression，简写为IT2-FCR)，进行T-S模糊模型参数辨识，得到上下超平面参数和α；
在步骤4中，根据(input,output)数据对，利用引力搜索算法(GravitationalSearch Algorithm，简写为：GSA)进行模型上下超平面参数优化，得到优化后上下超平面参数；
在步骤5中，得到完整的模型前件参数及后件参数，利用(input,output)数据对计算当前IMF计算模型预测输出；
在步骤6中，对于每一个本征模态函数IMF建立各自的T-S模糊预测模型，并计算相应IMF计算模型预测输出；进一步，反归一化后合并得到(input,output)数据对的模型预测输出。
本发明实施例所述方法进行短期风速预测时，IT2-FCRM可以辨识得到T-S模糊模型精确参数，更精确的预测模型辨识参数使得对应的模型预测结果精度更高，且过程简单，计算量小。
在本发明实施例中，通过变分模态分解VMD将原始风速观测数据分解为K个本征模态函数IMF，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵，具体包括：对每个IMF都执行下列内容：
对IMF前L间隔数据分量划分候选输入属性矩阵In，其中Output为风速实际观测数据矩阵，LL为风速观测时间序列长度，利用克莱姆-施密特正交化方法(Gram-Schmidt Orthogonal，简写为：GSO)进行候选输入属性对实际输出值的相关性排序，从中挑选前M个属性，建立风速预测T-S模糊模型的数据对(Input,Output)；其中，M≤L。
在本发明实施例中，针对于所述风速短期预测方法所适用的场合的不同，数据对(Input,Output)也将被做不同的拆分。具体的，在所述方法进行模型精确度的测试时，所述数据对(Input,Output)被分为三部分，包括：第一部分数据用作模型初始化(xini,yini)，第二部分数据用作模型优化(xopt,yopt)，第三部分数据用作模型测试(xtest,ytest)；
在所述方法进行实际模型使用时，所述数据对(Input,Output)被分为两部分，包括第一部分数据用作模型初始化(xini,yini)，第二部分数据用作模型优化(xopt,yopt)。
结合本发明实施例，存在一种可扩展的实现方案，其中，所述利用区间二型模糊C回归聚类IT2-FCR进行T-S模糊模型参数辨识，得到上下超平面参数和α，具体包括：
步骤3.1、设置算法参数，所述算法参数包括两个模糊因子m1和m2、聚类总数C、距离调节参数η、最大迭代次数Tmax和停止迭代阈值ε；
利用模糊C均值进行粗聚类获取一型隶属度矩阵U＝[uik]，令初始 和u ik是xk分别隶属于第i类上超平面和下超平面的隶属度，初始划分矩阵令当前迭代次数t＝0；
步骤3.2、利用下式计算每个聚类的上下超平面参数：
α i＝[XT P iX]-1XT P iy
其中，X＝[x 1]∈RN×(M+1)；
步骤3.3、利用类型还原并计算实际观测值与模糊模型预测输出的误差Eik(αi)＝yk-fi(xk,αi)，其中yk是实际观测矩阵y＝(y1,y2…yN)T的第k分量；
步骤3.4、按照下式计算上下隶属度函数：
步骤3.5、利用上下隶属度更新计算划分矩阵和P；
步骤3.6、令t＝t+1，转至步骤3.2，直到相邻两次迭代的差值||αi current-αi last||＞ε，或者t＞Tmax时，停止迭代；其中，αi current是当前迭代的第i类超平面参数，αi last是上一次迭代的第i类超平面参数；
步骤3.7、计算每个聚类的上下超平面参数。
结合本发明实施例，存在一种可扩展的实现方案，其中，所述采用引力搜索算法GSA进行模型超平面参数优化，具体包括：
步骤4.1、算法初始化，设置GSA算法参数，包括最大迭代次数Nitmax、群体规模N0、初始引力常数G0，衰减系数β；
步骤4.2、计算群体目标函数最小值，具有最小目标函数值的个体确定为当前最优个体Xbest；
步骤4.3、更新所有个体的引力常数Gi；
其中，G0为引力常数初始值，β为衰减系数，t为当前迭代次数，Nitmax为最大迭代次数；
步骤4.4、计算所有个体的万有引力Fi和加速度ai；
步骤4.5、更新所有个体的速度vi和位置Xi；
步骤4.6、t＝t+1，如果t>Nitmax，当前最优个体位置作为终解；否则，转入步骤4.2；
步骤4.7、输出最优个体位置即为模型上下超平面的最优参数、优化后的最优模型后件参数Α和模糊权重ωi。
结合本发明实施例，存在一种可扩展的实现方案，其中，所述计算群体目标函数最小值，具体包括：
步骤4.2.1、个体i位置向量Xi(t)解码得到控制参数，分别为模型上超平面参数和下超平面参数，基于超平面型隶属函数计算高斯模糊权重ωi；
步骤4.2.2、最小二乘法辨识模糊后件参数：利用数据对(Input,Output)的第一部分数据得到Α＝(ΦTΦ)-1XTy；
步骤4.2.3、分别以数据对(Input,Output)的第一部分数据和数据对(Input,Output)的第二部分数据计算模型预测输出和Φopt是按照第二部分数据计算得到；反归一化后分别计算模型预测值与实际观测值的均方根误差RMSEtrain和RMSEopt:
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值；
以两部分数据所计算得到的RMSE的加权平均作为个体i的目标函数值：
fiti＝t1*RMSEtrain+t2*RMSEopt
其中，t1+t2＝1。
结合本发明实施例，存在一种可扩展的实现方案，其中，所述基于超平面型隶属函数计算高斯模糊权重ωi，具体包括：
步骤4.2.1.1、计算输入数据点距离上下超平面的距离：
步骤4.2.1.2、依据下式计算数据点分别隶属于上下超平面型的高斯隶属度函数和
步骤4.2.1.3、计算模糊权重ωi(xk)表示第k个输入分量属于第i类的程度。
结合本发明实施例，存在一种可扩展的实现方案，其中，对所述方法进行模型精确度的测试时，所述方法还包括：
以模糊模型预测输出与实际观测数据均方根误差RMSE，平均绝对误差MAE以及平均绝对百分误差MAPE为指标，来评估模糊模型预测的准确度，其中，均方根误差RMSE，平均绝对误差MAE和平均绝对百分误差MAPE的表达式如下：
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值；
根据所述RMSE、MAE和MAPE的结果数据，评判当前模糊模型的准确度。其中，RMSE用来衡量模型预测值同实际观测值之间的偏差，对模型预测中特大或特小误差反映非常敏感，因此能够很好地反映预测精密度，RMSE的值越小说明预测模型具有更好的精确度。MAE有利于不同模型之间的比较。MAPE用来衡量模型预测结果好坏，值越小精度越高。
实施例2：
本发明实施例还提供了一种基于区间二型T-S模糊模型的短期风速预测系统，系统包括风速检测器和服务器，所述风速检测器用于将风速观测数据data发送给服务器，具体的：
所述服务器，用于风速观测数据data预处理，通过变分模态分解VMD将原始风速观测数据分解为K个本征模态函数IMF，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵(input,output)；
对应于每一个本征模态函数IMF，建立区间二型T-S模糊模型，其中，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个非线性系统为各个子系统的线性组合；所述IF-THEN模糊规则如下：
规则i：
i＝1…C,k＝1…N,j＝1…M
其中是线性隶属度函数的模糊集合、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,…xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，用于指代超平面；
再通过线性组合得到T-S模糊模型；其中，记 ωi为高斯模糊权重；
根据(input,output)数据对，通过区间二型模糊C回归聚类IT2-FCR，进行T-S模糊模型参数辨识，得到上下超平面参数和α；
根据(input,output)数据对，利用引力搜索算法GSA进行模型上下超平面参数优化，得到优化后上下超平面参数；
得到完整的模型前件参数及后件参数，利用(input,output)数据对计算当前IMF计算模型预测输出；
对于每一个本征模态函数IMF建立各自的T-S模糊预测模型，并计算相应IMF计算模型预测输出；进一步，反归一化后合并得到(input,output)数据对的模型预测输出。
值得说明的是，上述系统内的模块、单元之间的信息交互、执行过程等内容，由于与本发明的处理方法实施例基于同一构思，具体内容可参见本发明方法实施例1中的叙述，此处不再赘述。
实施例3：
为说明本发明效果，下面以西班牙西北部加利西亚自治区的Sotavento风电场在2014年1月5-11号每10分钟记录一次的风速观测数据作为本发明的实施对象对本发明方法进行详细说明，数据共计1008条。并且，本发明实施例以所述方法进行模型精确度的测试适应场景，具体过程如下：
步骤1、风速观测时间序列数据data预处理，建立T-S模糊预测模型的输入输出矩阵(input,output)：
步骤1.1、利用变分模态分解VMD将风速观测时间序列数据data分解为K＝5个本征模态函数IMF，对每个IMF都独立执行下列步骤；
步骤1.2、取IMF前L＝24间隔数据分量划分候选输入属性矩阵In，其中Output为风速实际观测数据矩阵，LL＝1008为风速观测时间序列长度，利用克莱姆-施密特正交化方法GSO进行候选输入属性对实际观测值的相关性排序，从中挑选前M＝9个属性，建立风速预测T-S模糊模型的数据对(Input,Output)；
步骤1.3、将数据对(Input,Output)进行数据归一化处理，归一化到区间[-1,1]内，得到归一化后的数据对(input,output)；
步骤2、建立区间二型T-S模糊模型，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个模糊系统即为各个子系统的线性组合。所述IF-THEN模糊规则如下：
规则i:
i＝1…C,k＝1…N,j＝1…M
其中是线性隶属度函数的模糊集合，缩写为前件参数、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,…xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，指代超平面。再通过线性组合得到T-S模糊模型。其中
记 其中ωi为高斯模糊权重；
步骤3、将数据对(input,output)分为三部分，第一部分用作模型初始化(xini∈R552×9,yini)，第二部分用作模型优化(xopt∈R144×9,yopt)，第三部分用作模型测试(xtest∈R288 ×9,ytest)。利用第一部分数据，通过区间二型模糊C回归聚类IT2-FCR进行T-S模糊模型前件参数辨识，得到上下超平面参数和α：
步骤3.1、算法初始化：设置算法参数，包括两个模糊因子m1＝1.5和m2＝7、聚类总数C＝4、距离调节参数η＝5、最大迭代次数Tmax＝100、停止迭代阈值ε＝1e-6。特别地，利用模糊C均值进行粗聚类获取一型隶属度矩阵U＝[uik]，令初始 和u ik是xk分别隶属于第i类上超平面和下超平面的隶属度，初始划分矩阵 令当前迭代次数t＝0；
步骤3.2、利用下式计算每个聚类的上下超平面参数：
α i＝[XT P iX]-1XT P iy
其中，X＝[x 1]∈RN×(M+1)；
步骤3.3、利用类型还原并计算实际观测值与模糊模型预测输出的误差Eik(αi)＝yk-fi(xk,αi)，其中yk是实际观测矩阵y＝(y1,y2…yN)T的第k分量；
步骤3.4、按照下式计算上下隶属度函数和u ik：
步骤3.5、利用上下隶属度函数更新计算划分矩阵和P；
步骤3.6、令t＝t+1，转至步骤3.2，直到相邻两次迭代的差值 是当前迭代的第i类超平面参数，αi last是上一次迭代的第i类超平面参数，ε是停止迭代阈值，或者t＞Tmax时，停止迭代；
步骤3.7、按照上述步骤计算每个聚类的上下超平面参数；
步骤4、利用引力搜索算法GSA进行模型上下超平面参数优化，得到优化后上下超平面参数：
步骤4.1、算法初始化，设置GSA算法参数，包括最大迭代次数Nitmax＝200、群体规模N0＝20、初始引力常数G0＝100，衰减系数β＝20等。选取模型上下超平面参数和α为待优化变量确定待优化变量的上下界[BL,BU]，和α max分别为上下超平面参数的最小值和最大值。在此区间随机初始化群体中所有个体的位置向量，个体位置向量代表一组控制参数；令当前迭代次数t＝0；
步骤4.2、计算群体目标函数最小值，具有最小目标函数值的个体确定为当前最优个体Xbest；
步骤4.2.1、个体i位置向量Xi(t)解码得到控制参数，分别为模型上超平面参数和下超平面参数，基于超平面型隶属函数计算高斯模糊权重ωi；
步骤4.2.1.1、计算输入数据点距离上下超平面的距离：
步骤4.2.1.2、依据下式计算数据点分别隶属于上下超平面型的高斯隶属度函数和
步骤4.2.1.3、计算模糊权重ωi(xk)表示第k个输入分量属于第i类的程度。
步骤4.2.2、最小二乘法辨识模糊后件参数：利用第一部分数据得到Α＝(ΦTΦ)- 1XTy；
步骤4.2.3、分别以第一和第二部分数据计算模型预测输出和Φopt是按照第二部分数据计算得到。反归一化后按照下式分别计算两部分数据的模型预测值与实际观测值的均方根误差RMSEtrain和RMSEopt：
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值。
以两部分数据所计算得到的RMSE的加权平均作为个体i的目标函数值：
fiti＝t1*RMSEtrain+t2*RMSEopt
其中，t1+t2＝1。
步骤4.2.4、重复步骤4.2.1-4.2.3，获取种群中N0个个体的目标函数值；
步骤4.3、更新所有个体的引力常数Gi；
其中，G0为引力常数初始值，β为衰减系数，t为当前迭代次数，Nitmax为最大迭代次数。
步骤4.4、计算所有个体的万有引力Fi和加速度ai；
依据牛顿引力定理，第i个粒子受到第j个粒子的作用力为：
其中，Maj是第j个粒子主动引力质量，Mpi是第i个粒子的被动引力质量，G(t)是引力时间常数，此时认为其为时变量。
对第i个粒子，受到来自其他粒子引力合力用引力的随机加权和表示为：
基于牛顿第二定理，粒子i产生的加速度为：
其中，Mii为粒子i的惯性质量；
引力及惯性质量则依据目标函数值计算。质量重的个体较质量轻的个体更为优秀。假设引力质量与惯性质量相等，根据目标函数给出粒子质量，其定义为：
Mai＝Mpi＝Mii＝Mi
对于极小化问题best＝min fitj，worst＝max fitj；
步骤4.5、更新所有个体的速度vi和位置Xi；
randi为(0，1)之间随机数，d表示位置向量的第d维；
步骤4.6、t＝t+1，如果t>Nitmax，当前最优个体位置作为终解；否则，转入步骤4.2。
步骤4.7、输出最优个体位置即为模型上下超平面的最优参数、优化后的最优模型后件参数Α和模糊权重ωi。
步骤5、得到完整的模型前件参数及后件参数，利用第三部分测试数据的当前IMF计算模型预测输出，Φtest是按照第三部分测试数据计算得到。
步骤6、每一个模态进行上述步骤1.2-步骤5，反归一化后合并得到第三部分数据的模型预测输出；
步骤7、最后以模糊模型预测输出与实际观测数据均方根误差RMSE，平均绝对值误差MAE以及平均绝对百分误差MAPE为指标，来评估模糊模型预测的准确度。其中，均方根误差RMSE，平均绝对误差MAE和平均绝对百分误差MAPE的表达式如下：
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值。
表1为本发明所述方法的超平面原型的上下超平面参数。表2为与ARIMA预测方法对比，基于第三部分测试数据所得的预测精度指标。从表3可以看出，对于短期风速预测，本发明方法的预测精度较高。
表1超平面聚类参数
表2基于测试数据与ARIMA预测方法的精度指标对比
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，方法包括：
风速观测数据data预处理，通过变分模态分解VMD将原始风速观测数据分解为K个本征模态函数IMF，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵(input,output)；
对应于每一个本征模态函数IMF，建立区间二型T-S模糊模型，其中，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个非线性系统为各个子系统的线性组合；所述IF-THEN模糊规则如下：
规则i：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msup> <mi>f</mi> <mi>i</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>,</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>a</mi> <mn>1</mn> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>a</mi> <mn>2</mn> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mn>2</mn> </mrow> </msub> <mo>+</mo> <mn>...</mn> <mo>+</mo> <msubsup> <mi>a</mi> <mi>M</mi> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mi>M</mi> </mrow> </msub> <mo>+</mo> <msubsup> <mi>a</mi> <mn>0</mn> <mi>i</mi> </msubsup> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>k</mi> </msub> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;alpha;</mi> <mrow> <mi>i</mi> <mi>T</mi> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>i</mi> <mo>=</mo> <mn>1...</mn> <mi>C</mi> <mo>,</mo> <mi>k</mi> <mo>=</mo> <mn>1...</mn> <mi>N</mi> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1...</mn> <mi>M</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中是线性隶属度函数的模糊集合、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,···xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，用于指代超平面；
再通过线性组合得到T-S模糊模型；其中，
记 ωi为高斯模糊权重；
根据(input,output)数据对，通过区间二型模糊C回归聚类IT2-FCR，进行T-S模糊模型参数辨识，得到上下超平面参数和α；
根据(input,output)数据对，利用引力搜索算法GSA进行模型上下超平面参数优化，得到优化后上下超平面参数；
得到完整的模型前件参数及后件参数，利用(input,output)数据对计算当前IMF计算模型预测输出；
对于每一个本征模态函数IMF建立各自的T-S模糊预测模型，并计算相应IMF计算模型预测输出；进一步，反归一化后合并得到(input,output)数据对的模型预测输出。
2.根据权利要求1所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，通过变分模态分解VMD将原始风速观测数据分解为K个本征模态函数IMF，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵，具体包括：对每个IMF都执行下列内容：
对IMF前L间隔数据分量划分候选输入属性矩阵In，IOnu(tt)p＝ut[(dt)at＝a(dta),tda(atta+(tL)+,t1＝).L..Ld-aLta(t+L-1)]，其中Output为风速实际观测数据矩阵，LL为风速观测时间序列长度，利用克莱姆-施密特正交化方法GSO进行候选输入属性对实际输出值的相关性排序，从中挑选前M个属性，建立风速预测T-S模糊模型的数据对(Input,Output)；其中，M≤L。
3.根据权利要求3所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，方法还包括：
在所述方法进行模型精确度的测试时，所述数据对(Input,Output)被分为三部分，包括：第一部分数据用作模型初始化(xini,yini)，第二部分数据用作模型优化(xopt,yopt)，第三部分数据用作模型测试(xtest,ytest)；
在所述方法进行实际模型使用时，所述数据对(Input,Output)被分为两部分，包括第一部分数据用作模型初始化(xini,yini)，第二部分数据用作模型优化(xopt,yopt)。
4.根据权利要求1所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，所述利用区间二型模糊C回归聚类IT2-FCR进行T-S模糊模型参数辨识，得到上下超平面参数和α，具体包括：
步骤3.1、设置算法参数，所述算法参数包括两个模糊因子m1和m2、聚类总数C、距离调节参数η、最大迭代次数Tmax和停止迭代阈值ε；
利用模糊C均值进行粗聚类获取一型隶属度矩阵U＝[uik]，令初始和u ik是xk分别隶属于第i类上超平面和下超平面的隶属度，初始划分矩阵令当前迭代次数t＝0；
步骤3.2、利用下式计算每个聚类的上下超平面参数：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>=</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <msup> <mi>X</mi> <mi>T</mi> </msup> <msup> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mi>X</mi> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>X</mi> <mi>T</mi> </msup> <msup> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mi>y</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msub> <mo>=</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <msup> <mi>X</mi> <mi>T</mi> </msup> <msub> <munder> <mi>P</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msub> <mi>X</mi> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>X</mi> <mi>T</mi> </msup> <msub> <munder> <mi>P</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msub> <mi>y</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，X＝[x 1]∈RN×(M+1)；
步骤3.3、利用类型还原并计算实际观测值与模糊模型预测输出的误差Eik(αi)＝yk-fi(xk,αi)，其中yk是实际观测矩阵y＝(y1,y2…yN)T的第k分量；
步骤3.4、按照下式计算上下隶属度函数：
<mrow> <msub> <mover> <mi>u</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mn>2</mn> <mo>/</mo> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mn>1</mn> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mfrac> </mtd> <mtd> <mrow> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> <mo>&lt;</mo> <mfrac> <mn>1</mn> <mi>C</mi> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mn>2</mn> <mo>/</mo> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mn>2</mn> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mfrac> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
<mrow> <msub> <munder> <mi>u</mi> <mo>&amp;OverBar;</mo> </munder> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mn>2</mn> <mo>/</mo> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mn>1</mn> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mfrac> </mtd> <mtd> <mrow> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> <mo>&amp;GreaterEqual;</mo> <mfrac> <mn>1</mn> <mi>C</mi> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mfrac> <mn>1</mn> <mrow> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>E</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>E</mi> <mrow> <mi>r</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mn>2</mn> <mo>/</mo> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mn>2</mn> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mfrac> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
步骤3.5、利用上下隶属度更新计算划分矩阵和P；
步骤3.6、令t＝t+1，转至步骤3.2，直到相邻两次迭代的差值||αi current-αi last||＞ε，或者t＞Tmax时，停止迭代；其中，αi current是当前迭代的第i类超平面参数，αi last是上一次迭代的第i类超平面参数；
步骤3.7、计算每个聚类的上下超平面参数。
5.根据权利要求1所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，所述采用引力搜索算法GSA进行模型超平面参数优化，具体包括：
步骤4.1、算法初始化，设置GSA算法参数，包括最大迭代次数Nitmax、群体规模N0、初始引力常数G0，衰减系数β；
步骤4.2、计算群体目标函数最小值，具有最小目标函数值的个体确定为当前最优个体Xbest；
步骤4.3、更新所有个体的引力常数Gi；
<mrow> <msub> <mi>G</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>G</mi> <mn>0</mn> </msub> <mo>&amp;CenterDot;</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mi>&amp;beta;</mi> <mo>&amp;CenterDot;</mo> <mfrac> <mi>t</mi> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mi>t</mi> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> </mfrac> <mo>)</mo> </mrow> </mrow>
其中，G0为引力常数初始值，β为衰减系数，t为当前迭代次数，Nitmax为最大迭代次数；
步骤4.4、计算所有个体的万有引力Fi和加速度ai；
步骤4.5、更新所有个体的速度vi和位置Xi；
步骤4.6、t＝t+1，如果t>Nitmax，当前最优个体位置作为终解；否则，转入步骤4.2；
步骤4.7、输出最优个体位置即为模型上下超平面的最优参数、优化后的最优模型后件参数Α和模糊权重ωi。
6.根据权利要求5所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，所述计算群体目标函数最小值，具体包括：
步骤4.2.1、个体i位置向量Xi(t)解码得到控制参数，分别为模型上超平面参数和下超平面参数，基于超平面型隶属函数计算高斯模糊权重ωi；
步骤4.2.2、最小二乘法辨识模糊后件参数：利用数据对(Input,Output)的第一部分数据得到Α＝(ΦTΦ)-1XTy；
步骤4.2.3、分别以数据对(Input,Output)的第一部分数据和数据对(Input,Output)的第二部分数据计算模型预测输出和Φopt是按照第二部分数据计算得到；反归一化后分别计算模型预测值与实际观测值的均方根误差RMSEtrain和RMSEopt:
<mrow> <mi>R</mi> <mi>M</mi> <mi>S</mi> <mi>E</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mi>N</mi> </mfrac> </msqrt> </mrow>
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值；
以两部分数据所计算得到的RMSE的加权平均作为个体i的目标函数值：
fiti＝t1*RMSEtrain+t2*RMSEopt
其中，t1+t2＝1。
7.根据权利要求6所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，所述基于超平面型隶属函数计算高斯模糊权重ωi，具体包括：
步骤4.2.1.1、计算输入数据点距离上下超平面的距离：
<mrow> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>&amp;CenterDot;</mo> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> </mrow>
<mrow> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>&amp;CenterDot;</mo> <msup> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msup> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msup> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msup> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> </mrow>
步骤4.2.1.2、依据下式计算数据点分别隶属于上下超平面型的高斯隶属度函数和
<mrow> <msub> <mover> <mi>&amp;mu;</mi> <mo>&amp;OverBar;</mo> </mover> <msub> <mi>A</mi> <mi>i</mi> </msub> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mi>&amp;eta;</mi> <mo>&amp;CenterDot;</mo> <mfrac> <mrow> <msub> <mi>d</mi> <mi>ik</mi> </msub> <mrow> <mo>(</mo> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mi>max</mi> <mo>{</mo> <msub> <mi>d</mi> <mi>ik</mi> </msub> <mrow> <mo>(</mo> <msup> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>}</mo> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>
<mrow> <msub> <munder> <mi>&amp;mu;</mi> <mo>&amp;OverBar;</mo> </munder> <msub> <mi>A</mi> <mi>i</mi> </msub> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mi>&amp;eta;</mi> <mo>&amp;CenterDot;</mo> <mfrac> <mrow> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <munder> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </munder> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>}</mo> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>
步骤4.2.1.3、计算模糊权重ωi(xk)表示第k个输入分量属于第i类的程度。
8.根据权利要求1所述的基于区间二型T-S模糊模型的短期风速预测方法，其特征在于，对所述方法进行模型精确度的测试时，所述方法还包括：
以模糊模型预测输出与实际观测数据均方根误差RMSE，平均绝对误差MAE以及平均绝对百分误差MAPE为指标，来评估模糊模型预测的准确度，其中，均方根误差RMSE，平均绝对误差MAE和平均绝对百分误差MAPE的表达式如下：
<mrow> <mi>R</mi> <mi>M</mi> <mi>S</mi> <mi>E</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mi>N</mi> </mfrac> </msqrt> </mrow>
<mrow> <mi>M</mi> <mi>A</mi> <mi>E</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mo>|</mo> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>|</mo> </mrow>
<mrow> <mi>M</mi> <mi>A</mi> <mi>P</mi> <mi>E</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mo>|</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> </mrow> <msub> <mi>y</mi> <mi>k</mi> </msub> </mfrac> <mo>|</mo> <mo>&amp;times;</mo> <mn>100</mn> </mrow>
其中N为采样总数，yk表示k时刻实际的观测数据，表示k时刻的预测结果值；
根据所述RMSE、MAE和MAPE的结果数据，评判当前模糊模型的准确度。
9.一种基于区间二型T-S模糊模型的短期风速预测系统，其特征在于，系统包括风速检测器和服务器，所述风速检测器用于将风速观测数据data发送给服务器，具体的：
所述服务器，用于风速观测数据data预处理，通过变分模态分解VMD将原始风速观测数据分解为K个本征模态函数IMF，并对于每一个本征模态函数IMF建立各自的T-S模糊预测模型的输入输出矩阵(input,output)；
对应于每一个本征模态函数IMF，建立区间二型T-S模糊模型，其中，T-S模糊模型是由一组“IF-THEN”模糊规则来描述的非线性系统，每一个规则代表一个子系统，整个非线性系统为各个子系统的线性组合；所述IF-THEN模糊规则如下：
规则i：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msup> <mi>f</mi> <mi>i</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>,</mo> <msup> <mi>&amp;alpha;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>a</mi> <mn>1</mn> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>a</mi> <mn>2</mn> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mn>2</mn> </mrow> </msub> <mo>+</mo> <mn>...</mn> <mo>+</mo> <msubsup> <mi>a</mi> <mi>M</mi> <mi>i</mi> </msubsup> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mi>M</mi> </mrow> </msub> <mo>+</mo> <msubsup> <mi>a</mi> <mn>0</mn> <mi>i</mi> </msubsup> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>k</mi> </msub> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;alpha;</mi> <mrow> <mi>i</mi> <mi>T</mi> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>i</mi> <mo>=</mo> <mn>1...</mn> <mi>C</mi> <mo>,</mo> <mi>k</mi> <mo>=</mo> <mn>1...</mn> <mi>N</mi> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1...</mn> <mi>M</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中是线性隶属度函数的模糊集合、C是聚类规则数、N是输入输出矩阵行数、M为输入矩阵的列数、xk＝[xk1,···xkM]是模糊模型输入矩阵x＝(x1,x2…xN)T的第k分量、是第k个模型预测输出中隶属于第i类的分量、αi是模型后件参数，用于指代超平面；
再通过线性组合得到T-S模糊模型；其中，
记 ωi为高斯模糊权重；
根据(input,output)数据对，通过区间二型模糊C回归聚类IT2-FCR，进行T-S模糊模型参数辨识，得到上下超平面参数和α；
根据(input,output)数据对，利用引力搜索算法GSA进行模型上下超平面参数优化，得到优化后上下超平面参数；
得到完整的模型前件参数及后件参数，利用(input,output)数据对计算当前IMF计算模型预测输出；
对于每一个本征模态函数IMF建立各自的T-S模糊预测模型，并计算相应IMF计算模型预测输出；进一步，反归一化后合并得到(input,output)数据对的模型预测输出。
CN201711463836.5A 2017-12-28 2017-12-28 一种基于区间二型t-s模糊模型的短期风速预测方法与系统 Pending CN108053077A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711463836.5A CN108053077A (zh) 2017-12-28 2017-12-28 一种基于区间二型t-s模糊模型的短期风速预测方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711463836.5A CN108053077A (zh) 2017-12-28 2017-12-28 一种基于区间二型t-s模糊模型的短期风速预测方法与系统

Publications (1)

Publication Number Publication Date
CN108053077A true CN108053077A (zh) 2018-05-18

Family

ID=62128594

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711463836.5A Pending CN108053077A (zh) 2017-12-28 2017-12-28 一种基于区间二型t-s模糊模型的短期风速预测方法与系统

Country Status (1)

Country Link
CN (1) CN108053077A (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460608A (zh) * 2018-11-09 2019-03-12 东北大学 一种基于模糊时间序列的高陡边坡形变预测的方法
CN109934386A (zh) * 2019-01-30 2019-06-25 山西河坡发电有限责任公司 热电联供系统热负荷预测方法
CN110058526A (zh) * 2019-05-20 2019-07-26 杭州电子科技大学 一种基于区间二型t-s模型的中立型系统的控制方法
CN111581746A (zh) * 2020-05-11 2020-08-25 中国矿业大学 一种新型三相圆筒开关磁阻直线发电机多目标优化方法
CN111898628A (zh) * 2020-06-01 2020-11-06 淮阴工学院 一种新型t-s模糊模型辨识方法
CN112465034A (zh) * 2020-11-30 2021-03-09 中国长江电力股份有限公司 一种基于水轮发电机的t-s模糊模型的建立方法及系统
CN113156811A (zh) * 2020-12-29 2021-07-23 上海理工大学 一种航天器姿态控制系统的自触发模糊模型预测控制器的设计方法
CN113515048A (zh) * 2021-08-13 2021-10-19 华中科技大学 一种模糊自适应pso-elm声品质预测模型的建立方法
CN113608435A (zh) * 2021-06-29 2021-11-05 河南科技大学 抑制高速飞行器散热片表面温度的弹性控制器设计方法
CN113625558A (zh) * 2021-07-23 2021-11-09 杭州电子科技大学 基于区间二型t-s模糊的欺骗攻击下的网络控制方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103400052A (zh) * 2013-08-22 2013-11-20 武汉大学 一种风电场短期风速组合预测方法
CN104533701A (zh) * 2014-12-23 2015-04-22 华中科技大学 一种水轮机调速系统控制参数的自动整定方法
CN105391083A (zh) * 2015-11-05 2016-03-09 河海大学 基于变分模态分解和相关向量机的风功率区间短期预测方法
CN106374465A (zh) * 2016-11-10 2017-02-01 南京信息工程大学 基于gsa‑lssvm模型的短期风电功率预测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103400052A (zh) * 2013-08-22 2013-11-20 武汉大学 一种风电场短期风速组合预测方法
CN104533701A (zh) * 2014-12-23 2015-04-22 华中科技大学 一种水轮机调速系统控制参数的自动整定方法
CN105391083A (zh) * 2015-11-05 2016-03-09 河海大学 基于变分模态分解和相关向量机的风功率区间短期预测方法
CN106374465A (zh) * 2016-11-10 2017-02-01 南京信息工程大学 基于gsa‑lssvm模型的短期风电功率预测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
CHAOSHUN LI等: "《Hydraulic turbine governing system identification using T–S fuzzy model optimized by chaotic gravitational search algorithm》", 《ENGINEERING APPLICATIONS OF ARTIFICIAL INTELLIGENCE》 *
WEN ZOU等: "《A T-S Fuzzy Model Identification Approach Based on a Modified Inter Type-2 FRCM Algorithm》", 《IEEE TRANSACTIONS ON FUZZY SYSTEMS》 *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460608A (zh) * 2018-11-09 2019-03-12 东北大学 一种基于模糊时间序列的高陡边坡形变预测的方法
CN109460608B (zh) * 2018-11-09 2022-11-25 东北大学 一种基于模糊时间序列的高陡边坡形变预测的方法
CN109934386A (zh) * 2019-01-30 2019-06-25 山西河坡发电有限责任公司 热电联供系统热负荷预测方法
CN109934386B (zh) * 2019-01-30 2022-06-28 山西河坡发电有限责任公司 热电联供系统热负荷预测方法
CN110058526B (zh) * 2019-05-20 2021-11-30 杭州电子科技大学 一种基于区间二型t-s模型的中立型系统的控制方法
CN110058526A (zh) * 2019-05-20 2019-07-26 杭州电子科技大学 一种基于区间二型t-s模型的中立型系统的控制方法
CN111581746A (zh) * 2020-05-11 2020-08-25 中国矿业大学 一种新型三相圆筒开关磁阻直线发电机多目标优化方法
CN111898628A (zh) * 2020-06-01 2020-11-06 淮阴工学院 一种新型t-s模糊模型辨识方法
CN111898628B (zh) * 2020-06-01 2023-10-03 淮阴工学院 一种新型t-s模糊模型辨识方法
CN112465034A (zh) * 2020-11-30 2021-03-09 中国长江电力股份有限公司 一种基于水轮发电机的t-s模糊模型的建立方法及系统
CN112465034B (zh) * 2020-11-30 2023-06-27 中国长江电力股份有限公司 一种基于水轮发电机的t-s模糊模型的建立方法及系统
CN113156811A (zh) * 2020-12-29 2021-07-23 上海理工大学 一种航天器姿态控制系统的自触发模糊模型预测控制器的设计方法
CN113608435A (zh) * 2021-06-29 2021-11-05 河南科技大学 抑制高速飞行器散热片表面温度的弹性控制器设计方法
CN113608435B (zh) * 2021-06-29 2024-05-10 河南科技大学 抑制高速飞行器散热片表面温度的弹性控制器设计方法
CN113625558A (zh) * 2021-07-23 2021-11-09 杭州电子科技大学 基于区间二型t-s模糊的欺骗攻击下的网络控制方法
CN113625558B (zh) * 2021-07-23 2024-09-06 杭州电子科技大学 基于区间二型t-s模糊的欺骗攻击下的网络控制方法
CN113515048A (zh) * 2021-08-13 2021-10-19 华中科技大学 一种模糊自适应pso-elm声品质预测模型的建立方法
CN113515048B (zh) * 2021-08-13 2023-04-07 华中科技大学 一种模糊自适应pso-elm声品质预测模型的建立方法

Similar Documents

Publication Publication Date Title
CN108053077A (zh) 一种基于区间二型t-s模糊模型的短期风速预测方法与系统
CN104881706B (zh) 一种基于大数据技术的电力系统短期负荷预测方法
Li et al. A light gradient boosting machine for remainning useful life estimation of aircraft engines
Wang et al. Short-term wind power prediction based on improved small-world neural network
CN106447098B (zh) 一种光伏超短期功率预测方法及装置
CN109146162B (zh) 一种基于集成循环神经网络的概率风速预测方法
Hoori et al. Electric load forecasting model using a multicolumn deep neural networks
Gong et al. Bottom-up load forecasting with Markov-based error reduction method for aggregated domestic electric water heaters
CN110942205A (zh) 一种基于himvo-svm的短期光伏发电功率预测方法
CN115995810A (zh) 一种考虑气象波动自适应匹配的风电功率预测方法及系统
Guoqiang et al. Study of RBF neural network based on PSO algorithm in nonlinear system identification
CN111008790A (zh) 一种水电站群发电调度规则提取方法
CN110738363B (zh) 一种光伏发电功率预测方法
CN117081063A (zh) 一种基于GCN-Crossformer模型的分布式充电负荷预测方法及系统
CN111027760A (zh) 一种基于最小二乘向量机的电力负荷预测方法
CN116842848B (zh) 基于混合搜索进化优化的工业过程软测量建模方法、系统
Alam et al. AI-based efficiency analysis technique for photovoltaic renewable energy system
Hu et al. Scenario reduction based on correlation sensitivity and its application in microgrid optimization
CN112330051A (zh) 一种基于Kmeans与FR-DBN的短期负荷预测方法
CN105956720A (zh) 一种基于t-s模糊模型的风电功率超短期预测方法
Sun et al. Short-term photovoltaic power prediction modeling based on AdaBoost algorithm and Elman
Fu et al. A spatial forecasting method for photovoltaic power generation combined of improved similar historical days and dynamic weights allocation
CN110059871A (zh) 光伏发电功率预测方法
Han et al. Swarm intelligent algorithm for re-entrant hybrid flow shop scheduling problems
Ogawa et al. PV output forecasting by deep Boltzmann machines with SS‐PPBSO

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
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20180518

WD01 Invention patent application deemed withdrawn after publication