CN110955865A - 一种基于粒子滤波的数据包络分析dea方法 - Google Patents

一种基于粒子滤波的数据包络分析dea方法 Download PDF

Info

Publication number
CN110955865A
CN110955865A CN201910997826.2A CN201910997826A CN110955865A CN 110955865 A CN110955865 A CN 110955865A CN 201910997826 A CN201910997826 A CN 201910997826A CN 110955865 A CN110955865 A CN 110955865A
Authority
CN
China
Prior art keywords
particle
particles
dea
value
fitness
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.)
Granted
Application number
CN201910997826.2A
Other languages
English (en)
Other versions
CN110955865B (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.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
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 Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN201910997826.2A priority Critical patent/CN110955865B/zh
Publication of CN110955865A publication Critical patent/CN110955865A/zh
Application granted granted Critical
Publication of CN110955865B publication Critical patent/CN110955865B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]

Landscapes

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

Abstract

一种基于粒子滤波的数据包络分析DEA方法，首先，将DEA非线性约束优化问题的求解过程转换为动态系统的状态估计过程，通过大量粒子的重要性采样、权值更新、重采样、状态估计等过程，以迭代搜索的方式，逐步寻找出对于决策单元本身最优的投入产出方案。本发明根据复杂系统投入的生产要素与产出的产品之间的相关性，建立非线性约束优化问题模型。采用粒子滤波的方法，通过迭代搜索的方式寻找出最优决策单元的线性组合系数以及相对效率。用于评价各部门间的相对有效性，从而找出对于决策单元本身最优的投入产出方案。

Description

