CN108733003B - 基于kmeans聚类算法的回转类零件工序工时预测方法及系统 - Google Patents

基于kmeans聚类算法的回转类零件工序工时预测方法及系统 Download PDF

Info

Publication number
CN108733003B
CN108733003B CN201710262953.9A CN201710262953A CN108733003B CN 108733003 B CN108733003 B CN 108733003B CN 201710262953 A CN201710262953 A CN 201710262953A CN 108733003 B CN108733003 B CN 108733003B
Authority
CN
China
Prior art keywords
historical
clustering
matrix
predicted
working
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
CN201710262953.9A
Other languages
English (en)
Other versions
CN108733003A (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.)
Yangzhou Jiangdu Yongjian Co ltd
Nanjing University of Science and Technology
Original Assignee
Yangzhou Jiangdu Yongjian Co ltd
Nanjing 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 Yangzhou Jiangdu Yongjian Co ltd, Nanjing University of Science and Technology filed Critical Yangzhou Jiangdu Yongjian Co ltd
Priority to CN201710262953.9A priority Critical patent/CN108733003B/zh
Publication of CN108733003A publication Critical patent/CN108733003A/zh
Application granted granted Critical
Publication of CN108733003B publication Critical patent/CN108733003B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B19/00—Programme-control systems
    • G05B19/02—Programme-control systems electric
    • G05B19/418—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM]
    • G05B19/41865—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM] characterised by job scheduling, process planning, material flow
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B2219/00—Program-control systems
    • G05B2219/30—Nc systems
    • G05B2219/32—Operator till task planning
    • G05B2219/32252—Scheduling production, machining, job shop
    • 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/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]

