Python实战用scipy.optimize.linprog构建超效率SBM模型评估环境效率当我们需要评估不同地区或企业的环境效率时传统的数据包络分析(DEA)方法往往难以处理非期望产出如污染物排放和效率值超过1的情况。这正是超效率SBM模型的用武之地——它不仅能处理期望产出和非期望产出还能对高效决策单元进行更精确的排序。本文将带你用Python的scipy.optimize.linprog函数从零开始实现这一复杂模型。1. 理解超效率SBM模型的核心概念超效率SBM模型是数据包络分析(DEA)框架下的一个重要扩展它解决了传统DEA模型的两个关键局限非期望产出处理在环境效率评估中我们通常有期望产出如GDP和非期望产出如CO2排放。传统径向DEA模型无法正确处理这种混合产出。效率值上限传统DEA模型的效率值被限制在0-1之间无法区分多个效率值为1的决策单元。模型的核心思想是通过引入松弛变量(slack variables)来捕捉投入过剩和产出不足的情况同时允许效率值超过1。其数学形式可以表示为min ρ (1 - (1/m)Σ(s_i^- / x_io)) / (1 (1/(s1s2))Σ(s_r^ / y_ro) Σ(s_t^b / b_to))其中m是投入指标数量s1和s2分别是期望和非期望产出数量s_i^-是投入松弛变量s_r^和s_t^b是产出松弛变量提示在实际编程实现时我们需要将这个非线性分式规划问题转化为线性规划问题这正是scipy.optimize.linprog的用武之地。2. 数据准备与预处理在开始建模前我们需要准备三类数据投入指标如劳动力、资本、能源消耗等期望产出如GDP、工业产值等非期望产出如CO2排放、废水排放等import numpy as np # 示例数据3个决策单元2个投入1个期望产出1个非期望产出 inputs np.array([ [5, 10], # DMU1的投入 [8, 12], # DMU2的投入 [6, 15] # DMU3的投入 ]) good_outputs np.array([ [20], # DMU1的期望产出 [25], # DMU2的期望产出 [18] # DMU3的期望产出 ]) bad_outputs np.array([ [8], # DMU1的非期望产出 [10], # DMU2的非期望产出 [12] # DMU3的非期望产出 ])数据标准化是重要的一步特别是当不同指标的数值尺度差异很大时def normalize_data(data): means np.mean(data, axis0) stds np.std(data, axis0) return (data - means) / stds norm_inputs normalize_data(inputs) norm_good normalize_data(good_outputs) norm_bad normalize_data(bad_outputs)3. 构建超效率SBM模型的线性规划问题超效率SBM模型的核心是将原始分式规划问题转化为线性规划问题。我们需要为每个决策单元构建以下要素目标函数系数向量c不等式约束矩阵A_ub和向量b_ub等式约束矩阵A_eq和向量b_eqfrom scipy.optimize import linprog def super_sbm(inputs, good_outputs, bad_outputs, dmu_index): num_dmu inputs.shape[0] num_input inputs.shape[1] num_good good_outputs.shape[1] num_bad bad_outputs.shape[1] # 目标函数系数 c np.zeros(1 num_dmu num_input num_good num_bad) c[0] 1 # 对应目标函数中的ρ # 构建约束矩阵 ## 投入约束 input_constraints np.hstack([ -inputs[dmu_index].reshape(-1, 1), # -x_0 inputs.T, # X -np.eye(num_input), # s_i^- np.zeros((num_input, num_good num_bad)) # 产出松弛变量 ]) ## 期望产出约束 good_output_constraints np.hstack([ good_outputs[dmu_index].reshape(-1, 1), # y_0 -good_outputs.T, # -Y np.zeros((num_good, num_input)), # 投入松弛变量 np.eye(num_good), # s_r^ np.zeros((num_good, num_bad)) # 非期望产出松弛变量 ]) ## 非期望产出约束 bad_output_constraints np.hstack([ -bad_outputs[dmu_index].reshape(-1, 1), # -b_0 bad_outputs.T, # B np.zeros((num_bad, num_input num_good)), # 其他松弛变量 np.eye(num_bad) # s_t^b ]) ## 组合所有不等式约束 A_ub np.vstack([input_constraints, good_output_constraints, bad_output_constraints]) b_ub np.zeros(A_ub.shape[0]) ## 等式约束λ之和1 (VRS假设) A_eq np.zeros(1 num_dmu num_input num_good num_bad) A_eq[1:1num_dmu] 1 A_eq A_eq.reshape(1, -1) b_eq np.array([1]) ## 变量边界 bounds [(0, None)] * (1 num_dmu) [(0, None)] * (num_input num_good num_bad) # 求解线性规划 res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds) return res.x[0] # 返回效率值ρ4. 批量计算与结果分析有了单个决策单元的计算函数后我们可以批量计算所有决策单元的效率值def calculate_all_super_sbm(inputs, good_outputs, bad_outputs): num_dmu inputs.shape[0] efficiencies [] for i in range(num_dmu): eff super_sbm(inputs, good_outputs, bad_outputs, i) efficiencies.append(eff) return np.array(efficiencies) # 计算所有DMU的效率值 efficiencies calculate_all_super_sbm(norm_inputs, norm_good, norm_bad) # 打印结果 for i, eff in enumerate(efficiencies): print(fDMU{i1}的超效率SBM效率值: {eff:.4f})结果解读要点效率值1表示该决策单元比生产前沿面更高效效率值1表示位于生产前沿面上效率值1表示存在效率改进空间我们可以进一步将结果可视化import matplotlib.pyplot as plt plt.figure(figsize(10, 6)) plt.bar(range(1, len(efficiencies)1), efficiencies, colorskyblue) plt.axhline(y1, colorr, linestyle--) plt.xlabel(决策单元) plt.ylabel(效率值) plt.title(超效率SBM模型评估结果) plt.grid(axisy, alpha0.3) plt.show()5. 模型优化与实用技巧在实际应用中我们可能会遇到各种问题。以下是几个常见问题的解决方案问题1大规模数据计算速度慢解决方案使用稀疏矩阵存储约束条件并行计算各决策单元的效率值from multiprocessing import Pool def parallel_super_sbm(args): inputs, good_outputs, bad_outputs, i args return super_sbm(inputs, good_outputs, bad_outputs, i) def calculate_parallel(inputs, good_outputs, bad_outputs): num_dmu inputs.shape[0] args [(inputs, good_outputs, bad_outputs, i) for i in range(num_dmu)] with Pool() as p: efficiencies p.map(parallel_super_sbm, args) return np.array(efficiencies)问题2模型无可行解可能原因和解决方法数据存在异常值 → 检查并清洗数据变量尺度差异过大 → 标准化数据约束条件矛盾 → 检查模型设定问题3解释效率驱动因素我们可以通过计算松弛变量来分析效率低下的原因def analyze_slacks(res, num_input, num_good, num_bad): slacks res.x[1num_input:1num_inputnum_goodnum_bad] input_slacks slacks[:num_input] good_slacks slacks[num_input:num_inputnum_good] bad_slacks slacks[num_inputnum_good:] print(投入过剩情况:, input_slacks) print(期望产出不足情况:, good_slacks) print(非期望产出过剩情况:, bad_slacks)6. 实际应用案例中国各省工业碳排放效率评估让我们用一个更实际的案例来演示模型的应用。假设我们有中国30个省份的工业部门数据投入指标工业从业人员(万人)、工业资本存量(亿元)、能源消费量(万吨标煤)期望产出工业增加值(亿元)非期望产出工业CO2排放量(万吨)# 模拟数据 - 实际应用中应从CSV或数据库读取 province_names [f省份{i} for i in range(1, 31)] industrial_labor np.random.uniform(50, 500, 30) industrial_capital np.random.uniform(1000, 10000, 30) energy_consumption np.random.uniform(1000, 8000, 30) industrial_output np.random.uniform(500, 5000, 30) co2_emission industrial_output * np.random.uniform(0.5, 2.0, 30) # 构建输入输出矩阵 inputs np.column_stack([industrial_labor, industrial_capital, energy_consumption]) good_outputs industrial_output.reshape(-1, 1) bad_outputs co2_emission.reshape(-1, 1) # 计算效率 province_efficiencies calculate_all_super_sbm( normalize_data(inputs), normalize_data(good_outputs), normalize_data(bad_outputs) ) # 按效率值排序 ranking np.argsort(-province_efficiencies) print(省份效率排名:) for i, idx in enumerate(ranking): print(f{i1}. {province_names[idx]}: {province_efficiencies[idx]:.3f})结果分析表排名省份效率值主要改进方向1省份X1.215-2省份Y1.103能源利用效率3省份Z0.987资本产出效率............28省份A0.652能源和劳动力效率29省份B0.589全面改进30省份C0.521碳排放强度过高这个排名可以帮助政策制定者识别哪些省份的工业部门效率较高哪些需要改进以及改进的具体方向如减少能源消耗或降低碳排放。
用Python的scipy.optimize.linprog搞定超效率SBM模型：一个环境效率评估的实战案例
Python实战用scipy.optimize.linprog构建超效率SBM模型评估环境效率当我们需要评估不同地区或企业的环境效率时传统的数据包络分析(DEA)方法往往难以处理非期望产出如污染物排放和效率值超过1的情况。这正是超效率SBM模型的用武之地——它不仅能处理期望产出和非期望产出还能对高效决策单元进行更精确的排序。本文将带你用Python的scipy.optimize.linprog函数从零开始实现这一复杂模型。1. 理解超效率SBM模型的核心概念超效率SBM模型是数据包络分析(DEA)框架下的一个重要扩展它解决了传统DEA模型的两个关键局限非期望产出处理在环境效率评估中我们通常有期望产出如GDP和非期望产出如CO2排放。传统径向DEA模型无法正确处理这种混合产出。效率值上限传统DEA模型的效率值被限制在0-1之间无法区分多个效率值为1的决策单元。模型的核心思想是通过引入松弛变量(slack variables)来捕捉投入过剩和产出不足的情况同时允许效率值超过1。其数学形式可以表示为min ρ (1 - (1/m)Σ(s_i^- / x_io)) / (1 (1/(s1s2))Σ(s_r^ / y_ro) Σ(s_t^b / b_to))其中m是投入指标数量s1和s2分别是期望和非期望产出数量s_i^-是投入松弛变量s_r^和s_t^b是产出松弛变量提示在实际编程实现时我们需要将这个非线性分式规划问题转化为线性规划问题这正是scipy.optimize.linprog的用武之地。2. 数据准备与预处理在开始建模前我们需要准备三类数据投入指标如劳动力、资本、能源消耗等期望产出如GDP、工业产值等非期望产出如CO2排放、废水排放等import numpy as np # 示例数据3个决策单元2个投入1个期望产出1个非期望产出 inputs np.array([ [5, 10], # DMU1的投入 [8, 12], # DMU2的投入 [6, 15] # DMU3的投入 ]) good_outputs np.array([ [20], # DMU1的期望产出 [25], # DMU2的期望产出 [18] # DMU3的期望产出 ]) bad_outputs np.array([ [8], # DMU1的非期望产出 [10], # DMU2的非期望产出 [12] # DMU3的非期望产出 ])数据标准化是重要的一步特别是当不同指标的数值尺度差异很大时def normalize_data(data): means np.mean(data, axis0) stds np.std(data, axis0) return (data - means) / stds norm_inputs normalize_data(inputs) norm_good normalize_data(good_outputs) norm_bad normalize_data(bad_outputs)3. 构建超效率SBM模型的线性规划问题超效率SBM模型的核心是将原始分式规划问题转化为线性规划问题。我们需要为每个决策单元构建以下要素目标函数系数向量c不等式约束矩阵A_ub和向量b_ub等式约束矩阵A_eq和向量b_eqfrom scipy.optimize import linprog def super_sbm(inputs, good_outputs, bad_outputs, dmu_index): num_dmu inputs.shape[0] num_input inputs.shape[1] num_good good_outputs.shape[1] num_bad bad_outputs.shape[1] # 目标函数系数 c np.zeros(1 num_dmu num_input num_good num_bad) c[0] 1 # 对应目标函数中的ρ # 构建约束矩阵 ## 投入约束 input_constraints np.hstack([ -inputs[dmu_index].reshape(-1, 1), # -x_0 inputs.T, # X -np.eye(num_input), # s_i^- np.zeros((num_input, num_good num_bad)) # 产出松弛变量 ]) ## 期望产出约束 good_output_constraints np.hstack([ good_outputs[dmu_index].reshape(-1, 1), # y_0 -good_outputs.T, # -Y np.zeros((num_good, num_input)), # 投入松弛变量 np.eye(num_good), # s_r^ np.zeros((num_good, num_bad)) # 非期望产出松弛变量 ]) ## 非期望产出约束 bad_output_constraints np.hstack([ -bad_outputs[dmu_index].reshape(-1, 1), # -b_0 bad_outputs.T, # B np.zeros((num_bad, num_input num_good)), # 其他松弛变量 np.eye(num_bad) # s_t^b ]) ## 组合所有不等式约束 A_ub np.vstack([input_constraints, good_output_constraints, bad_output_constraints]) b_ub np.zeros(A_ub.shape[0]) ## 等式约束λ之和1 (VRS假设) A_eq np.zeros(1 num_dmu num_input num_good num_bad) A_eq[1:1num_dmu] 1 A_eq A_eq.reshape(1, -1) b_eq np.array([1]) ## 变量边界 bounds [(0, None)] * (1 num_dmu) [(0, None)] * (num_input num_good num_bad) # 求解线性规划 res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds) return res.x[0] # 返回效率值ρ4. 批量计算与结果分析有了单个决策单元的计算函数后我们可以批量计算所有决策单元的效率值def calculate_all_super_sbm(inputs, good_outputs, bad_outputs): num_dmu inputs.shape[0] efficiencies [] for i in range(num_dmu): eff super_sbm(inputs, good_outputs, bad_outputs, i) efficiencies.append(eff) return np.array(efficiencies) # 计算所有DMU的效率值 efficiencies calculate_all_super_sbm(norm_inputs, norm_good, norm_bad) # 打印结果 for i, eff in enumerate(efficiencies): print(fDMU{i1}的超效率SBM效率值: {eff:.4f})结果解读要点效率值1表示该决策单元比生产前沿面更高效效率值1表示位于生产前沿面上效率值1表示存在效率改进空间我们可以进一步将结果可视化import matplotlib.pyplot as plt plt.figure(figsize(10, 6)) plt.bar(range(1, len(efficiencies)1), efficiencies, colorskyblue) plt.axhline(y1, colorr, linestyle--) plt.xlabel(决策单元) plt.ylabel(效率值) plt.title(超效率SBM模型评估结果) plt.grid(axisy, alpha0.3) plt.show()5. 模型优化与实用技巧在实际应用中我们可能会遇到各种问题。以下是几个常见问题的解决方案问题1大规模数据计算速度慢解决方案使用稀疏矩阵存储约束条件并行计算各决策单元的效率值from multiprocessing import Pool def parallel_super_sbm(args): inputs, good_outputs, bad_outputs, i args return super_sbm(inputs, good_outputs, bad_outputs, i) def calculate_parallel(inputs, good_outputs, bad_outputs): num_dmu inputs.shape[0] args [(inputs, good_outputs, bad_outputs, i) for i in range(num_dmu)] with Pool() as p: efficiencies p.map(parallel_super_sbm, args) return np.array(efficiencies)问题2模型无可行解可能原因和解决方法数据存在异常值 → 检查并清洗数据变量尺度差异过大 → 标准化数据约束条件矛盾 → 检查模型设定问题3解释效率驱动因素我们可以通过计算松弛变量来分析效率低下的原因def analyze_slacks(res, num_input, num_good, num_bad): slacks res.x[1num_input:1num_inputnum_goodnum_bad] input_slacks slacks[:num_input] good_slacks slacks[num_input:num_inputnum_good] bad_slacks slacks[num_inputnum_good:] print(投入过剩情况:, input_slacks) print(期望产出不足情况:, good_slacks) print(非期望产出过剩情况:, bad_slacks)6. 实际应用案例中国各省工业碳排放效率评估让我们用一个更实际的案例来演示模型的应用。假设我们有中国30个省份的工业部门数据投入指标工业从业人员(万人)、工业资本存量(亿元)、能源消费量(万吨标煤)期望产出工业增加值(亿元)非期望产出工业CO2排放量(万吨)# 模拟数据 - 实际应用中应从CSV或数据库读取 province_names [f省份{i} for i in range(1, 31)] industrial_labor np.random.uniform(50, 500, 30) industrial_capital np.random.uniform(1000, 10000, 30) energy_consumption np.random.uniform(1000, 8000, 30) industrial_output np.random.uniform(500, 5000, 30) co2_emission industrial_output * np.random.uniform(0.5, 2.0, 30) # 构建输入输出矩阵 inputs np.column_stack([industrial_labor, industrial_capital, energy_consumption]) good_outputs industrial_output.reshape(-1, 1) bad_outputs co2_emission.reshape(-1, 1) # 计算效率 province_efficiencies calculate_all_super_sbm( normalize_data(inputs), normalize_data(good_outputs), normalize_data(bad_outputs) ) # 按效率值排序 ranking np.argsort(-province_efficiencies) print(省份效率排名:) for i, idx in enumerate(ranking): print(f{i1}. {province_names[idx]}: {province_efficiencies[idx]:.3f})结果分析表排名省份效率值主要改进方向1省份X1.215-2省份Y1.103能源利用效率3省份Z0.987资本产出效率............28省份A0.652能源和劳动力效率29省份B0.589全面改进30省份C0.521碳排放强度过高这个排名可以帮助政策制定者识别哪些省份的工业部门效率较高哪些需要改进以及改进的具体方向如减少能源消耗或降低碳排放。