一种基于粒子滤波的数据包络分析DEA方法
技术领域
本发明涉及运筹学、管理科学与数理经济学交叉研究领域，具体涉及一种基于粒子滤波的数据包络分析DEA方法。
技术背景
数据包络分析(Data Envelopment Analysis)简称DEA，是运筹学、管理科学的一个新的研究领域。该方法是由著名的运筹学家A.Charmes、W.W.Cooper及 E.Rhodes首先提出来的，它主要用于评价相同部门间的相对有效性。该方法能充分考虑对于决策单元本身最优的投入产出方案，因而能够更理想地反映评价对象自身的信息和特点；同时对于评价复杂系统的多投入多产出分析具有独到之处。目前，DEA的优点吸引众多的应用者。它的应用范围从军用飞机、基地维修保养、银行、交通，证券投资、公共事业到企业规模、项目评价等。
DEA方法的一般模型为C2R模型，该模型是指对多个决策单元(DMU)通过“投入一定数量的生产要素，并产出一定数量的产品”的经济系统来判断各个单元的投入产出相对于其他决策单元的相对合理性和相对有效性。
C2R模型可以写成如下形式：
Figure BDA0002238558770000011
其中，m为整数，表示系统投入生产要素的数量；n为整数，表示DMU决策单元数；s为整数，表示系统产出产品的数量；xij,(i＝1,2,…,m；j＝1,2,…,n) 表示DMU-j决策单元对第i种输入的投入量；yrj,(r＝1,2,…,s；j＝1,2,…,n)表示DMU-j决策单元对第r种输出的产出量；
Figure BDA0002238558770000012
表示第i种输入的参考投入量；
Figure BDA0002238558770000013
表示第r种输出的参考产出量。E表示投入产出比率，wj(j＝1,2,…,n)表示DMU-j决策单元的线性组合系数。
C2R模型是用于评价DMU的规模和技术有效性的。当需要研究各部门之间技术有效性时，这时需要用到它的改进模型“C2G S2”模型。这个模型只不过是在C2R模型的基础上增加了约束条件
Figure BDA0002238558770000021
这样就构成了“C2G S2”模型：
Figure BDA0002238558770000022
目前，主要运用线性规划软件(QSB、Lindo)、统计软件SAS，但是，这些软件多为专业人士使用，并不为更多非专业人员熟知。另外比较常用的一种方法是Excel数据处理软件，该方法需要用户对Excel数据处理软件、求解步骤有比较入的了解，并且还要使用户参与到求解过程中，操作过程比较烦琐。因此，如何能够高效精确地求解数据包络分析DEA问题，仍是一个值得研究的课题。
发明内容
针对多投入多产出的复杂系统，为了克服目前多投入多产出系统中DEA问题求解过程复杂、精度较低的问题，本发明提出一种基于粒子滤波的数据包络分析方法。根据复杂系统投入的生产要素与产出的产品之间的相关性，建立非线性约束优化问题模型。采用粒子滤波的方法，通过迭代搜索的方式寻找出最优决策单元的线性组合系数以及相对效率。用于评价各部门间的相对有效性，从而找出对于决策单元本身最优的投入产出方案。
为了解决上述技术问题，本发明提供如下的技术方案：
一种基于粒子滤波的数据包络分析DEA方法，包括以下步骤：
步骤一，将“C2G S2”数学模型转换为非线性约束优化问题，定义投入生产要素的约束函数和产出产品的约束函数分别为公式(1)和公式(2)：
Figure BDA0002238558770000023
Figure BDA0002238558770000031
其中，w＝{w1,w2,…,wn}∈R1×n表示系统DMU决策单元的线性组合系数向量，且有0≤wj≤1,(j＝1,2,…,n)，m为整数，表示系统投入生产要素的数量；n为整数，表示DMU决策单元数；s为整数，表示系统产出产品的数量；
xij,(i＝1,2,…,m；j＝1,2,…,n)表示DMU-j决策单元对第i种输入的投入量；
yrj(r＝1,2,…,s；j＝1,2,…,n)表示DMU-j决策单元对第r种输出的产出量；
xij0(i＝1,2,…,m)表示第i种输入的参考投入量；yrj0(r＝1,2,…,s)表示第r种输出的参考产出量；E表示投入产出比率，wj(j＝1,2,…,n)表示DMU-j决策单元的线性组合系数；那么，“C2G S2”数学模型转换为以下非线性约束优化问题：
Figure BDA0002238558770000032
这样，DEA分析问题就转换为一个形式简单的非线性约束优化问题，优化适应度函数为fitness(E)；
步骤二，用动态时变系统来描述数据包络分析DEA非线性约束优化问题的求解过程，并进行建模；如果采用迭代优化的方式来求解如公式(3)所示的DEA 非线性约束优化问题，那么该求解过程即可视为是一个动态时变系统：用离散的时间量来表示迭代搜索次数，用系统状态值来表示每次迭代的局部最优解，那么该动态时变系统的运动模型即描述了数据包络分析DEA非线性约束优化问题的求解过程，该动态时变系统的观测模型即描述了DEA非线性约束优化问题中局部最优解的更新过程，系统的运动模型和观测模型分别用以下公式(4)和公式(5) 来描述：
vk＝fk(vk-1,uk) (4)
zk＝fitness(vk) (5)
其中，vk为系统在k,(k＝1,2,…,K)时刻的状态量，K为系统总时长；uk为系统在k,(k＝1,2,…,K)时刻的过程噪声，fk(·)是k时刻状态量vk与k-1时刻状态量vk-1之间的关系函数，是一种非确定性函数，函数的形式取决于数据包络分析DEA问题的优化过程与收敛速度；
步骤三，系统状态初始化，设定粒子群中粒子总数为P，那么第 p,(p＝1,2,…,P)个粒子在k,(k＝1,2,…,K)时刻表示为
Figure BDA0002238558770000041
将各粒子的初始值均初始化为
Figure BDA0002238558770000042
将系统的状态值初始化为
Figure BDA0002238558770000043
系统的观测值初始化为 z1＝fitness(v1)，系统的最优解初始化为vbest＝v1，最优适应度值初始化为 zbest＝z1；
步骤四，重要性采样：对于第p,(p＝1,2,…,P)个粒子，根据概率密度分布
Figure BDA0002238558770000044
采集新粒子
Figure BDA0002238558770000045
其中k,(k＝1,2,…,K)表示系统离散的时间量；这里的概率密度分布
Figure BDA0002238558770000046
是由公式六所示的系统状态方程中的非确定性函数fk(vk-1,uk)决定的；根据系统假设，函数的寻优过程应该是一个搜索范围逐渐缩小的过程，采用均匀分布的方式进行搜索，假设概率密度分布
Figure BDA0002238558770000047
用均匀分布
Figure BDA0002238558770000048
来近似，其中ck是一个待定的参量，该参量的取值应当满足以下原则：随着时间的推移，迭代次数k递增，而参量ck的值递减，参量ck在每次迭代过程中减小的幅度将直接影响系统寻找最优解的速度和精度；
步骤五，更新全局最优解。根据系统观测方程中的适应度函数fitness(·)，对所有重要性采样后的粒子
Figure BDA0002238558770000049
进行评价，计算其适应度值(或观测值)
Figure BDA00022385587700000410
对于任意第p,(p＝1,2,…,P)粒子
Figure BDA00022385587700000411
都进行判断，当
Figure BDA00022385587700000412
时，更新全局最优解
Figure BDA00022385587700000413
和全局最优适应度值
Figure BDA00022385587700000414
而当
Figure BDA00022385587700000415
系统全局最优解vbest和最优适应度值zbest保持不变；
步骤六，更新粒子权值，首先，判断粒子是否有效，即对于不满足公式(3) 所示的非线性约束条件的粒子，其权值直接置零，即
Figure BDA00022385587700000416
然后，计算有效粒子的权值，对于满足公式(3)所示的非线性约束条件的粒子，若粒子的适应度值(或观测值)
Figure BDA0002238558770000051
比系统的当前状态值xk大，即
Figure BDA0002238558770000052
同样将其权值置零，即
Figure BDA0002238558770000053
若粒子的适应度值(或观测值)
Figure BDA0002238558770000054
小于或等于系统的当前状态值xk，即
Figure BDA0002238558770000055
则衡量粒子的适应度值(或观测值)与系统的当前状态值之间的欧式距离，给予欧式距离小的粒子一个较小的权值，而给予欧式距离大的粒子一个较大的权值；
当粒子数足够多的时候，粒子的适应度值(或观测值)服从正态分布
Figure BDA0002238558770000056
其中S2为样本方差，据此，根据以下公式(6)和公式(7) 计算各粒子的权值；
Figure BDA0002238558770000057
Figure BDA0002238558770000058
在得到所有P个粒子的权值
Figure BDA0002238558770000059
之后，对其进行归一化处理，即根据以下公式(8)归一化权值；
Figure BDA00022385587700000510
步骤七，粒子重采样：首先，设定一个粒子匮乏判断阈值Nth，取Nth＝2P/3，假设当有效粒子数占总粒子数2/3以上时，判定粒子数量充足，有效粒子数的计算方式如下：
Figure BDA00022385587700000511
如果Neff＜Nth，则启动粒子重采样过程，采用轮盘赌的方式重新抽样；如果 Neff≥Nth，则不需要进行粒子重采样过程；
步骤八，系统状态更新，当前k,(k＝1,2,…,K)时刻系统的状态量表示为各粒子的加权平均，即：
Figure BDA0002238558770000061
步骤九，判断迭代终止条件，输出最优解；如未满足迭代次数，即当前迭代次数k＜K，那么返回步骤四，粒子重要性采样；当满足迭代次数，即当前迭代次数k＝K，那么粒子滤波过程结束，输出系统最优解vbest及其对应的适应度值 zbest；至此，该数据包络分析DEA问题得以求解，最优的投入产出比率 Ebest＝vbest[1]。
本发明的有益效果表现在：将DEA问题“C2G S2”模型的优化过程描述为一个动态时变系统的状态估计过程。在有效定义区间内，随着离散时间量(或迭代次数)的逐步递增，系统的观测值(或适应度值)逐步朝着最优的投入产出比率靠近。通过大量粒子的重要性采样、权值更新、重采样、状态估计等过程，寻找出系统最优的投入产出比率。通过理论分析及实验结果得到，本发明方法在求解数据包络分析DEA问题上具有良好的稳定性和较高的求解精度，是一种行之有效的方法。
具体实施方式
下面对本发明做进一步说明。
一种基于粒子滤波的数据包络分析DEA方法，包括以下步骤：
步骤一，将“C2G S2”数学模型转换为非线性约束优化问题，定义投入生产要素的约束函数和产出产品的约束函数分别为公式(1)和公式(2)：
Figure BDA0002238558770000062
Figure BDA0002238558770000063
其中，w＝{w1,w2,…,wn}∈R1×n表示系统DMU决策单元的线性组合系数向量，且有0≤wj≤1,(j＝1,2,…,n)，其他参数定义如下：m为整数，表示系统投入生产要素的数量；n为整数，表示DMU决策单元数；s为整数，表示系统产出产品的数量；xij,(i＝1,2,…,m；j＝1,2,…,n)表示DMU-j决策单元对第i种输入的投入量；yrj(r＝1,2,…,s；j＝1,2,…,n)表示DMU-j决策单元对第r种输出的产出量；
Figure BDA0002238558770000064
表示第i种输入的参考投入量；
Figure BDA0002238558770000065
表示第r种输出的参考产出量；E表示投入产出比率，wj(j＝1,2,…,n)表示DMU-j 决策单元的线性组合系数；那么，“C2G S2”数学模型转换为以下非线性约束优化问题：
Figure BDA0002238558770000071
这样，DEA分析问题就转换为一个形式简单的非线性约束优化问题，优化适应度函数为fitness(E)；
步骤二，用动态时变系统来描述数据包络分析DEA非线性约束优化问题的求解过程，并进行建模；如果采用迭代优化的方式来求解如公式(3)所示的DEA 非线性约束优化问题，那么该求解过程即可视为是一个动态时变系统：用离散的时间量来表示迭代搜索次数，用系统状态值来表示每次迭代的局部最优解，那么该动态时变系统的运动模型即描述了数据包络分析DEA非线性约束优化问题的求解过程，该动态时变系统的观测模型即描述了DEA非线性约束优化问题中局部最优解的更新过程，系统的运动模型和观测模型分别用以下公式(4)和公式(5) 来描述。
vk＝fk(vk-1,uk) (4)
zk＝fitness(vk) (5)
其中，vk为系统在k,(k＝1,2,…,K)时刻的状态量，K为系统总时长；uk为系统在k,(k＝1,2,…,K)时刻的过程噪声，fk(·)是k时刻状态量vk与k-1时刻状态量vk-1之间的关系函数，是一种非确定性函数，函数的形式取决于数据包络分析DEA问题的优化过程与收敛速度；
步骤三，系统状态初始化，设定粒子群中粒子总数为P，那么第 p,(p＝1,2,…,P)个粒子在k,(k＝1,2,…,K)时刻表示为
Figure BDA0002238558770000072
将各粒子的初始值均初始化为
Figure BDA0002238558770000081
将系统的状态值初始化为
Figure BDA0002238558770000082
系统的观测值初始化为 z1＝fitness(v1)，系统的最优解初始化为vbest＝v1，最优适应度值初始化为 zbest＝z1；
步骤四，重要性采样：对于第p,(p＝1,2,…,P)个粒子，根据概率密度分布
Figure BDA0002238558770000083
采集新粒子
Figure BDA0002238558770000084
其中k,(k＝1,2,…,K)表示系统离散的时间量，这里的概率密度分布
Figure BDA0002238558770000085
是由公式六所示的系统状态方程中的非确定性函数fk(vk-1,uk)决定的，根据系统假设，函数的寻优过程应该是一个搜索范围逐渐缩小的过程，采用均匀分布的方式进行搜索，假设概率密度分布
Figure BDA0002238558770000086
用均匀分布
Figure BDA0002238558770000087
来近似，其中ck是一个待定的参量，该参量的取值应当满足以下原则：随着时间的推移，迭代次数k递增，而参量ck的值递减，参量ck在每次迭代过程中减小的幅度将直接影响系统寻找最优解的速度和精度；
步骤五，更新全局最优解，根据系统观测方程中的适应度函数fitness(·)，对所有重要性采样后的粒子
Figure BDA0002238558770000088
进行评价，计算其适应度值(或观测值)
Figure BDA0002238558770000089
对于任意第p,(p＝1,2,…,P)粒子
Figure BDA00022385587700000810
都进行判断，当
Figure BDA00022385587700000811
时，更新全局最优解
Figure BDA00022385587700000812
和全局最优适应度值
Figure BDA00022385587700000813
而当
Figure BDA00022385587700000814
系统全局最优解vbest和最优适应度值zbest保持不变；
步骤六，更新粒子权值，首先，判断粒子是否有效，即对于不满足公式(3) 所示的非线性约束条件的粒子，其权值直接置零，即
Figure BDA00022385587700000815
然后，计算有效粒子的权值，对于满足公式(3)所示的非线性约束条件的粒子，若粒子的适应度值(或观测值)
Figure BDA00022385587700000816
比系统的当前状态值xk大，即
Figure BDA00022385587700000817
同样将其权值置零，即
Figure BDA00022385587700000818
若粒子的适应度值(或观测值)
Figure BDA00022385587700000819
小于或等于系统的当前状态值xk，即
Figure BDA00022385587700000820
则衡量粒子的适应度值(或观测值)与系统的当前状态值之间的欧式距离，给予欧式距离小的粒子一个较小的权值，而给予欧式距离大的粒子一个较大的权值，具体实施方式如下：当粒子数足够多的时候，粒子的适应度值(或观测值)服从正态分布
Figure BDA0002238558770000091
其中S2为样本方差，据此，根据以下公式(6)和公式(7)计算各粒子的权值。
Figure BDA0002238558770000092
Figure BDA0002238558770000093
在得到所有P个粒子的权值
Figure BDA0002238558770000094
之后，对其进行归一化处理，即根据以下公式(8)归一化权值。
Figure BDA0002238558770000095
步骤七，粒子重采样，为了减小粒子匮乏对系统收敛性的影响，需要启动粒子重采样过程，首先，设定一个粒子匮乏判断阈值Nth，取Nth＝2P/3，假设当有效粒子数占总粒子数2/3以上时，判定粒子数量充足，有效粒子数的计算方式如下：
Figure BDA0002238558770000096
如果Neff＜Nth，则启动粒子重采样过程，可以采用轮盘赌的方式重新抽样；如果Neff≥Nth，则不需要进行粒子重采样过程；
步骤八，系统状态更新。当前k,(k＝1,2,…,K)时刻系统的状态量可以表示为各粒子的加权平均，即：
Figure BDA0002238558770000097
步骤九，判断迭代终止条件，输出最优解，如未满足迭代次数，即当前迭代次数k＜K，那么返回步骤四，粒子重要性采样；当满足迭代次数，即当前迭代次数k＝K，那么粒子滤波过程结束，输出系统最优解vbest及其对应的适应度值 zbest，至此，该数据包络分析DEA问题得以求解，最优的投入产出比率 Ebest＝vbest[1]。
实验对比：为了验证本发明方法求解数据包络分析问题的效果，这里采用文献中所使用的快餐连锁店问题来进行测试。该问题是对连锁店经营业绩的评价问题，有关10家快餐连锁店效率的数据如表1所示。
Figure BDA0002238558770000101
表1
按照数据包络法的“C2G S2”模型就可以得到各个连锁店的效率评价模型。这里以连锁店4为例，可以得到下面的评价模型：
Figure 1
其他连锁店也仿照此模型进行处理，最后通过求解该线性约束优化问题模型就可以求得各个连锁店的相对效率E。
为了考察本发明方法求解数据包络分析问题的性能，采用Matlab 7.1实验平台作为仿真环境，对10家连锁店模型分别进行30次测试，并与其他文献中采用的线性规划方法得到的结果相比较，仿真实验结果如表2、表3所示。
连锁店序号 最好解 最差解 样本均值 样本方差
1 0.94405718 0.944057688 0.944057383 5.6E-14
2 1 1 1 0
3 1 1 1 0
4 0.885030477 0.885037857 0.885033001 2.82E-12
5 0.759749435 0.759754596 0.75975011 1.997E-12
6 1 1 1 0
7 0.966608323 0.966613484 0.966609011 3.184E-12
8 1 1 1 0
9 1 1 1 0
10 0.945087952 0.945098411 0.945089147 1.0069E-11
表2
Figure BDA0002238558770000112
Figure BDA0002238558770000121
表3
表2给出了本发明方法对连锁店DEA分析的30次仿真实验所得结果的统计分析，从表中数据可以看出，这些样本的方差都很小，接近于0，即这30次仿真实验的结果差异性较小，可见本发明方法的稳定性较高。表3将本发明方法 30次仿真的结果均值与文献中采用的线性规划方法得到的结果进行比较，并计算了绝对误差，从表中的数据可以看出，利用本发明方法得到的结果误差较小，可见本发明方法的求解精度较高。因此，本发明方法在求解该数据包络分析DEA 问题上具有良好的稳定性和较高的求解精度，是一种行之有效的方法。