Landscapes

  • Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Probability & Statistics with Applications (AREA)
  • Manufacturing & Machinery (AREA)
  • Quality & Reliability (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于kmeans聚类算法的回转类零件工序工时预测方法及系统，工序数据采集与处理模块用于对历史工序的数据采集与处理；订单导入模块用于将该企业的零件订单导入该回转类零件工序工时预测系统；工序编码模块用于对历史工序及待预测工序进行工序编码；聚类模块用于对历史零件工序编码进行定期聚类成组；工时预测模块对待预测工序进行工时预测，其中该模块包含的聚类组划分单元用于将待预测工序划分至与其最相近的聚类组，工时计算单元用于对工时进行预测计算；数据存储模块用于对该系统中的数据进行存储；该方法与系统具有对回转类零件工序工时预测快速、准确、少数据量时预测效果较好的优点。

Description

基于kmeans聚类算法的回转类零件工序工时预测方法及系统
技术领域
本发明属于离散制造领域，特别是一种基于聚类kmeans算法的回转类零件工序工时预测方法及系统。
背景技术
液压基础零部件大多为回转类零件，其按订单生产的特点决定了其生产类型为多品种小批量，其生产企业为典型的离散加工企业。多品种小批量的生产方式导致企业管理难度的加大，制造过程中工时数据不及时、不准确，是导致其制造执行系统中生产计划编排产生偏差的主要原因。柔性车间作业调度是制造执行系统中的一个重要问题，其准确实施的基础是能够快速、准确地获得工序工时。
目前，企业对于工时的预测主要采用两种方式：(1)程序仿真法(2)神经网络法。程序仿真法必须在零件三维模型的基础上才能对三维模型进行自动编程，仿真实际加工工时，而车间对于大部分图纸的管理仅限于CAD二维模型，因此该方法实施过程中耗时耗力；神经网络法主要是通过从大量历史数据中学习、训练，得到输出变量即工时与输入变量之间的关系，在该方法中，训练过程需要大量的历史数据，大多数企业对工时的管理还不够完善，其数据量也不足以满足神经网络训练的要求。
发明内容
本发明为解决现有回转类零件工时预测响应慢且数据需求量大的问题，本发明提出了一种基于kmeans聚类算法的回转类零件工时预测方法及系统。
实现本发明目的的技术解决方案为：
一种基于kmeans聚类算法的回转类零件工序工时预测方法，包括以下步骤：
步骤1、采集回转类零件工序历史数据：从企业MES、ERP、PLM系统中对回转类零件工序历史数据进行采集或计算；
步骤2、根据工艺信息及零件图纸信息将历史零件工序进行编码，将各历史零件工序编码转换为历史工序编码矩阵Xn×m：
整理出企业已加工回转类零件工序的工艺信息及零件图纸信息，对所有回转类零件工艺信息及图纸中所出现的信息及结构进行标识；对不同工序的零件图纸信息和工艺信息进行组合，构成历史零件工序编码；各个历史零件工序编码完成后，将各历史零件工序编码转换为历史工序编码矩阵Xn×m；
步骤3、对历史工序编码矩阵依据kmeans聚类算法进行聚类：采用初始聚类中心计算方法、距离计算方法及更新聚类中心方法，得到不同K值下的聚类结果及聚类中心，计算不同K值下的SP/DB指标，取SP/DB取最小时对应的K值、聚类中心及聚类结果，即为最终K值、最终聚类结果及最终聚类中心；
步骤4、订单拆分成工序，对待预测工序进行编码：将订单拆分成各个工序并判断各工序是否为外协加工，如果不是外协加工则认定为待预测工序，需进行工时预测，预测时首先对待预测的工序进行编码；
步骤5、将待预测工序划分至最终聚类组：
结合权重系数矩阵公式、加权欧式距离公式计算所有待预测工序到K个聚类中心的加权欧式距离，相互比较后将所有待预测工序划分至与自身距离最小的最终聚类组，并筛选出步骤1中采集的所有数据中属于该聚类组的历史零件图号、工艺信息、零件工序历史特征参数矩阵及零件工序历史工时数据；
步骤6、根据零件图号及工艺信息，将待预测工序与历史工序进行对比，若预测工序与历史工序的零件图号及工艺信息均相同，直接将历史理论标准工时作为预测理论工时并输出预测结果；若预测工序与历史工序不相同，则跳过该步骤直接进行下一步骤；
步骤7、得到历史工时z与PCA降至二维后的特征参数(x，y)的对应关系：用MatLab画出待预测工序所属聚类组中历史工时与特征参数的三维散点图，对三维散点图运用MatLab中用最小二乘法拟合生成相应的函数关系式z＝f(x，y)；
步骤8、预测工时计算：以新矩阵中s个待预测数据的第一维为x坐标，新矩阵中s个待预测数据的第二维为y坐标，代入关系式z＝f(x，y)得到待预测工序的理论预测工时公式并输出预测结果。
一种回转类零件工序工时预测系统，包括工序数据采集与处理模块、订单导入模块、工序编码模块、聚类模块、工时预测模块、数据存储模块；
所述工序数据采集与处理模块用于对历史工序的数据采集与处理；
所述订单导入模块用于将该企业的零件订单导入该回转类零件工序工时预测系统；
所述工序编码模块用于对历史工序及待预测工序进行工序编码，同时输出历史工序编码结果及待预测工序编码结果；
所述聚类模块，用于对历史零件工序编码进行定期聚类成组；
所述聚类单元主要用于对历史工序编码进行kmeans聚类算法将历史工序编码聚类成组；
所述工时预测模块用于对待预测工序进行工时预测；
所述数据存储模块主要用以存储所述回转类工序工时预测系统的所有关键数据。
本发明与现有技术相比，其显著优点为：
(1)聚类中心更新过程中，计算潜在孤立点的接受概率，并采用随机判断其接受与否，与传统kmeans聚类算法中直接采用均值更新聚类中心相比，降低了孤立点的不良影响，同时以随机方式对潜在孤立点进行选择与否，能够以扰动的形式避免其过早陷入局部极小值；
(2)距离计算中引入权重系数，相比于原有kmeans聚类算法能有效对各属性重要程度进行评定；
(3)基于待预测工序所属的聚类组分析历史相似工序，减少计算量提高相应速度的同时提高了预测的准确性与可靠性。
附图说明
图1为本发明基于kmeans聚类算法的回转类零件工序工时预测流程图；
图2为本发明基于kmeans聚类算法的回转类零件工序工时预测整体模块图；
图3为本发明步骤3所述kmeans聚类算法与传统kmeans聚类算法在相同实施方式下的收敛曲线图。
具体实施方式
下面结合附图对本发明作进一步详细描述：
结合图1，一种基于kmeans聚类算法的回转类零件工序工时预测方法,包括以下步骤：
步骤1、采集回转类零件工序历史数据：从企业MES、ERP、PLM系统中对回转类零件工序历史数据进行采集或计算；
回转类零件工序历史数据包括零件图号、工艺信息、零件工序历史工时及历史加工设备、设备效率系数、零件工序历史特征参数矩阵；
1.1对零件图号、工艺信息、零件工序历史工时及历史加工设备的采集及设备效率系数的计算：
零件图号、工艺信息、零件工序历史工时及历史加工设备直接从企业MES系统中读取；设备效率系数由MES系统中的工序实际加工时间经计算得到，计算方式为该设备上各工序的历史标准工时与设备实际加工时间商的平均值；
1.2对零件工序历史特征参数矩阵采集：
对ERP系统中的工艺流程及PLM系统中不同零件的CAD图纸进行统计，统计工艺流程中每个零件中各道工序对应CAD图纸的零件加工特征，如余量、长度、外径、外表面台阶个数、内外表面槽个数、内外螺纹个数及长度等零件工序历史特征参数，并将每个统计量作为一个维度从前往后依次组成零件工序历史特征参数矩阵；
步骤2、根据工艺信息及零件图纸信息将历史零件工序进行编码，将各历史零件工序编码转换为历史工序编码矩阵Xn×m：
整理出企业已加工回转类零件工序的工艺信息及零件图纸信息，对所有回转类零件工艺信息及图纸中所出现的信息及结构用数字0、1、2……进行标识；对不同工序的零件图纸信息和工艺信息进行组合，构成历史零件工序编码；各个历史零件工序编码完成后，将各历史零件工序编码转换为历史工序编码矩阵Xn×m。
2.1对外圆特征信息进行标识：
如对外圆特征信息进行标识时，如表1所示，具体步骤为：(1)PLM系统中导出所有回转类零件的CAD图纸；(2)对所有图纸中外圆表面出现的特征进行分类：如光滑外圆面、单向台阶外圆、双向台阶外圆、单向台阶外圆带功能槽、单向台阶外圆带螺纹、单向台阶外圆带表面孔、单向台阶外圆带功能槽及螺纹、单向台阶外圆带功能槽及表面孔、单向台阶外圆带螺纹及表面孔、双向台阶外圆带功能槽、双向台阶外圆带螺纹、双向台阶外圆带表面孔、双向台阶外圆带功能槽及螺纹、双向台阶外圆带功能槽及表面孔、双向台阶外圆带螺纹及表面孔；(3)用数字对描述的所有特征进行标识：如(2)中所分类的外圆特征顺序对应标识为0、1、2、3、4、5、6、7、8、9、10、11…，如表1所示；
表1外圆特征对应数字标识
Figure GDA0002585601030000041
Figure GDA0002585601030000051
2.2对其他特征进行标识，其中其他特征包括零件类别特征信息、工艺信息、材料种类、最大外径、长度、内孔特征信息以及端面特征信息等：
标识方法与上述外圆表面特征标识方法类似，对回转类零件所有图纸中零件类别特征信息、工艺信息、材料种类、最大外径、长度、内孔特征信息、端面特征信息等出现的特征进行如下表2所示的标识；
表2其他特征对应数字标识
Figure GDA0002585601030000052
Figure GDA0002585601030000061
上述表1、表2的实施例仅为了更清楚的说明本发明的实施方式，表中的外圆特征及其他特征可根据实际需要进行调整。
2.3对不同工序的零件图纸信息和工艺信息进行组合，构成历史零件工序编码：
所有工序编码信息的特征(共m个)均用数字标识完成后，根据不同工序的图纸信息与工艺信息，将该工序划分至特定的组别所对应的数字，并按照零件类别特征信息、工艺信息、材料种类、最大外径、长度、外圆特征信息、内孔特征信息以及端面特征信息等的顺序将数字进行组合，构成历史零件工序编码；
例如,若历史零件工序编码为0 0 1 7 4 12 4 0…表示该零件特征为：零件类别特征为：导向套、工艺信息为：精车、材料种类为：中碳钢、最大外径范围为：
Figure GDA0002585601030000062
长度范围为：(200，240]、外圆特征信息为：双向台阶外圆带功能槽及螺纹、内孔特征信息为：单向台阶内孔带功能槽、端面特征信息为：光滑端面，……；
本实施例中采用表3的方式对零件图纸信息和工艺信息进行组合；在一些实施方式中，组合方式可根据需要做调整。
表3历史零件工序编码组合方式
Figure GDA0002585601030000071
2.4各个历史零件工序编码完成后，将各历史零件工序编码转换为历史工序编码矩阵Xn×m，矩阵中每一行代表一个历史零件工序样本，每一列代表一个零件特征。
步骤3、对历史工序编码矩阵依据kmeans聚类算法进行聚类：采用初始聚类中心计算方法、距离计算方法及更新聚类中心方法，得到不同K值下的聚类结果及聚类中心，计算不同K值下的SP/DB指标，取SP/DB取最小时对应的K值、聚类中心及聚类结果，即为最终K值、最终聚类结果及最终聚类中心；
3.1确定聚类组数K的取值范围是
Figure GDA0002585601030000072
随机选择一个K值；
3.2为满足初始聚类中心尽可能距离较远的要求，由以下公式确定初始聚类中心：
(1)根据历史工序编码矩阵Xn×m求得历史工序样本均值矩阵A1×m和历史工序样本极差矩阵B1×m：
Figure GDA0002585601030000073
B1×m＝{bj|bj＝xjmax-xjmin，j＝1，2，…，m} (2)
(2)根据历史工序样本均值矩阵A1×m和历史工序样本极差矩阵B1×m，计算聚类中心矩阵C：
Figure GDA0002585601030000074
式中：xij-历史工序编码矩阵Xn×m中第i个历史工序样本第j个特征的数值；
xjmax-历史工序编码矩阵Xn×m中第j个特征的最大值；
xjmin-历史工序编码矩阵Xn×m中第j特征的最小值；
aj-历史工序编码矩阵Xn×m中所有历史工序样本第j个特征的均值；
bj-历史工序编码矩阵Xn×m中第j个特征的极差值；
crj-第r个聚类中心的第j个特征数值；
RK×m-随机正太分布矩阵；
C-聚类中心矩阵；
3.3以K个聚类中心为基准，计算所有历史工序样本到K个聚类中心的距离，相互比较后将历史工序样本划分至距离最小的那一组，对所有历史工序样本划分完成后，得到历史工序样本的一次聚类结果；
在计算距离时，引入权重系数对历史工序样本各零件特征进行重要程度划分，权重系数的计算采用层次分析法构造判断矩阵，经归一化处理及一致性检验后得到其最终的权重系数矩阵S：
S＝{sj，j＝1，2，…，m} (4)
Figure GDA0002585601030000081
式中：sj-权重系数矩阵中第j列的数值；
dir-第i个样本到第r个聚类中心的加权欧式距离；
3.4根据3.3得到的聚类结果，更新原有聚类中心，得到K个新的聚类中心：
为降低历史工序样本中孤立点的负面影响，同时改善传统kmeans算法易陷入局部极小值的弊端，针对每个聚类组，分别采用下述步骤更新原有聚类中心：
3.4.1计算该聚类组下各历史工序样本两两之间的加权欧式距离diw,Ω为属于该聚类组的所有历史工序样本按Xn×m行递增的方式构成的矩阵，如若聚类组1中的数据有nΩ条，则聚类组1下的两两加权欧式距离应该有
Figure GDA0002585601030000082
个；
Figure GDA0002585601030000083
式中，diw-第i个历史工序样本到第w个历史工序样本的加权欧式距离；
xwj-历史工序编码矩阵Xn×m中第w个历史工序样本第j个特征的数值；
3.4.2求出3.4.1中所有加权欧式距离diw的平均值，作为平均距离dis_avrg；
Figure GDA0002585601030000091
3.4.3对该聚类组下的历史工序样本进行潜在孤立点的判断，得到密集点历史工序样本矩阵Φ及潜在孤立点历史工序样本矩阵Ψ：
判断规则为：对3.4.1中计算的所有加权欧式距离diw与平均距离dis_avrg比较，如若3.4.1中计算的diw＜dis_avrg，则认为历史工序样本i、w均为密集点样本；所有加权欧式距离均比较完成后，将得到的所有密集点样本去除重复项并取出密集点历史工序样本的数据构成矩阵Φ；由Φ∪Ψ＝Ω计算得到矩阵Ψ，矩阵Ψ即为所有潜在孤立点历史工序样本；求出矩阵Φ中各特征的平均值作为基础聚类中心cetr，同时统计矩阵Φ与矩阵Ψ中历史工序样本数量分别为nΦ，nΨ；
3.4.4对于潜在孤立点历史工序样本矩阵Ψ中的历史工序样本β，以一定的概率Pβ接受该样本数据β的存在；同时采用随机方式对潜在孤立点历史工序样本矩阵Ψ中的样本的数据进行接受与舍弃，构成接受孤立点样本历史工序样本矩阵F：
Figure GDA0002585601030000092
式中：dβcetr-潜在孤立点历史工序样本矩阵Ψ中的β样本与基础聚类中心cetr的加权欧式距离；
MatLab中采用randn函数产成nΨ个标准正太分布的伪随机数，与3.4.4中的Pβ进行对应比较，若小于Pβ，则接受该潜在孤立点历史工序样本，若大于Pβ，则舍弃该潜在孤立点历史工序样本；统计出所有被接受的潜在孤立点历史工序样本，构成接受孤立点样本历史工序样本矩阵
Figure GDA0002585601030000093
接受孤立点样本历史工序样本矩阵F中历史工序样本数量为nF；
3.4.5由下式计算新聚类中心：
Figure GDA0002585601030000094
3.5对新聚类中心，按照步骤3.3和3.4进行迭代，直到聚类中心不再变化则完成聚类，得到该K值下聚类结果及聚类中心；注意在迭代的过程中权重系数矩阵S不发生变化。
3.6根据3.1步骤确定的K值范围，选择不同的K值，按照步骤3.2-3.5迭代得到不同K值下的聚类结果及聚类中心。
3.7计算不同聚类组数K值下的分类评价指标SP、DB，并取SP/DB取最小时对应的K值、聚类结果及聚类中心为最终K值、最终聚类结果及最终聚类中心：
(1)根据不同K值下的聚类结果及聚类中心计算分类指标SP与DB：
Figure GDA0002585601030000101
Figure GDA0002585601030000102
式中：Chs-第h组与第s组聚类中心的加权欧式距离；
Wh-第h组内所有历史样本数据到该组聚类中心的平均加权欧式距离；
(2)取SP/DB取最小时对应的K值、聚类中心及聚类结果为本发明最终K值、最终聚类结果及最终聚类中心。
图3为按步骤3所述kmeans聚类算法进行聚类计算与传统kmeans聚类算法进行聚类计算时，各历史工序样本到聚类中心的平均距离与迭代次数之间的关系。由图可以看出，本实施例所述kmeans聚类算法跟传统kmeans聚类算法比较，具有收敛效果较好不易于陷入局部极小值的优点。
步骤4、订单拆分成工序，对待预测工序进行编码：
将订单拆分成各个工序并判断各工序是否为外协加工，如果不是外协加工则认定为待预测工序，需进行工时预测，预测时首先对待预测的工序进行编码，编码过程如上述步骤2.3；
步骤5、将待预测工序划分至最终聚类组：
结合公式(4)、(5)计算所有待预测工序到K个聚类中心的加权欧式距离，相互比较后将所有待预测工序划分至与自身距离最小的最终聚类组，并筛选出步骤1中采集的所有数据中属于该聚类组的历史零件图号、工艺信息、零件工序历史特征参数矩阵及零件工序历史工时数据；
步骤6、根据零件图号及工艺信息，将待预测工序与历史工序进行对比，若预测工序与历史工序的零件图号及工艺信息均相同，直接将历史理论标准工时作为预测理论工时并输出预测结果；若预测工序与历史工序不相同，则跳过该步骤直接进行下一步骤；
历史理论标准工时为：
Figure GDA0002585601030000111
式中：tp实—设备p加工的实做工时；
αp—设备p效率系数，由步骤1.1得到；
g—零件图号及工艺信息均相同的历史工序的个数；
p—与待预测工序零件图号及工艺信息均相同的第p个历史工序；
步骤7、得到历史工时z与PCA降至二维后的特征参数(x，y)的对应关系：用MatLab画出待预测工序所属聚类组中历史工时与特征参数的三维散点图，对三维散点图运用MatLab中用最小二乘法拟合生成相应的函数关系式z＝f(x，y)；
7.1用MatLab画出待预测工序所属聚类组中历史工时与特征参数的三维散点图：
按照步骤1.2的方法统计待预测工序的工序历史特征参数矩阵，假设共有q个待预测工序，其中有s个待预测工序由步骤5被划分至同一个聚类组如聚类组1，将待预测工序的工序历史特征参数矩阵附加至聚类组1的工序历史特征参数矩阵的最后s行，构成附加矩阵，并对附加后的矩阵进行PCA降维，降至二维形成新矩阵，假设聚类组1中有l个历史样本数据，则新矩阵为(l+s)×2维；
以新矩阵中l个历史样本数据的第一维为x坐标，新矩阵中l个历史样本数据的第二维为y坐标，l个历史样本数据的历史工时为z坐标，利用Matlab画出三维散点图；
7.2对三维散点图运用MatLab中用最小二乘法拟合生成相应的函数关系式z＝f(x，y)；
步骤8、预测工时计算：
以新矩阵中s个待预测数据的第一维为x坐标，新矩阵中s个待预测数据的第二维为y坐标，代入第7.2求得的关系式得到待预测工序的理论预测工时公式并输出预测结果。
下面结合具体实施例，对一批回转类零件进行历史工序编码并聚类，聚类后在聚类结果取中取一个分组对其中的待预测工序按照上述步骤进行工时预测并输出结果，结果如下表所示，其中所有工序的特征参数中余量相同，为简化分析余量不列入表4中。
表4该聚类组下待预测工时工序与零件工序历史特征参数矩阵
Figure GDA0002585601030000112
Figure GDA0002585601030000121
MatLab中用最小二乘法拟合生成相应的函数关系式：
z＝f(x，y)＝p00+p10*x+p01*y+p20*x2+p11*x*y+p02*y2+p30*x3+p21*x2*y+p12*x*y2+p03*y3+p40*x4+p31*x3*y+p22*x2*y2+p13*x*y3
其中，p00＝12.64，p10＝0.5314，p01＝-1.052，p20＝-0.01747，p11＝-0.0424，p02＝0.286，p30＝0.0001266，p21＝-0.001025，p12＝-0.004595，p03=0.01802，p40＝2.309e-5，p31＝4.65e-5，p22＝-0.00023，p13＝-0.0001118。
代入降二维后的待预测数据(-46.3942,3.0449)与(-26.6253,8.1068)，得到待预测工时为14.8294和8.8893，经MES系统采集的实际工时乘设备系数计算得到理论工时为14.63636和9.232222，预测误差为1.32％与3.71％。采用本方法的工序工时预测具有很好的有效性。
结合图2，本发明基于上述实施方法提出了一种回转类零件工序工时预测系统，包括工序数据采集与处理模块、订单导入模块、工序编码模块、聚类模块、工时预测模块、数据存储模块；
所述工序数据采集与处理模块用于对历史工序的数据采集与处理；历史工序数据包括：零件图号、工艺信息、零件工序历史工时及历史加工设备、零件工序历史特征参数矩阵采集、设备效率系数；具体在方法的实施例中已介绍，此处不再赘述；
所述订单导入模块用于将该企业的零件订单导入该回转类零件工序工时预测系统；
所述工序编码模块用于对历史工序及待预测工序进行工序编码，同时输出历史工序编码结果及待预测工序编码结果；
其中工序编码方式为：整理出企业已加工回转类零件工序的工艺信息及零件图纸信息，对所有回转类零件工艺信息及图纸中所出现的信息及结构用数字0、1、2……进行标识；对待编码工序的m个特征分别进行编码，找出其各个特征对应的编码值，并按照零件类别特征信息、工艺信息、材料种类、最大外径、长度、外圆特征信息、内孔特征信息以及端面特征信息等的顺序将数字进行组合，构成工序编码；
所述聚类模块，用于对历史零件工序编码进行定期聚类成组，其中该模块包括定期更新单元与聚类单元；
所述定期更新单元主要用于定期更新聚类结果信息，该单元将工序数据采集模块中采集的信息定期传送至聚类模块中的聚类单元；
所述聚类单元主要用于对历史工序编码进行kmeans聚类算法将历史工序编码聚类成组；
聚类过程为：确定聚类组数K的取值范围是
Figure GDA0002585601030000131
随机选择K值；由上述方法实施例中的公式(1)、(2)、(3)确定初始聚类中心；以每个聚类中心为基准，由上述方法实施例中的公式(4)、(5)计算所有历史样本数据到K个聚类中心的加权欧式距离，相互比较后将历史样本数据划分至距离最小的那一组，对所有历史样本数据划分完成后，得到历史样本数据的一次聚类结果；由上述方法实施例中的公式(6)、(7)、(8)、(9)计算新聚类中心，更新原有聚类中心；对新聚类中心，重新进行聚类组成员的划分，直到聚类中心不再变化则完成聚类，得到聚类结果及聚类中心；选择不同K值，计算不同K值下的聚类结果及聚类中心；由上述方法中的公式(10)、(11)计算不同K值下的分类评价指标SP、DB，选择SP/DB取最小时对应的K值、聚类结果及聚类中心为最终K值、最终聚类结果及最终聚类中心；具体过程见上述方法实施例；
所述工时预测模块用于对待预测工序进行工时预测，其中该模块包括订单拆分单元、聚类组划分单元、历史匹配单元与工时计算单元；
所述订单拆分单元用于对订单导入模块导入的零件订单进行工序拆分，保留拆分的非外协加工工序为待预测工序；
所述聚类组划分单元主要用于将待预测工序划分至合适的历史聚类组中，划分过程为：对工序编码模块中输出的待预测工序编码，结合上述方法实施例中的公式(4)、(5)计算所有待预测工序到K个聚类中心的加权欧式距离，相互比较后将待预测工序编码一一划分至距离最小那一聚类组，同时结合工序编码模块及数据采集模块，输出待预测工序及其所属聚类组的工序编码、零件图号、工艺信息及所属聚类组的工序编码、零件图号、工艺信息、零件工序历史特征参数矩阵及零件工序历史工时数据信息；
所述历史匹配单元主要用于根据聚类组划分单元输出信息中的零件图号及工艺信息对待预测工序与所属聚类组的历史工序进行匹配，并将匹配结果输出至工时计算单元；
所述工时计算单元主要用于对预测待加工工时计算并输出预测结果，其预测过程为：判断历史匹配单元的输出结果，若匹配成功，则按照上述方法实施例中的公式(12)计算，并将计算结果作为理论预测工时；若匹配不成功，则统计待预测工序的工序历史特征参数矩阵，将其与聚类组划分单元输出的工序历史特征参数矩阵进行附加，并对附加矩阵进行PCA降二维处理形成新矩阵，以新矩阵中历史样本数据的第一维为x坐标，新矩阵中历史样本数据的第二维为y坐标，历史样本数据的历史工时为z坐标，画出三维散点图，对三维散点图运用MatLab中用最小二乘法拟合生成相应的函数关系式z＝f(x，y)，以新矩阵中待预测数据的第一维为x坐标，新矩阵中待预测数据的第二维为y坐标，代入关系式z＝f(x，y)得到待预测工序的理论预测工时并输出预测结果。
所述数据存储模块主要用以存储所述回转类工序工时预测系统的所有关键数据，包括：数据采集与处理模块采集与处理的所有零件图号、工艺信息、零件工序历史工时及历史加工设备、零件工序历史特征参数矩阵、设备效率系数数据，订单导入模块中导入的订单数据，工序编码模块输出的工序编码数据，聚类模块的聚类结果数据，工时预测模块的待预测工序所属聚类组数据及其预测工时数据。
以该种方式对回转类零件工序工时进行预测，采用本文所述kmeans聚类算法能够快速、简单、准确的将编码信息聚类成组，划分聚类组能够有效减少分析量提高分析速度；对组内相似成员结构分析，由于其相似性的存在，在数据量不大时也能有较好的预测效果，有效解决了回转类零件工时预测响应慢或数据需求量大的问题。

Claims (10)

1.一种基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，包括以下步骤：
步骤1、采集回转类零件工序历史数据：从企业MES、ERP、PLM系统中对回转类零件工序历史数据进行采集或计算；
步骤2、根据工艺信息及零件图纸信息将历史零件工序进行编码，将各历史零件工序编码转换为历史工序编码矩阵Xn×m：
整理出企业已加工回转类零件工序的工艺信息及零件图纸信息，对所有回转类零件工艺信息及图纸中所出现的信息及结构进行标识；对不同工序的零件图纸信息和工艺信息进行组合，构成历史零件工序编码；各个历史零件工序编码完成后，将各历史零件工序编码转换为历史工序编码矩阵Xn×m；
步骤3、对历史工序编码矩阵依据kmeans聚类算法进行聚类：采用初始聚类中心计算方法、距离计算方法及更新聚类中心方法，得到不同K值下的聚类结果及聚类中心，根据不同K值下的聚类结果及聚类中心计算分类指标SP与DB，取分类指标SP/DB取最小时对应的K值、聚类中心及聚类结果，即为最终K值、最终聚类结果及最终聚类中心；
步骤4、订单拆分成工序，对待预测工序进行编码：将订单拆分成各个工序并判断各工序是否为外协加工，如果不是外协加工则认定为待预测工序，需进行工时预测，预测时首先对待预测的工序进行编码；
步骤5、将待预测工序划分至最终聚类组：
结合权重系数矩阵公式、加权欧式距离公式计算所有待预测工序到K个聚类中心的加权欧式距离，相互比较后将所有待预测工序划分至与自身距离最小的最终聚类组，并筛选出步骤1中采集的所有数据中属于该聚类组的历史零件图号、工艺信息、零件工序历史特征参数矩阵及零件工序历史工时数据；
步骤6、根据零件图号及工艺信息，将待预测工序与历史工序进行对比，若预测工序与历史工序的零件图号及工艺信息均相同，直接将历史理论标准工时作为预测理论工时并输出预测结果；若预测工序与历史工序不相同，则跳过该步骤直接进行下一步骤；
步骤7、得到历史工时z与PCA降至二维后的特征参数(x，y)的对应关系：用MatLab画出待预测工序所属聚类组中历史工时与特征参数的三维散点图，对三维散点图运用MatLab中用最小二乘法拟合生成相应的函数关系式z＝f(x，y)；
步骤8、预测工时计算：以新矩阵中s个待预测数据的第一维为x坐标，新矩阵中s个待预测数据的第二维为y坐标，代入关系式z＝f(x，y)得到待预测工序的理论预测工时公式并输出预测结果。
2.如权利要求1所述的基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，步骤3对历史工序编码矩阵依据kmeans聚类算法进行聚类，具体步骤包括：
3.1确定聚类组数K的取值范围是
Figure FDA0002631671290000021
随机选择一个K值；
3.2确定初始聚类中心，以满足初始聚类中心尽可能距离较远的要求；
3.3以K个聚类中心为基准，计算所有历史工序样本到K个聚类中心的距离，相互比较后将历史工序样本划分至距离最小的那一组，对所有历史工序样本划分完成后，得到历史工序样本的一次聚类结果；
在计算距离时，引入权重系数对历史工序样本各零件特征进行重要程度划分，权重系数的计算采用层次分析法构造判断矩阵，经归一化处理及一致性检验后得到其最终的权重系数矩阵S：
S＝{sj，j＝1，2，…，m} (4)
Figure FDA0002631671290000022
式中：sj-权重系数矩阵中第j列的数值；xij-历史工序编码矩阵Xn×m中第i个历史工序样本第j个特征的数值；crj-第r个聚类中心的第j个特征数值；dir-第i个样本到第r个聚类中心的加权欧式距离；
3.4根据3.3得到的聚类结果，更新原有聚类中心，得到K个新的聚类中心：
为降低历史工序样本中孤立点的负面影响，同时改善传统kmeans算法易陷入局部极小值的弊端，针对每个聚类组，分别更新原有聚类中心；
3.5对新聚类中心，按照步骤3.3和3.4进行迭代，直到聚类中心不再变化则完成聚类，得到该K值下聚类结果及聚类中心；注意在迭代的过程中权重系数矩阵S不发生变化；
3.6根据3.1步骤确定的K值范围，选择不同的K值，按照步骤3.2-3.5迭代得到不同K值下的聚类结果及聚类中心；
3.7计算不同聚类组数K值下的分类评价指标SP、DB，并取SP/DB取最小时对应的K值、聚类结果及聚类中心为最终K值、最终聚类结果及最终聚类中心。
3.如权利要求2所述的一种基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，确定初始聚类中心的公式为：
(1)根据历史工序编码矩阵Xn×m求得历史工序样本均值矩阵A1×m和历史工序样本极差矩阵B1×m：
Figure FDA0002631671290000031
B1×m＝{bj|bj＝xjmax-xjmin，j＝1，2，…，m} (2)
(2)根据历史工序样本均值矩阵A1×m和历史工序样本极差矩阵B1×m，计算聚类中心矩阵C：
Figure FDA0002631671290000032
式中：xij-历史工序编码矩阵Xn×m中第i个历史工序样本第j个特征的数值；
xjmax-历史工序编码矩阵Xn×m中第j个特征的最大值；
xjmin-历史工序编码矩阵Xn×m中第j特征的最小值；
aj-历史工序编码矩阵Xn×m中所有历史工序样本第j个特征的均值；
bj-历史工序编码矩阵Xn×m中第j个特征的极差值；
crj-第r个聚类中心的第j个特征数值；
RK×m-随机正太分布矩阵；
C-聚类中心矩阵。
4.如权利要求2所述的一种基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，步骤3.4更新原有聚类中心的具体步骤为：
3.4.1计算该聚类组下各历史工序样本两两之间的加权欧式距离diw,Ω为属于该聚类组的所有历史工序样本按Xn×m行递增的方式构成的矩阵，如若聚类组1中的数据有nΩ条，则聚类组1下的两两加权欧式距离应该有
Figure FDA0002631671290000041
个；
Figure FDA0002631671290000042
式中，diw-第i个历史工序样本到第w个历史工序样本的加权欧式距离；
xwj-历史工序编码矩阵Xn×m中第w个历史工序样本第j个特征的数值；
xij-历史工序编码矩阵Xn×m中第i个历史工序样本第j个特征的数值；
3.4.2求出3.4.1中所有加权欧式距离diw的平均值，作为平均距离dis_avrg；
Figure FDA0002631671290000043
3.4.3对该聚类组下的历史工序样本进行潜在孤立点的判断，得到密集点历史工序样本矩阵Φ及潜在孤立点历史工序样本矩阵Ψ：
判断规则为：对3.4.1中计算的所有加权欧式距离diw与平均距离dis_avrg比较，如若3.4.1中计算的diw＜dis_avrg，则认为历史工序样本i、w均为密集点样本；所有加权欧式距离均比较完成后，将得到的所有密集点样本去除重复项并取出密集点历史工序样本的数据构成矩阵Φ；由Φ∪Ψ＝Ω计算得到矩阵Ψ，矩阵Ψ即为所有潜在孤立点历史工序样本；求出矩阵Φ中各特征的平均值作为基础聚类中心cetr，同时统计矩阵Φ与矩阵Ψ中历史工序样本数量分别为nΦ，nΨ；
3.4.4对于潜在孤立点历史工序样本矩阵Ψ中的历史工序样本β，以一定的概率Pβ接受历史工序样本β的存在；同时采用随机方式对潜在孤立点历史工序样本矩阵Ψ中的样本的数据进行接受与舍弃，构成接受孤立点样本历史工序样本矩阵F：
Figure FDA0002631671290000051
式中：dβcetr-潜在孤立点历史工序样本矩阵Ψ中的历史工序样本β与基础聚类中心cetr的加权欧式距离；
MatLab中采用randn函数产生nΨ个标准正太分布的伪随机数，与3.4.4中的Pβ进行对应比较，若小于Pβ，则接受该潜在孤立点历史工序样本，若大于Pβ，则舍弃该潜在孤立点历史工序样本；统计出所有被接受的潜在孤立点历史工序样本，构成接受孤立点样本历史工序样本矩阵F，
Figure FDA0002631671290000052
接受孤立点样本历史工序样本矩阵F中历史工序样本数量为nF；
3.4.5由下式计算新聚类中心：
Figure FDA0002631671290000053
5.如权利要求2所述的基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，步骤3.7中分类评价指标SP、DB分别为：
根据不同K值下的聚类结果及聚类中心计算分类指标SP与DB：
Figure FDA0002631671290000061
Figure FDA0002631671290000062
式中：Chs-第h组与第s组聚类中心的加权欧式距离；
Wh-第h组内所有历史样本数据到该组聚类中心的平均加权欧式距离。
6.如权利要求1所述的基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，步骤6中的历史理论标准工时为：
Figure FDA0002631671290000063
式中：tp实—设备p加工的实做工时；
αp—设备p效率系数，由步骤1.1得到；
g—零件图号及工艺信息均相同的历史工序的个数；
p—与待预测工序零件图号及工艺信息均相同的第p个历史工序。
7.如权利要求1所述的基于kmeans聚类算法的回转类零件工序工时预测方法，其特征在于，步骤7用MatLab画出待预测工序所属聚类组中历史工时与特征参数的三维散点图，具体过程为：
统计待预测工序的工序历史特征参数矩阵，假设共有q个待预测工序，其中有s个待预测工序由步骤5被划分至同一个聚类组：聚类组1，将待预测工序的工序历史特征参数矩阵附加至聚类组1的工序历史特征参数矩阵的最后s行，构成附加矩阵，并对附加后的矩阵进行PCA降维，降至二维形成新矩阵，假设聚类组1中有l个历史样本数据，则新矩阵为(l+s)×2维；
以新矩阵中l个历史样本数据的第一维为x坐标，新矩阵中l个历史样本数据的第二维为y坐标，l个历史样本数据的历史工时为z坐标，利用Matlab画出三维散点图。
8.一种回转类零件工序工时预测系统，其特征在于，包括工序数据采集与处理模块、订单导入模块、工序编码模块、聚类模块、工时预测模块、数据存储模块；
所述工序数据采集与处理模块用于对历史工序的数据采集与处理；
所述订单导入模块用于将企业的零件订单导入该回转类零件工序工时预测系统；
所述工序编码模块用于对历史工序及待预测工序进行工序编码，同时输出历史工序编码结果及待预测工序编码结果；
所述聚类模块用于对历史工序编码进行kmeans聚类算法将历史工序编码聚类成组；
所述工时预测模块用于对待预测工序进行工时预测；
所述数据存储模块用以存储所述回转类零件工序工时预测系统的所有关键数据。
9.如权利要求8所述的回转类零件工序工时预测系统，其特征在于，所述聚类模块包括定期更新单元与聚类单元；
所述定期更新单元用于定期更新聚类结果信息，该定期更新单元将工序数据采集与处理模块中采集的信息定期传送至聚类模块中的聚类单元；
所述聚类单元用于对历史工序编码进行kmeans聚类算法将历史工序编码聚类成组。
10.如权利要求8所述的回转类零件工序工时预测系统，其特征在于，所述工时预测模块包括订单拆分单元、聚类组划分单元、历史匹配单元与工时计算单元；
所述订单拆分单元用于对订单导入模块导入的零件订单进行工序拆分，保留拆分的非外协加工工序为待预测工序；
所述聚类组划分单元用于将待预测工序划分至合适的历史聚类组中；
所述历史匹配单元用于根据聚类组划分单元输出信息中的零件图号及工艺信息对待预测工序与所属聚类组的历史工序进行匹配，并将匹配结果输出至工时计算单元；
所述工时计算单元用于对预测待加工工时计算并输出预测结果。
CN201710262953.9A 2017-04-20 2017-04-20 基于kmeans聚类算法的回转类零件工序工时预测方法及系统 Active CN108733003B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710262953.9A CN108733003B (zh) 2017-04-20 2017-04-20 基于kmeans聚类算法的回转类零件工序工时预测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710262953.9A CN108733003B (zh) 2017-04-20 2017-04-20 基于kmeans聚类算法的回转类零件工序工时预测方法及系统

Publications (2)

Publication Number Publication Date
CN108733003A CN108733003A (zh) 2018-11-02
CN108733003B true CN108733003B (zh) 2020-11-13

Family

ID=63933274

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710262953.9A Active CN108733003B (zh) 2017-04-20 2017-04-20 基于kmeans聚类算法的回转类零件工序工时预测方法及系统

Country Status (1)

Country Link
CN (1) CN108733003B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109902954B (zh) * 2019-02-27 2020-11-13 浙江工业大学 一种基于工业大数据的柔性作业车间动态调度方法
CN111489032B (zh) * 2020-04-09 2022-09-06 北京理工大学 一种预测航天产品装配工时的处理方法及装置
CN111833012B (zh) * 2020-06-19 2024-06-21 联想(北京)有限公司 一种工业数据处理方法和装置
CN112069317A (zh) * 2020-09-07 2020-12-11 北京理工大学 一种装配工时的获取方法及处理器
CN113298349A (zh) * 2021-04-25 2021-08-24 阿里巴巴（中国）有限公司 一种订单处理方法、设备及存储介质
CN113341883B (zh) * 2021-08-05 2021-11-09 山东豪泉软件技术有限公司 一种用于机床加工工时预测的方法及设备

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5438629A (en) * 1992-06-19 1995-08-01 United Parcel Service Of America, Inc. Method and apparatus for input classification using non-spherical neurons
CN102945517B (zh) * 2012-10-19 2016-02-10 利诚服装集团股份有限公司 一种基于聚类分析的服装标准工时的数据挖掘方法
CN103095534B (zh) * 2013-02-28 2016-01-20 南京理工大学常熟研究院有限公司 一种二维k均值熵的网络攻击效果评估方法
CN104376057A (zh) * 2014-11-06 2015-02-25 南京邮电大学 一种基于最大最小距离和K-means的自适应聚类方法
CN104574141A (zh) * 2014-12-23 2015-04-29 中国移动通信集团广东有限公司 一种业务影响度分析方法
CN105915451A (zh) * 2016-05-19 2016-08-31 东华大学 无线传感器网络的多sink部署与容错方法

Also Published As

Publication number Publication date
CN108733003A (zh) 2018-11-02

Similar Documents

Publication Publication Date Title
CN108733003B (zh) 基于kmeans聚类算法的回转类零件工序工时预测方法及系统
CN107220734B (zh) 基于决策树的数控车床车削过程能耗预测系统
CN107451747B (zh) 基于自适应非支配遗传算法的车间调度系统及其工作方法
CN106875094A (zh) 一种基于多色集合遗传算法的多目标车间调度方法
CN110689195A (zh) 一种电力日负荷预测方法
EP4080429A1 (en) Technology readiness level determination method and system based on science and technology big data
CN113454661A (zh) 产品不良成因分析的系统和方法、计算机可读介质
US20220277330A1 (en) Method and system for product processing price prediction based on multiple regression model
CN114239989A (zh) 一种物料需求计划的计算方法、系统、设备及存储介质
CN117055513B (zh) 基于人工智能的设备自动化控制方法及装置
CN117576823B (zh) 一种排队叫号系统终端
CN114548494A (zh) 一种可视化造价数据预测智能分析系统
CN113112121A (zh) 一种基于多目标非支配排序的车间布局调度优化方法
CN115841184A (zh) 一种基于工序聚类的工时预测方法
CN116502943A (zh) 一种熔模铸造产品质量追溯方法
CN115759455A (zh) 一种基于时序高斯混合密度网络的负荷概率密度预测方法
CN114926075A (zh) 基于工时预测的机械零件生产调度方法
CN109639283A (zh) 基于决策树的工件编码方法
CN113205274A (zh) 一种施工质量定量化排名方法
CN112070336A (zh) 一种基于层次分析法的制造业信息量化分析方法及装置
CN117556339B (zh) 一种网络违法行为风险危险等级评估方法
CN117371861B (zh) 基于数字化的家政服务质量智能分析方法及系统
CN113487080B (zh) 一种基于风速分类的风速动态场景生成方法、系统及终端
CN116823057B (zh) 一种基于工效分析的用户体验效果评价方法及系统
CN116070929B (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
GR01 Patent grant
GR01 Patent grant