超效率SBM模型Python实战用scipy.optimize处理含非期望产出的政府数据效率排名在公共管理领域政府开放数据的效率评估一直是学术界和实务界关注的焦点。传统的DEA模型虽然被广泛应用但在处理包含非期望产出如污染排放、能源消耗的复杂场景时往往显得力不从心。超效率SBMSlack-Based Measure模型作为一种非径向、非角度的效率评估方法不仅能够有效处理期望产出与非期望产出并存的情况还能对效率值为1的决策单元进行进一步区分排序为政策制定者提供更精细的决策依据。本文将带您深入超效率SBM模型的数学原理和Python实现细节使用scipy.optimize从零构建完整的求解流程。不同于简单的库调用我们将重点关注如何将数学模型精确转化为线性规划问题处理多类型产出数据时的reshape和拼接技巧约束条件矩阵(Aeq)的构建逻辑与边界条件(bounds)设置实际政府数据评估中的常见陷阱与解决方案1. 超效率SBM模型的核心原理1.1 从传统DEA到SBM的演进传统DEA模型如CCR、BCC存在两个主要局限径向性限制要求所有投入或产出按相同比例增减无法区分有效单元当多个DMU效率值为1时无法进一步排序超效率SBM模型通过引入松弛变量解决了这些问题# 关键改进对比 traditional_dea { 径向性: 是, 效率值范围: [0,1], 非期望产出处理: 困难 } super_sbm { 径向性: 否, 效率值范围: [0,∞), 非期望产出处理: 原生支持 }1.2 数学建模要点考虑有n个决策单元(DMU)每个DMU有m种投入x ∈ R^ms1种期望产出y^g ∈ R^s1s2种非期望产出y^b ∈ R^s2超效率SBM模型的目标函数为min ρ (1 - (1/m)∑(s_i^-/x_ik)) / (1 (1/(s1s2))(∑(s_r^g/y_rk^g) ∑(s_r^b/y_rk^b)))其中s_i^- 为投入松弛量s_r^g 为期望产出松弛量s_r^b 为非期望产出松弛量注意非期望产出的处理需要特别注意方向性其松弛量应被视为需要减少的部分2. 数据准备与预处理2.1 政府开放数据的典型结构以16个城市的政府数据为例常见指标包括指标类型具体指标示例投入指标财政投入、人力投入、IT基础设施期望产出数据开放数量、API调用次数、公众满意度非期望产出数据泄露事件、响应延迟、能源消耗import pandas as pd import numpy as np # 模拟数据加载 data pd.read_excel(government_data.xlsx, index_col0) raw_values data.values # 数据拆分 x raw_values[:, 0:2].T # 2个投入指标 y_g raw_values[:, 2:4].T # 2个期望产出 y_b raw_values[:, 4:5].T # 1个非期望产出2.2 非期望产出的特殊处理非期望产出需要反向处理常见方法包括取倒数法1/y_b线性变换max(y_b) 1 - y_b方向性距离函数直接纳入模型约束# 方法1取倒数处理 y_b_processed 1 / y_b # 方法2线性变换 y_b_max np.max(y_b) y_b_processed y_b_max 1 - y_b提示实际应用中建议测试不同处理方法对结果的影响3. scipy.optimize实现详解3.1 线性规划问题构建将超效率SBM转化为标准LP问题min c^T * v s.t. Aeq * v beq v 0其中决策变量v包含n个λ权重变量m个投入松弛量s1个期望产出松弛量s2个非期望产出松弛量1个ρ效率变量from scipy import optimize m, n x.shape s1 y_g.shape[0] s2 y_b.shape[0] theta [] for i in range(n): # 对每个DMU单独求解 # 目标函数系数 c np.concatenate([ np.zeros(n), # λ -1/(m * x[:, i]), # s_i^- (投入松弛) np.zeros(s1 s2), # s_r^g, s_r^b (产出松弛) np.array([1]) # ρ ]) # 等式约束矩阵 Aeq1 np.hstack([ x, # λ*x np.eye(m), # s_i^- np.zeros((m, s1 s2)), # 产出松弛占位 -x[:, i, None] # -x_k ]) Aeq2 np.hstack([ y_g, # λ*y^g np.zeros((s1, m)), # 投入松弛占位 -np.eye(s1), # s_r^g np.zeros((s1, s2)), # 非期望产出松弛占位 -y_g[:, i, None] # -y_k^g ]) Aeq3 np.hstack([ y_b, # λ*y^b np.zeros((s2, m)), # 投入松弛占位 np.zeros((s2, s1)), # 期望产出松弛占位 np.eye(s2), # s_r^b -y_b[:, i, None] # -y_k^b ]) Aeq4 np.hstack([ np.zeros(n), # λ np.zeros(m), # s_i^- 1/((s1s2)*y_g[:, i]), # s_r^g系数 1/((s1s2)*y_b[:, i]), # s_r^b系数 np.array([1]) # ρ ]).reshape(1, -1) Aeq np.vstack([Aeq1, Aeq2, Aeq3, Aeq4]) beq np.concatenate([ np.zeros(m s1 s2), # 前三个约束 np.array([1]) # 最后一个约束 ]) bounds [(0, None)] * (n m s1 s2 1) result optimize.linprog(cc, A_eqAeq, b_eqbeq, boundsbounds) theta.append(result.fun)3.2 关键实现技巧矩阵拼接技巧使用np.hstack进行水平拼接使用np.vstack进行垂直拼接单列向量需用[:, None]确保二维结构稀疏矩阵优化 对于大规模问题可转换为稀疏矩阵from scipy.sparse import csr_matrix # 将Aeq转换为稀疏矩阵 Aeq_sparse csr_matrix(Aeq) result optimize.linprog(cc, A_eqAeq_sparse, b_eqbeq, boundsbounds)并行计算加速 使用joblib并行化DMU循环from joblib import Parallel, delayed def solve_dmu(i): # 构建并求解单个DMU的问题 return result.fun theta Parallel(n_jobs4)(delayed(solve_dmu)(i) for i in range(n))4. 结果分析与可视化4.1 效率排名解读计算完成后我们可以得到各DMU的超效率值efficiency_results pd.DataFrame({ City: data.index, SuperEfficiency: theta }).sort_values(SuperEfficiency, ascendingFalse)典型结果特征效率值1超高效单元效率值1有效但非超高效效率值1无效单元4.2 改进方向分析对于无效单元可以从松弛变量找出改进空间# 获取最终解的所有变量 full_solution result.x # 提取关键部分 lambda_weights full_solution[:n] # λ权重 input_slacks full_solution[n:nm] # 投入松弛 good_output_slacks full_solution[nm:nms1] # 期望产出松弛 bad_output_slacks full_solution[nms1:nms1s2] # 非期望产出松弛改进建议表改进类型计算方式政策含义投入过剩s_i^- / x_ik可减少的投入比例期望产出不足s_r^g / y_rk^g可增加的产出比例非期望产出过多s_r^b / y_rk^b需减少的负面产出4.3 可视化呈现使用matplotlib绘制效率分布import matplotlib.pyplot as plt plt.figure(figsize(10, 6)) efficiency_results[SuperEfficiency].plot(kindbarh) plt.axvline(x1, colorr, linestyle--) plt.xlabel(Super-efficiency Score) plt.title(Government Data Openness Efficiency Ranking) plt.tight_layout() plt.show()对于地理数据可使用geopandas绘制空间分布import geopandas as gpd # 假设有城市边界shp文件 cities gpd.read_file(city_boundaries.shp) results_geo cities.merge(efficiency_results, onCity) fig, ax plt.subplots(1, 1, figsize(12, 8)) results_geo.plot(columnSuperEfficiency, legendTrue, schemequantiles, cmapRdYlGn, axax) ax.set_title(Spatial Distribution of Efficiency Scores) plt.show()5. 高级应用与扩展5.1 窗口分析技术对于面板数据可以引入窗口分析window_size 3 # 3年窗口 for year in range(len(years) - window_size 1): window_data data[(data[Year] years[year]) (data[Year] years[year window_size - 1])] # 对每个窗口计算效率5.2 网络SBM模型当存在多阶段生产过程时可扩展为网络SBM# 阶段1投入产出 x1 data[:, 0:2].T y1 data[:, 2:4].T # 阶段2投入产出 x2 data[:, 4:6].T y2 data[:, 6:8].T # 构建网络约束矩阵 Aeq_network np.block([ [Aeq1_stage1, np.zeros_like(Aeq1_stage1)], [np.zeros_like(Aeq1_stage2), Aeq1_stage2], [linkage_constraints] # 阶段间连接约束 ])5.3 敏感性分析评估指标选择对结果的影响sensitivity_results {} for indicator in optional_indicators: modified_data data.drop(indicator, axis1) # 重新计算效率 sensitivity_results[indicator] calculate_efficiency(modified_data)6. 实际应用中的挑战与解决方案6.1 数据质量问题常见问题及处理方法问题类型检测方法解决方案零值数据describe()查看min微小值替换(如1e-6)极端离群值箱线图分析Winsorize处理指标负相关相关系数矩阵重新审视指标合理性# Winsorize处理示例 from scipy.stats import mstats winsorized mstats.winsorize(data, limits[0.05, 0.05])6.2 模型选择验证验证SBM模型适用性相关性检验from scipy.stats import spearmanr # 比较SBM与传统DEA结果 corr, p_value spearmanr(sbm_scores, dea_scores)超效率必要性检验检查有效DMU数量占比若30%则需要超效率模型6.3 计算效率优化大规模问题加速策略问题分解# 按地区分组并行计算 regions data[Region].unique() regional_results {} for region in regions: regional_data data[data[Region] region] regional_results[region] solve_sbm(regional_data)热启动技巧# 使用上一个DMU的解作为初始值 initial_guess None for i in range(n): result optimize.linprog(..., x0initial_guess) initial_guess result.xGPU加速# 使用cupy替代numpy import cupy as cp x_gpu cp.asarray(x) # ...后续矩阵运算在GPU执行7. 完整案例省级政府数据开放评估以某省16个地市为例演示端到端分析流程数据收集投入信息化投资额、专职人员数期望产出开放数据集数量、API调用量非期望产出数据安全事件数预处理脚本def preprocess_data(raw_df): # 缺失值处理 df raw_df.fillna(raw_df.median()) # 非期望产出处理 df[security_incidents] 1 / (df[security_incidents] 1) # 标准化 return (df - df.mean()) / df.std()效率计算# 使用前文实现的SBM求解器 efficiency_scores solve_sbm(preprocessed_data)结果解读效率最高城市A市(1.25)效率最低城市B市(0.68)改进建议B市应重点减少信息化投资冗余(松弛量23%)提高API开放程度(可增加35%)政策模拟# 模拟投入增加10%的效果 modified_input data[investment] * 1.1 new_efficiency solve_sbm(modified_input, data[outputs])8. 前沿扩展方向8.1 结合机器学习方法效率预测模型from sklearn.ensemble import RandomForestRegressor # 使用城市特征预测效率值 X data[[GDP, Population, Education]] y efficiency_scores model RandomForestRegressor().fit(X, y)异常DMU检测from sklearn.svm import OneClassSVM # 检测异常效率模式 clf OneClassSVM().fit(efficiency_features) anomalies clf.predict(efficiency_features)8.2 动态SBM模型引入Malmquist指数分析效率变化def malmquist_index(data_year1, data_year2): eff1 solve_sbm(data_year1) eff2 solve_sbm(data_year2) # 计算技术效率变化和技术进步 return mi_results8.3 空间SBM模型考虑地理空间效应# 构建空间权重矩阵 from libpysal.weights import Queen w Queen.from_dataframe(geo_data) spatial_lag weights.spatial_lag.lag_spatial(w, efficiency_scores)
超效率SBM模型Python实战：用scipy.optimize处理含非期望产出的政府数据效率排名
发布时间：2026/5/24 2:53:32
超效率SBM模型Python实战用scipy.optimize处理含非期望产出的政府数据效率排名在公共管理领域政府开放数据的效率评估一直是学术界和实务界关注的焦点。传统的DEA模型虽然被广泛应用但在处理包含非期望产出如污染排放、能源消耗的复杂场景时往往显得力不从心。超效率SBMSlack-Based Measure模型作为一种非径向、非角度的效率评估方法不仅能够有效处理期望产出与非期望产出并存的情况还能对效率值为1的决策单元进行进一步区分排序为政策制定者提供更精细的决策依据。本文将带您深入超效率SBM模型的数学原理和Python实现细节使用scipy.optimize从零构建完整的求解流程。不同于简单的库调用我们将重点关注如何将数学模型精确转化为线性规划问题处理多类型产出数据时的reshape和拼接技巧约束条件矩阵(Aeq)的构建逻辑与边界条件(bounds)设置实际政府数据评估中的常见陷阱与解决方案1. 超效率SBM模型的核心原理1.1 从传统DEA到SBM的演进传统DEA模型如CCR、BCC存在两个主要局限径向性限制要求所有投入或产出按相同比例增减无法区分有效单元当多个DMU效率值为1时无法进一步排序超效率SBM模型通过引入松弛变量解决了这些问题# 关键改进对比 traditional_dea { 径向性: 是, 效率值范围: [0,1], 非期望产出处理: 困难 } super_sbm { 径向性: 否, 效率值范围: [0,∞), 非期望产出处理: 原生支持 }1.2 数学建模要点考虑有n个决策单元(DMU)每个DMU有m种投入x ∈ R^ms1种期望产出y^g ∈ R^s1s2种非期望产出y^b ∈ R^s2超效率SBM模型的目标函数为min ρ (1 - (1/m)∑(s_i^-/x_ik)) / (1 (1/(s1s2))(∑(s_r^g/y_rk^g) ∑(s_r^b/y_rk^b)))其中s_i^- 为投入松弛量s_r^g 为期望产出松弛量s_r^b 为非期望产出松弛量注意非期望产出的处理需要特别注意方向性其松弛量应被视为需要减少的部分2. 数据准备与预处理2.1 政府开放数据的典型结构以16个城市的政府数据为例常见指标包括指标类型具体指标示例投入指标财政投入、人力投入、IT基础设施期望产出数据开放数量、API调用次数、公众满意度非期望产出数据泄露事件、响应延迟、能源消耗import pandas as pd import numpy as np # 模拟数据加载 data pd.read_excel(government_data.xlsx, index_col0) raw_values data.values # 数据拆分 x raw_values[:, 0:2].T # 2个投入指标 y_g raw_values[:, 2:4].T # 2个期望产出 y_b raw_values[:, 4:5].T # 1个非期望产出2.2 非期望产出的特殊处理非期望产出需要反向处理常见方法包括取倒数法1/y_b线性变换max(y_b) 1 - y_b方向性距离函数直接纳入模型约束# 方法1取倒数处理 y_b_processed 1 / y_b # 方法2线性变换 y_b_max np.max(y_b) y_b_processed y_b_max 1 - y_b提示实际应用中建议测试不同处理方法对结果的影响3. scipy.optimize实现详解3.1 线性规划问题构建将超效率SBM转化为标准LP问题min c^T * v s.t. Aeq * v beq v 0其中决策变量v包含n个λ权重变量m个投入松弛量s1个期望产出松弛量s2个非期望产出松弛量1个ρ效率变量from scipy import optimize m, n x.shape s1 y_g.shape[0] s2 y_b.shape[0] theta [] for i in range(n): # 对每个DMU单独求解 # 目标函数系数 c np.concatenate([ np.zeros(n), # λ -1/(m * x[:, i]), # s_i^- (投入松弛) np.zeros(s1 s2), # s_r^g, s_r^b (产出松弛) np.array([1]) # ρ ]) # 等式约束矩阵 Aeq1 np.hstack([ x, # λ*x np.eye(m), # s_i^- np.zeros((m, s1 s2)), # 产出松弛占位 -x[:, i, None] # -x_k ]) Aeq2 np.hstack([ y_g, # λ*y^g np.zeros((s1, m)), # 投入松弛占位 -np.eye(s1), # s_r^g np.zeros((s1, s2)), # 非期望产出松弛占位 -y_g[:, i, None] # -y_k^g ]) Aeq3 np.hstack([ y_b, # λ*y^b np.zeros((s2, m)), # 投入松弛占位 np.zeros((s2, s1)), # 期望产出松弛占位 np.eye(s2), # s_r^b -y_b[:, i, None] # -y_k^b ]) Aeq4 np.hstack([ np.zeros(n), # λ np.zeros(m), # s_i^- 1/((s1s2)*y_g[:, i]), # s_r^g系数 1/((s1s2)*y_b[:, i]), # s_r^b系数 np.array([1]) # ρ ]).reshape(1, -1) Aeq np.vstack([Aeq1, Aeq2, Aeq3, Aeq4]) beq np.concatenate([ np.zeros(m s1 s2), # 前三个约束 np.array([1]) # 最后一个约束 ]) bounds [(0, None)] * (n m s1 s2 1) result optimize.linprog(cc, A_eqAeq, b_eqbeq, boundsbounds) theta.append(result.fun)3.2 关键实现技巧矩阵拼接技巧使用np.hstack进行水平拼接使用np.vstack进行垂直拼接单列向量需用[:, None]确保二维结构稀疏矩阵优化 对于大规模问题可转换为稀疏矩阵from scipy.sparse import csr_matrix # 将Aeq转换为稀疏矩阵 Aeq_sparse csr_matrix(Aeq) result optimize.linprog(cc, A_eqAeq_sparse, b_eqbeq, boundsbounds)并行计算加速 使用joblib并行化DMU循环from joblib import Parallel, delayed def solve_dmu(i): # 构建并求解单个DMU的问题 return result.fun theta Parallel(n_jobs4)(delayed(solve_dmu)(i) for i in range(n))4. 结果分析与可视化4.1 效率排名解读计算完成后我们可以得到各DMU的超效率值efficiency_results pd.DataFrame({ City: data.index, SuperEfficiency: theta }).sort_values(SuperEfficiency, ascendingFalse)典型结果特征效率值1超高效单元效率值1有效但非超高效效率值1无效单元4.2 改进方向分析对于无效单元可以从松弛变量找出改进空间# 获取最终解的所有变量 full_solution result.x # 提取关键部分 lambda_weights full_solution[:n] # λ权重 input_slacks full_solution[n:nm] # 投入松弛 good_output_slacks full_solution[nm:nms1] # 期望产出松弛 bad_output_slacks full_solution[nms1:nms1s2] # 非期望产出松弛改进建议表改进类型计算方式政策含义投入过剩s_i^- / x_ik可减少的投入比例期望产出不足s_r^g / y_rk^g可增加的产出比例非期望产出过多s_r^b / y_rk^b需减少的负面产出4.3 可视化呈现使用matplotlib绘制效率分布import matplotlib.pyplot as plt plt.figure(figsize(10, 6)) efficiency_results[SuperEfficiency].plot(kindbarh) plt.axvline(x1, colorr, linestyle--) plt.xlabel(Super-efficiency Score) plt.title(Government Data Openness Efficiency Ranking) plt.tight_layout() plt.show()对于地理数据可使用geopandas绘制空间分布import geopandas as gpd # 假设有城市边界shp文件 cities gpd.read_file(city_boundaries.shp) results_geo cities.merge(efficiency_results, onCity) fig, ax plt.subplots(1, 1, figsize(12, 8)) results_geo.plot(columnSuperEfficiency, legendTrue, schemequantiles, cmapRdYlGn, axax) ax.set_title(Spatial Distribution of Efficiency Scores) plt.show()5. 高级应用与扩展5.1 窗口分析技术对于面板数据可以引入窗口分析window_size 3 # 3年窗口 for year in range(len(years) - window_size 1): window_data data[(data[Year] years[year]) (data[Year] years[year window_size - 1])] # 对每个窗口计算效率5.2 网络SBM模型当存在多阶段生产过程时可扩展为网络SBM# 阶段1投入产出 x1 data[:, 0:2].T y1 data[:, 2:4].T # 阶段2投入产出 x2 data[:, 4:6].T y2 data[:, 6:8].T # 构建网络约束矩阵 Aeq_network np.block([ [Aeq1_stage1, np.zeros_like(Aeq1_stage1)], [np.zeros_like(Aeq1_stage2), Aeq1_stage2], [linkage_constraints] # 阶段间连接约束 ])5.3 敏感性分析评估指标选择对结果的影响sensitivity_results {} for indicator in optional_indicators: modified_data data.drop(indicator, axis1) # 重新计算效率 sensitivity_results[indicator] calculate_efficiency(modified_data)6. 实际应用中的挑战与解决方案6.1 数据质量问题常见问题及处理方法问题类型检测方法解决方案零值数据describe()查看min微小值替换(如1e-6)极端离群值箱线图分析Winsorize处理指标负相关相关系数矩阵重新审视指标合理性# Winsorize处理示例 from scipy.stats import mstats winsorized mstats.winsorize(data, limits[0.05, 0.05])6.2 模型选择验证验证SBM模型适用性相关性检验from scipy.stats import spearmanr # 比较SBM与传统DEA结果 corr, p_value spearmanr(sbm_scores, dea_scores)超效率必要性检验检查有效DMU数量占比若30%则需要超效率模型6.3 计算效率优化大规模问题加速策略问题分解# 按地区分组并行计算 regions data[Region].unique() regional_results {} for region in regions: regional_data data[data[Region] region] regional_results[region] solve_sbm(regional_data)热启动技巧# 使用上一个DMU的解作为初始值 initial_guess None for i in range(n): result optimize.linprog(..., x0initial_guess) initial_guess result.xGPU加速# 使用cupy替代numpy import cupy as cp x_gpu cp.asarray(x) # ...后续矩阵运算在GPU执行7. 完整案例省级政府数据开放评估以某省16个地市为例演示端到端分析流程数据收集投入信息化投资额、专职人员数期望产出开放数据集数量、API调用量非期望产出数据安全事件数预处理脚本def preprocess_data(raw_df): # 缺失值处理 df raw_df.fillna(raw_df.median()) # 非期望产出处理 df[security_incidents] 1 / (df[security_incidents] 1) # 标准化 return (df - df.mean()) / df.std()效率计算# 使用前文实现的SBM求解器 efficiency_scores solve_sbm(preprocessed_data)结果解读效率最高城市A市(1.25)效率最低城市B市(0.68)改进建议B市应重点减少信息化投资冗余(松弛量23%)提高API开放程度(可增加35%)政策模拟# 模拟投入增加10%的效果 modified_input data[investment] * 1.1 new_efficiency solve_sbm(modified_input, data[outputs])8. 前沿扩展方向8.1 结合机器学习方法效率预测模型from sklearn.ensemble import RandomForestRegressor # 使用城市特征预测效率值 X data[[GDP, Population, Education]] y efficiency_scores model RandomForestRegressor().fit(X, y)异常DMU检测from sklearn.svm import OneClassSVM # 检测异常效率模式 clf OneClassSVM().fit(efficiency_features) anomalies clf.predict(efficiency_features)8.2 动态SBM模型引入Malmquist指数分析效率变化def malmquist_index(data_year1, data_year2): eff1 solve_sbm(data_year1) eff2 solve_sbm(data_year2) # 计算技术效率变化和技术进步 return mi_results8.3 空间SBM模型考虑地理空间效应# 构建空间权重矩阵 from libpysal.weights import Queen w Queen.from_dataframe(geo_data) spatial_lag weights.spatial_lag.lag_spatial(w, efficiency_scores)