Claims (5)

1.一种基于粒子滤波的数据包络分析DEA方法，其特征在于，所述方法包括以下步骤：
步骤一，将“C2G S2”数学模型转换为非线性约束优化问题，定义投入生产要素的约束函数和产出产品的约束函数分别为公式(1)和公式(2)：
Figure FDA0002238558760000011
Figure FDA0002238558760000012
其中，w＝{w1,w2,…,wn}∈R1×n表示系统DMU决策单元的线性组合系数向量，且有0≤wj≤1,(j＝1,2,…,n)，m为整数，表示系统投入生产要素的数量；n为整数，表示DMU决策单元数；s为整数，表示系统产出产品的数量；
xij,(i＝1,2,…,m；j＝1,2,…,n)表示DMU-j决策单元对第i种输入的投入量；
yrj(r＝1,2,…,s；j＝1,2,…,n)表示DMU-j决策单元对第r种输出的产出量；
Figure FDA0002238558760000013
表示第i种输入的参考投入量；
Figure FDA0002238558760000014
表示第r种输出的参考产出量；E表示投入产出比率，wj(j＝1,2,…,n)表示DMU-j决策单元的线性组合系数；那么，“C2G S2”数学模型转换为以下非线性约束优化问题：
Figure FDA0002238558760000015
这样，DEA分析问题就转换为一个形式简单的非线性约束优化问题，优化适应度函数为fitness(E)；
步骤二，用动态时变系统来描述数据包络分析DEA非线性约束优化问题的求解过程，并进行建模，如果采用迭代优化的方式来求解如公式(3)所示的DEA非线性约束优化问题，那么该求解过程即可视为是一个动态时变系统：用离散的时间量来表示迭代搜索次数，用系统状态值来表示每次迭代的局部最优解，那么该动态时变系统的运动模型即描述了数据包络分析DEA非线性约束优化问题的求解过程，该动态时变系统的观测模型即描述了DEA非线性约束优化问题中局部最优解的更新过程，系统的运动模型和观测模型分别用以下公式(4)和公式(5)来描述：
vk＝fk(vk-1,uk) (4)
zk＝fitness(vk) (5)
其中，vk为系统在k,(k＝1,2,…,K)时刻的状态量，K为系统总时长；uk为系统在k,(k＝1,2,…,K)时刻的过程噪声，fk(·)是k时刻状态量vk与k-1时刻状态量vk-1之间的关系函数，是一种非确定性函数，函数的形式取决于数据包络分析DEA问题的优化过程与收敛速度；
步骤三，系统状态初始化，设定粒子群中粒子总数为P，那么第p,(p＝1,2,…,P)个粒子在k,(k＝1,2,…,K)时刻表示为
Figure FDA0002238558760000021
将各粒子的初始值
Figure FDA0002238558760000022
初始化，将系统的状态值v1初始化，系统的观测值初始化为z1＝fitness(v1)，系统的最优解初始化为vbest＝v1，最优适应度值初始化为zbest＝z1；
步骤四，重要性采样：对于第p,(p＝1,2,…,P)个粒子，根据概率密度分布
Figure FDA0002238558760000023
采集新粒子
Figure FDA0002238558760000024
其中k,(k＝1,2,…,K)表示系统离散的时间量，这里的概率密度分布
Figure FDA0002238558760000025
是由公式(4)所示的系统状态方程中的非确定性函数fk(vk-1,uk)决定的，根据系统假设，函数的寻优过程应该是一个搜索范围逐渐缩小的过程，采用均匀分布的方式进行搜索；
步骤五，更新全局最优解，根据系统观测方程中的适应度函数fitness(·)，对所有重要性采样后的粒子
Figure FDA0002238558760000026
进行评价，计算其适应度值
Figure FDA0002238558760000027
对于任意第p,(p＝1,2,…,P)粒子
Figure FDA0002238558760000028
都进行判断，当
Figure FDA0002238558760000029
时，更新全局最优解
Figure FDA00022385587600000210
和全局最优适应度值
Figure FDA00022385587600000211
而当
Figure FDA00022385587600000212
系统全局最优解vbest和最优适应度值zbest保持不变；
步骤六，更新粒子权值，首先，判断粒子是否有效，即对于不满足公式(3)所示的非线性约束条件的粒子，其权值直接置零，即
Figure FDA00022385587600000213
然后，计算有效粒子的权值，对于满足公式(3)所示的非线性约束条件的粒子，若粒子的适应度值
Figure FDA00022385587600000214
比系统的当前状态值xk大，即
Figure FDA00022385587600000215
同样将其权值置零，即
Figure FDA00022385587600000216
若粒子的适应度值
Figure FDA00022385587600000217
小于或等于系统的当前状态值xk，即
Figure FDA00022385587600000218
则衡量粒子的适应度值与系统的当前状态值之间的欧式距离，给予欧式距离小的粒子一个较小的权值，而给予欧式距离大的粒子一个较大的权值；
步骤七，粒子重采样；
步骤八，系统状态更新，当前k,(k＝1,2,…,K)时刻系统的状态量表示为各粒子的加权平均，即：
Figure FDA00022385587600000219
步骤九，判断迭代终止条件，输出最优解，如未满足迭代次数，即当前迭代次数k＜K，那么返回步骤四，粒子重要性采样；当满足迭代次数，即当前迭代次数k＝K，那么粒子滤波过程结束，输出系统最优解vbest及其对应的适应度值zbest，至此，该数据包络分析DEA问题得以求解，最优的投入产出比率Ebest＝vbest[1]。
2.根据权利要求1所述的一种基于粒子滤波的数据包络分析DEA方法，其特征在于，所述步骤三中，所述的系统状态初始化中，各粒子的初始值
Figure FDA0002238558760000031
应在区间范围
Figure FDA0002238558760000032
内选取，设置为
Figure FDA0002238558760000033
而系统的状态值初始值应设置为
Figure FDA0002238558760000034
3.根据权利要求1或2所述的一种基于粒子滤波的数据包络分析DEA方法，其特征在于，所述步骤四中，所述的概率密度函数
Figure FDA0002238558760000035
应当采用均匀分布
Figure FDA0002238558760000036
来近似，其中ck是一个待定的参量，该参量的取值应当满足以下原则：随着时间的推移，迭代次数k递增，而参量ck的值递减，参量ck在每次迭代过程中减小的幅度将直接影响系统寻找最优解的速度和精度。
4.根据权利要求1或2所述的一种基于粒子滤波的数据包络分析DEA方法，其特征在于，所述步骤六中，所述的更新粒子权值的具体实施方式如下：当粒子数足够多的时候，粒子的适应度值服从正态分布
Figure FDA0002238558760000037
其中S2为样本方差，据此，根据以下公式(6)和公式(7)计算各粒子的权值：
Figure FDA0002238558760000038
Figure FDA0002238558760000039
在得到所有P个粒子的权值
Figure FDA00022385587600000310
之后，对其进行归一化处理，即根据以下公式(8)归一化权值：
Figure FDA00022385587600000311
5.根据权利要求4所述的一种基于粒子滤波的数据包络分析DEA方法，其特征在于，所述步骤七中，所述的粒子重采样的具体实施方式如下：首先，设定一个粒子匮乏判断阈值Nth，取Nth＝2P/3，假设当有效粒子数占总粒子数2/3以上时，判定粒子数量充足，有效粒子数的计算方式如下：
Figure FDA0002238558760000041
如果Neff＜Nth，则启动粒子重采样过程，可以采用轮盘赌的方式重新抽样；如果Neff≥Nth，则不需要进行粒子重采样过程。
CN201910997826.2A 2019-10-18 2019-10-18 一种基于粒子滤波的数据包络分析dea方法 Active CN110955865B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910997826.2A CN110955865B (zh) 2019-10-18 2019-10-18 一种基于粒子滤波的数据包络分析dea方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910997826.2A CN110955865B (zh) 2019-10-18 2019-10-18 一种基于粒子滤波的数据包络分析dea方法

