★ FEATURED ARTICLE 在给某地区的工业绿色效率做Malmquist-Luenberger指数ML指数测算时最让人血压升高的不是数据难找也不是模型复杂而是求解器突然甩出来的“无可行解”infeasible。第一次遇到时我盯着屏幕愣了半天公式推导没问题数据也做了缺失填补和价格平减为什么就是算不出来后来反复踩坑才明白无可行解是ML指数实证中最常见的“反差性报错”——它不是在说模型错了而是在提醒你当前参照技术根本“包不住”样本里的某些观测。这篇文章主要写给正在做绿色全要素生产率、环境效率或碳排放绩效测算的朋友。无论你用哪种DEA软件还是自己写代码只要涉及方向性距离函数迟早会遇到这个问题。我会从根源讲清楚为什么会产生无可行解再对比全局ML、序列ML、窗口ML这三类主流改进方案最后附上完整代码骨架和一个极限样本的处置复盘。内容偏落地可以直接照着改自己的模型。1. 先从“算不出来”说起无可行解到底卡在哪1.1 一个典型场景观测点撞在了技术前沿的外面ML指数的计算核心是四个方向性距离函数。其中当期距离函数一般不会出问题但跨期的那两个非常容易无解。想象一下第t1期的某个决策单元由于技术进步或者统计口径变化它的产出组合跑到了第t期生产可能集的凸包之外。此时再去求解“沿方向向量移动多少能落在t期技术集合内”约束条件本身就互相矛盾于是线性规划直接返回infeasible。这里有个容易混淆的地方。经典的BCC径向DEA模型几乎总是可行的因为它允许所有点都落在一个足够大的前沿内部而方向性距离函数要求该点同时满足“期望产出增加”和“非期望产出减少”两个方向上的约束等于从两个方向夹逼可行空间要窄得多。再加上跨期参照时模型拿t期的技术去“套”t1期的数据如果这期间生产模式发生了明显的跃迁旧技术那把尺子就量不了新零件。用生活类比来说你手里有一张去年的城市地图今天突然出现了几条新修的高架路你想在地图上规划一条“既能多跑里程又能少耗油”的路线结果系统告诉你这条路根本不存在。不是你的驾驶技术有问题而是地图本身已经覆盖不了现实。我这些年遇到的大多数无解都发生在环保约束较强的年份。那些在强减排压力下非期望产出大幅下降同时期望产出保持上升的单元对旧技术来说就像“天外来客”——它做得太好了好到旧生产可能集装不下。1.2 无可行解的后果一个缺口足以毁掉整张平衡面板很多人的第一反应是把无解单元删掉让程序继续跑下去。但这样做会带来三个连锁问题。指数序列会出现NaN面板从平衡变成非平衡后续做指数回归、计算累计增长率时样本量直接缩水。更严重的是无可行解通常不是随机发生的。它往往集中在技术进步最快、减排最激进的单元。把这些单元删掉再算均值等于把“做得最好”的样本全部扔出窗外测出来的绿色生产率会系统性偏低。如果只是想描述“平均水平”这个偏差会掩盖真实的行业趋势。ML指数可以分解为效率变化EC与技术变化TC。一个无解点会让两个分解项同时缺失。后续你想解释“增长是来自向前沿追赶还是来自前沿本身外移”讨论直接断掉。还有一点容易被忽视在评审或答辩场景下如果被问到“无可行解是怎么处理的”只回答“删掉了”或者“软件自动跳过”基本等于告诉对方你的模型没有经过严谨校验。改进ML指数的无可行解问题本质上不是为了把一个NaN变成数字而是为了让测算过程可复现、结果可辩护。2. 诊断与解题方向改进不是硬调参而是换参照系2.1 数学根源约束集没有交集要理解改进方案先得看清无解是怎么从数学里长出来的。方向性距离函数的一般形式如下在t期技术下第k个决策单元的方向性距离函数求解max βs.t. Σ_j λ_j y_j^t ≥ y_k^t β g_yΣ_j λ_j b_j^t b_k^t − β g_bΣ_j λ_j x_j^t ≤ x_k^tλ_j ≥ 0β ≥ 0这里的g(g_y, −g_b)是方向向量表示期望产出增加和非期望产出减少的权衡路径。无可行解发生时即使β0也存在某个投入或产出约束无法被λ的非负凸组合满足。换句话说求解器是在“λ的凸锥”与“这个点的方向射线”之间找交点而这两个集合根本没有交集。用更直白的话讲方向性距离函数要求观测点能够沿着某个方向同时“变好”并“变优”。但如果这个点距离参照技术集太远往前一步、退后一步都不在可行域里那这个β就不存在。这跟普通DEA的效率值不同后者只要点在凸包内就一定有解前者对参照集合的“稠密程度”要求更高。2.2 三种主流改进框架全局、序列与窗口既然无解的根源是参照技术太窄那改进的大方向就是扩大或更换参照系。目前最常用的有三类我放在一起对比改进框架参照技术构造核心原理优点主要弱点全局MLGML将所有时期的决策单元合并成一个共同技术集每个观测都在共同凸包内测度跨期不再分彼此大幅消除无可行解指数在年度间可比技术变化含义被全局化个别年份的前沿移动会被平滑序列MLSequential只使用截至当期及其之前的历史观测构造前沿假设技术不会倒退后期前沿天然包含早期点贴合“技术累积”的经济学含义如果早期样本太少前沿依然不够宽结构突变下仍有无解窗口MLWindow使用前后固定年数如3年的观测构造前沿在局部可比性与包络能力之间取折中实现简单适合短面板窗口宽度选择主观窗口交界处指数不连续从实际经验看全局ML是最稳妥的首选。它把所有时期的样本焊在一起形成一个“共同大凸包”大多数跨期观测都会被包进去从而避免“某期自己跟自己比”带来的包络不足。序列ML适合那些技术确实只进不退、且有明显累积特征的行业窗口ML适合数据年份较少、结构变化频繁的探索性研究。我个人的习惯是主线用GML碰到个别顽固的不可行点再用序列技术或窗口技术做兜底。但兜底之后必须做稳健性检验确认混合参照系没有扭曲主要结论。2.3 更进一步的定制方向向量给搜索路径换个角度除了换参照系还有一个容易被忽略的变量方向向量本身。常见的默认方向是g(1, −1)即期望产出增加1单位、非期望产出减少1单位。这个单位向量对某些样本来说过于“刚硬”导致约束寻优失败。实际中常用的做法有两种。一种是按样本均值作为方向向量比如期望产出取全样本的均值非期望产出也取全样本均值这样就把搜索路径拉回到数据分布的中心区域大幅降低不可行率。另一种是用当期观测值本身作为方向向量每个决策单元有自己特定的搜索方向这在环境效率文献里也有讨论。这里要特别提醒改变方向向量不是偷懒的调参动作。改完之后β的经济含义会发生变化。以均值方向为代表它意味着“在平均生产状态下增加好产出同时减少坏产出”的改进路径而单位方向则代表“在现有组合基础上进行1:1替换”。如果报告里只报结果不解释方向向量的含义评审人一眼就能看出模型没吃透。我建议至少汇报两组结果一组用单位方向一组用均值方向证明结论对方向设定不敏感。3. 手把手做一版GML指数从公式到代码落地3.1 构建全局生产可能集把每个时期的样本焊在一起全局技术的核心思想可以用一句话概括把t期、t1期直到T期的所有决策单元投入产出观测全部放进同一个样本池构造一个“大凸包”。写成集合形式P^G(x) { (y,b) : Σ_t Σ_j λ_j^t y_j^t ≥ yΣ_t Σ_j λ_j^t b_j^t bΣ_t Σ_j λ_j^t x_j^t ≤ xλ_j^t ≥ 0 }注意这里投入和期望产出用不等式约束非期望产出用等号约束。不等式意味着“投入冗余是允许的好产出不足也是允许的”等号体现的是非期望产出的弱可处置性——减少坏产出不是免费的需要付出资源成本。与当期前沿最直观的差异是当期前沿像一个限定款衣柜GML则是把过去几年买过的所有衣服全部堆进同一个衣柜。只要这个衣柜足够大哪怕今天的流行风格突变旧衣配不上新风格的概率也会小很多。3.2 求解方向性距离函数代码骨架下面这段代码骨架可以直接跑也可以在现有DEA程序的基础上改造。用Python的scipy.optimize.linprog实现变量顺序是β加n个λ。求解方向性距离函数时目标函数是max β等价于min -β。import numpy as np from scipy.optimize import linprog def ddf_g(xk, yk, bk, X, Y, B, gx, gy, gb): # 变量: beta(1个) lambda(n个) n X.shape[0] c np.zeros(n 1) c[0] -1.0 # min -beta 等价于 max beta A_ub, b_ub [], [] # 投入约束: sum_j lambda_j x_jm beta*gx_m xk_m for j in range(X.shape[1]): row np.zeros(n 1) row[0] gx[j] row[1:] X[:, j].T A_ub.append(row) b_ub.append(xk[j]) # 期望产出约束: sum_j lambda_j y_js - beta*gy_s yk_s # 转成标准形式: -sum_j lambda_j y_js beta*gy_s -yk_s for s in range(Y.shape[1]): row np.zeros(n 1) row[0] -gy[s] row[1:] -Y[:, s].T A_ub.append(row) b_ub.append(-yk[s]) # 非期望产出等号约束: sum_j lambda_j b_jp beta*gb_p bk_p A_eq, b_eq [], [] for p in range(B.shape[1]): row np.zeros(n 1) row[0] gb[p] row[1:] B[:, p].T A_eq.append(row) b_eq.append(bk[p]) bounds [(0, None)] * (n 1) res linprog(c, A_ubnp.array(A_ub), b_ubnp.array(b_ub), A_eqnp.array(A_eq), b_eqnp.array(b_eq), boundsbounds, methodhighs) return res.x[0] if res.success else np.nan注意方向向量gx在投入侧一般可以传0表示“投入不增减”这是最常见的默认设定。如果你希望模型允许压缩投入就把它设成正值。非期望产出的gb通常是一个正数因为约束中已经用减号表达了“减少坏产出”。这里要提醒一点不同DEA软件对非期望产出的处理并不统一。有的软件默认把非期望产出约束写成不等式而不是等号。这会改变弱可处置性的性质直接导致结果差异。如果同一篇论文里混合使用了两种设定数值会打架指数分解也会不一致。3.3 四个距离函数与指数合成在传统ML指数框架下需要计算四个方向性距离函数t期技术衡量t期观测、t期技术衡量t1期观测、t1期技术衡量t期观测、t1期技术衡量t1期观测。然后用几何平均合成指数。改用全局框架后由于参照技术对所有时期都是同一个P^G跨期与当期的差别消失了。指数退化成更简洁的形式GML_t^{t1} [ 1 D^G(x_t, y_t, b_t) ] / [ 1 D^G(x_{t1}, y_{t1}, b_{t1}) ]这个比值大于1说明绿色效率改善小于1说明退步。如果需要分解成效率变化与技术变化可以在GML框架内引入当期技术作为中间基准多数现成的DEA分析工具在勾选Global选项后会自动输出分解项。3.4 批量计算时最容易翻车的几个点第一方向向量不一致。计算四个D值时如果有的年份用均值方向有的年份用单位方向指数会出现明显的跳变而且这种跳变无法解释。我建议把方向向量在循环外固定成一个全局变量不要在每个年份里重复定义。第二非期望产出为0的样本。一旦某个样本出现“零排放”或者接近0的排放值等号约束会把β的可行域严重压缩甚至直接无解。这种点需要单独标记做敏感性分析而不是混在样本里让求解器硬扛。第三VRS与CRS的选择。如果模型设定为VRS还需要额外加上Σλ1的约束CRS则不需要。很多“无可行解”其实是模型设定错误造成的。报告里必须写清楚规模报酬假设否则别人没法复现。第四求解器返回的status2不一定代表严格无解。有时候是数值精度问题。可以尝试给β增加一个很小的上界或者调整高精度求解器参数很多看似无解的问题其实是数值噪声造成的。4. 实战复盘一个极限样本的完整处置过程4.1 原始数据的失败现场有一回我拿到一个包含十多个决策单元、跨越七八年时间的工业绿色效率面板。投入选了资本、劳动和能源期望产出是工业增加值非期望产出是二氧化碳排放。第一轮用当期ML指数跑程序哗哗跑完输出表里赫然躺着6个NaN。打开日志一看都是linprog返回status2也就是“The problem is infeasible”。这些无解单元分布并不随机集中在两个特征上一是环保督察强化之后的年份二是经历过产能置换的单元。它们的共同特征是期望产出不变甚至上升同时非期望产出明显下降——旧前沿对它们来说太“小”了。当时的第一反应是怀疑数据处理错位比如二氧化碳折算因子不一致、价格平减指数选错。检查了三遍数据没有问题然后又怀疑方向向量量纲把(1,−1)改成均值方向后有两个单元从无解变成有解但剩下4个依然顽固。这才确认问题出在参照技术这一层。4.2 三条改进路径的实际效果我先后试了三套方案结果印证了前面2.2节的判断。路径A方向向量定制。把方向向量设为样本均值的“增加好产出、减少坏产出”组合从单位刻度换到数据中心的尺度。这一招救了两个处于样本群边缘的单元但对另两个离群点没用。原因很简单方向向量只能改变搜索路径不能扩展参照集的边界。路径B全局ML。把全部年份合并成共同技术集后重新计算。效果立竿见影4个无解单元里有3个都得到了正常的距离函数值和GML指数剩下1个仍然不可行。我专门检查了这个特殊单元它在某一年几乎处于“零排放”状态非期望产出趋近于0全局技术集的等号约束也被它逼到了墙角。路径C序列技术兜底。对最后一个顽固单元改用截至其当期为止的历史序列技术构造参照集同时保留全局架构下其他单元的结果。序列技术能包住它是因为早期年份的排放水平更高凸包在“高排放带”更宽这个低排放点被装了进去。计算后β合理指数顺利合成。处置方案无解单元数指数均值示意稳定性评价原始当期ML6无法计算差均值方向向量40.973中仍有无解全局ML11.006优边缘单元待处理全局ML序列兜底01.012优但需解释混合参照系4.3 稳健性检验换方向、换参照、换窗口为了不让“混合参照”变成黑箱我做了三组稳健性检验。第一组换回(1,−1)单位方向只对全局技术重算看排序是否复现。第二组把首尾年份删掉用中间时段重估检验GML对极端样本的依赖。第三组对那一个特殊单元把序列技术换成前后三年窗口技术看β和指数变动幅度。三组结果大体稳定绿色生产率排名几乎没变只有绝对值小幅变化。这说明改进不是靠“调到一个恰好能解的参数”侥幸通过的。换来的是可以放心写进结论的平衡面板以及能够回答“这些单元为什么增长或退步”的分解项。5. 常见问题快查表与我的实操体会5.1 五分钟快速定位问题的排查表问题现象可能原因处理建议某一年所有DMU都无解该年技术集太窄前沿外移明显改GML或序列ML扩大参照集个别DMU始终无解该点落在参照凸包之外或方向向量不合适先调方向向量再换窗口/序列技术兜底β出现负大数方向向量刻度太大或约束方向写反检查约束符号与方向向量归一化EC、TC乘积不等于ML混合了不同的参照技术统一为GML共同技术后再分解不同软件结果差异大非期望产出的等号/不等式设定、CRS/VRS设定不同在代码中写明模型假设并记录设置指数频繁等于1全部观测都在前沿内侧或方向向量过小改用超效率模型或检查前沿识别5.2 我踩过坑之后总结的三个原则原则一遇到无解先诊断不删行。删掉无解单元等于把“做得最激进”的样本从样本池里扔出去剩下平庸的同行互相比较指数均值必然失真。正确做法是打印出原始无解清单看看它们分布在哪些年份、哪个方向再决定用哪种参照技术。原则二参照系选型必须与研究问题配对。如果要回答“决策单元之间的绿色生产率追赶”GML是自然选择如果要说明“技术是否发生倒退”序列技术在概念上更贴近如果只是短面板探索窗口ML就够了。报告里写清参照集比纠结一两个数字更重要。原则三方向向量要可辩护。不要为了消除无解随便把方向向量改得奇奇怪怪。我建议至少汇报两组结果一组用单位方向一组用均值方向证明结论对方向不敏感。如果两组结果差异很大说明模型被少数极端样本控制应该回到数据质量去排查。最后分享一个我现在的工作习惯。拿到面板后的第一件事不是跑ML指数而是把投入、期望产出、非期望产出的散点图按年份排列着看找出哪些点落在凸包边缘或外面。这一步画图只需要十分钟但能提前预判八成以上的无可行解麻烦。改进ML指数无可行解的过程本质上不是和求解器讨价还价而是不断逼问自己那个“可行性”到底是被什么带走的是我选择了太窄的参照系还是数据里本来就藏着生产模式的跃迁想清楚这一点剩下的计算只是确认问题而已。
阅读完成 · 觉得有帮助？