Publications (2)

Publication Number Publication Date
CN110955865A true CN110955865A (zh) 2020-04-03
CN110955865B CN110955865B (zh) 2023-12-29

Family

ID=69975736

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910997826.2A Active CN110955865B (zh) 2019-10-18 2019-10-18 一种基于粒子滤波的数据包络分析dea方法

Country Status (1)

Country Link
CN (1) CN110955865B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6334099B1 (en) * 1999-05-25 2001-12-25 Digital Gene Technologies, Inc. Methods for normalization of experimental data
CN106338651A (zh) * 2016-08-31 2017-01-18 长沙理工大学 应用于电力系统低频振荡模式识别的粒子滤波分析方法
WO2018072351A1 (zh) * 2016-10-20 2018-04-26 北京工业大学 一种基于粒子群优化算法对支持向量机的优化方法
CN108564245A (zh) * 2018-02-27 2018-09-21 国网冀北电力有限公司电力科学研究院 一种基于dea分析的配网投资效率评价方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6334099B1 (en) * 1999-05-25 2001-12-25 Digital Gene Technologies, Inc. Methods for normalization of experimental data
CN106338651A (zh) * 2016-08-31 2017-01-18 长沙理工大学 应用于电力系统低频振荡模式识别的粒子滤波分析方法
WO2018072351A1 (zh) * 2016-10-20 2018-04-26 北京工业大学 一种基于粒子群优化算法对支持向量机的优化方法
CN108564245A (zh) * 2018-02-27 2018-09-21 国网冀北电力有限公司电力科学研究院 一种基于dea分析的配网投资效率评价方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
吴新杰;黄国兴;: "利用粒子滤波求解旅行商问题", 计算机应用, no. 08 *

Also Published As

Publication number Publication date
CN110955865B (zh) 2023-12-29

Similar Documents

Publication Publication Date Title
CN112382352B (zh) 基于机器学习的金属有机骨架材料结构特征快速评估方法
CN108228716B (zh) 基于加权极限学习机的SMOTE_Bagging集成污水处理故障诊断方法
CN108595916B (zh) 基于生成对抗网络的基因表达全谱推断方法
Gong et al. Multi-objective parameter optimization of common land model using adaptive surrogate modeling
Kuptametee et al. A review of resampling techniques in particle filtering framework
CN106503867A (zh) 一种遗传算法最小二乘风电功率预测方法
CN110705029B (zh) 一种基于迁移学习的振荡扑翼能量采集系统流场预测方法
CN109558893B (zh) 基于重采样池的快速集成污水处理故障诊断方法
CN108846261B (zh) 基于可视图算法的基因表达时序数据分类方法
CN113486591B (zh) 一种卷积神经网络结果的重力多参量数据密度加权反演方法
CN111242867B (zh) 基于截断泰勒级数近似的图信号分布式在线重构方法
CN106296434B (zh) 一种基于pso-lssvm算法的粮食产量预测方法
CN112613658A (zh) 逐日降雨量预测方法、装置、电子设备及存储介质
CN116993548A (zh) 基于增量学习的LightGBM-SVM的教育培训机构信用评估方法及系统
CN115237878A (zh) 基于增材制造的工艺数据库构建方法及介质
CN109919374A (zh) 基于apso-bp神经网络的股票价格预测方法
CN107577896B (zh) 基于混合Copula理论的风电场多机聚合等值方法
Ao et al. Entropy estimation via normalizing flow
CN107766887A (zh) 一种局部加权的不完整数据混杂聚类方法
Niemi et al. Adaptive mixture modeling Metropolis methods for Bayesian analysis of nonlinear state-space models
CN116865255A (zh) 基于改进熵权法和seceemd的短期风电功率预测方法
CN110955865A (zh) 一种基于粒子滤波的数据包络分析dea方法
Li et al. Time series clustering based on relationship network and community detection
CN114819107B (zh) 基于深度学习的混合数据同化方法
CN116680969A (zh) 一种pso-bp算法的充填体评估参数预测方法及装置

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