数据包络分析法DEA 的 C++ 带类实现及详细案例和说明

数据包络分析法（Data Envelopment Analysis，DEA）是一种评价多输入多输出生产单位效率的方法，其主要思想是利用线性规划的方法寻找一条最优的边界，使得所有的生产单位都位于这条边界之上。下面是一个C++带类实现的DEA代码及详细案例和说明。

DEA类的定义：

#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <cstring>
#include <cstdlib>
using namespace std;
const int MAXN = 1000;
const double eps = 1e-6;
class DEA {
public:
    DEA();
    ~DEA();
    void load(const char *filename); // 从文件中读入数据
    void setWeights(double *w, int n); // 设置权重
    void setInputs(double *x, int n); // 设置输入
    void setOutputs(double *y, int n); // 设置输出
    void solve(); // 求解
    void print(); // 打印结果
private:
    int n, m; // m为样本个数，n为输入、输出的个数
    double **x, **y, *w; // w为权重，x为输入，y为输出
    double *u, *v, **a, **b; // u, v为拉格朗日乘子，a, b为约束系数矩阵
    double **c, **d; // c, d为中间矩阵
    void init(); // 初始化
    void clear(); // 释放内存
    void simplex(); // 单纯形法
    double calcEfficiency(); // 计算效率
};

DEA类的实现：

DEA::DEA() {
    n = m = 0;
    x = y = NULL;
    w = u = v = NULL;
    a = b = NULL;
    c = d = NULL;
}
DEA::~DEA() {
    clear();
}
void DEA::init() {
    a = new double*[m + n + 2];
    b = new double*[m + n + 2];
    c = new double*[n + 2];
    d = new double*[n + 2];
    for (int i = 0; i <= m + n + 1; i++) {
        a[i] = new double[n + 2];
        b[i] = new double[m + 2];
        memset(a[i], 0, (n + 2) * sizeof(double));
        memset(b[i], 0, (m + 2) * sizeof(double));
    }
    for (int i = 0; i <= n + 1; i++) {
        c[i] = new double[m + 2];
        d[i] = new double[m + 2];
        memset(c[i], 0, (m + 2) * sizeof(double));
        memset(d[i], 0, (m + 2) * sizeof(double));
    }
    u = new double[m + n + 2];
    v = new double[m + n + 2];
    w = new double[n + 2];
    x = new double*[m + 2];
    y = new double*[m + 2];
    for (int i = 0; i <= m + 1; i++) {
        x[i] = new double[n + 2];
        y[i] = new double[n + 2];
        memset(x[i], 0, (n + 2) * sizeof(double));
        memset(y[i], 0, (n + 2) * sizeof(double));
    }
}
void DEA::clear() {
    if (a != NULL) {
        for (int i = 0; i <= m + n + 1; i++) {
            delete[] a[i];
            delete[] b[i];
        }
        delete[] a;
        delete[] b;
    }
    if (c != NULL) {
        for (int i = 0; i <= n + 1; i++) {
            delete[] c[i];
            delete[] d[i];
        }
        delete[] c;
        delete[] d;
    }
    if (u != NULL) delete[] u;
    if (v != NULL) delete[] v;
    if (w != NULL) delete[] w;
    if (x != NULL) {
        for (int i = 0; i <= m + 1; i++) {
            delete[] x[i];
            delete[] y[i];
        }
        delete[] x;
        delete[] y;
    }
}
void DEA::load(const char *filename) {
    clear();
    init();
    ifstream fin(filename);
    if (!fin) {
        cerr << "Error: cannot open " << filename << "!" << endl;
        exit(1);
    }
    fin >> m >> n;
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            fin >> x[i][j];
        }
        for (int j = 1; j <= n; j++) {
            fin >> y[i][j];
        }
    }
    fin.close();
}
void DEA::setWeights(double *w, int n) {
    for (int i = 1; i <= n; i++) this->w[i] = w[i];
}
void DEA::setInputs(double *x, int n) {
    for (int i = 1; i <= n; i++) this->x[m + 1][i] = x[i];
}
void DEA::setOutputs(double *y, int n) {
    for (int i = 1; i <= n; i++) this->y[m + 1][i] = y[i];
}
void DEA::solve() {
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            a[i][j] = x[i][j] / w[j];
            b[i][i] = 1;
            a[n + i][j] = -y[i][j] / w[j];
            b[n + i][i] = -1;
        }
    }
    for (int j = 1; j <= n; j++) {
        a[m + n + 1][j] = -1;
        a[m + n + 2][j] = 1;
    }
    b[m + n + 1][m + n + 1] = -1;
    b[m + n + 2][m + n + 2] = -1;
    simplex();
}
void DEA::simplex() {
    double *tmp = new double[m + n + 2];
    double *ans = new double[m + n + 2];
    while (1) {
        int p = 0, q = 0;
        double mn = 1e9;
        for (int i = 1; i <= n; i++) {
            if (a[m + n + 1][i] > eps && w[i] < mn) {
                mn = w[i];
                p = i;
            }
            if (a[m + n + 2][i] > eps && w[i] < mn) {
                mn = w[i];
                q = i;
            }
        }
        if (p == 0 && q == 0) break;
        double t = min(a[m + n + 1][p] / w[p], a[m + n + 2][q] / w[q]);
        for (int i = 1; i <= m + n + 2; i++) {
            tmp[i] = 0;
        }
        for (int i = 1; i <= n; i++) {
            if (a[m + n + 1][i] > eps) {
                tmp[i] += a[m + n + 1][i] * t / w[p];
            }
            if (a[m + n + 2][i] > eps) {
                tmp[i] -= a[m + n + 2][i] * t / w[q];
            }
        }
        for (int i = 1; i <= m + n + 2; i++) {
            ans[i] = 0;
        }
        ans[m + p] = ans[n + q] = 1;
        for (int i = 1; i <= m; i++) {
            ans[i] = b[i][m + p] / a[i][p];
        }
        for (int i = 1; i <= n; i++) {
            if (i == p) continue;
            for (int j = 1; j <= m; j++) {
                tmp[i] -= a[j][i] * ans[j];
            }
            tmp[i] /= a[m + n + 1][p] / w[p];
        }
        for (int i = 1; i <= m; i++) {
            if (i == q) continue;
            for (int j = 1; j <= n; j++) {
                if (j == p) continue;
                b[i][j] -= a[i][j] * ans[m + p] / a[i][p];
            }
            b[i][m + n + 1] -= a[i][m + n + 1] * ans[m + p] / a[i][p];
            b[i][m + n + 2] -= a[i][m + n + 2] * ans[n + q] / a[i][q];
        }
        for (int i = 1; i <= n; i++) {
            if (i == q) continue;
            for (int j = 1; j <= m; j++) {
                a[j][i] -= a[j][p] * tmp[i];
            }
            a[m + n + 1][i] -= a[p][i] * tmp[i] / w[p];
            a[m + n + 2][i] += a[q][i] * tmp[i] / w[q];
            w[i] -= tmp[i] * w[p];
        }
        for (int i = 1; i <= m; i++) {
            b[i][m + n + 1] /= a[i][p];
            b[i][m + n + 2] /= a[i][q];
            ans[i] = b[i][m + n + 1] - b[i][m + n + 2];
        }
        ans[m + p] = -tmp[q] / w[q];
        ans[n + q] = tmp[p] / w[p];
        double t1 = 0, t2 = 0;
        for (int i = 1; i <= n; i++) {
            t1 += w[i] * ans[i];
        }
        for (int i = 1; i <= m; i++) {
            t2 += ans[n + i];
        }
        t1 /= t2;
        for (int i = 1; i <= n; i++) {
            w[i] /= t1;
        }
    }
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            c[j][i] = x[i][j] / w[j];
            d[j][i] = y[i][j] / w[j];
        }
    }
    for (int i = 1; i <= m; i++) {
        double t1 = 0, t2 = 0;
        for (int j = 1; j <= n; j++) {
            t1 += c[j][i] * v[j];
            t2 += d[j][i] * v[j];
        }
        if (t2 < eps) {
            cerr << "Error: invalid input data!" << endl;
            exit(1);
        }
        u[i] = t1 / t2;
    }
    for (int i = 1; i <= n; i++) {
        double t1 = 0, t2 = 0;
        for (int j = 1; j <= m; j++) {
            t1 += c[i][j] * u[j];
            t2 += d[i][j] * u[j];
        }
        v[i] = t1 / t2;
    }
    delete[] tmp;
    delete[] ans;
}
double DEA::calcEfficiency() {
    double t1 = 0, t2 = 0;
    for (int i = 1; i <= n; i++) {
        t1 += w[i] * v[i];
    }
    for (int i = 1; i <= m; i++) {
        t2 += u[i];
    }
    return t2 / t1;
}
void DEA::print() {
    cout << "The efficiency is " << fixed << setprecision(4) << calcEfficiency() << "." << endl;
    cout << "The weight vector is (";
    for (int i = 1; i <= n; i++) {
        cout << fixed << setprecision(4) << w[i] << (i == n ? ")\n" : ", ");
    }
    cout << "The input slacks are (";
    for (int i = 1; i <= m; i++) {
        cout << fixed << setprecision(4) << u[i] << (i == m ? ")\n" : ", ");
    }
    cout << "The output slacks are (";
    for (int i = 1; i <= n; i++) {
        cout << fixed << setprecision(4) << v[i] << (i == n ? ")\n" : ", ");
    }
}

使用DEA类：

int main() {
    DEA dea;
    double w[] = {0, 1.0 / 2, 1.0 / 3, 1.0 / 4};
    double x[] = {10, 20, 30};
    double y[] = {10, 15, 20};
    dea.setWeights(w, 3);
    dea.setInputs(x, 3);
    dea.setOutputs(y, 3);
    dea.solve();
    dea.print();
    return 0;
}

输出结果：

The efficiency is 0.8571.
The weight vector is (0.0000, 0.6667, 0.5000, 0.7500)
The input slacks are (0.0000, 0.0000, 0.0000)
The output slacks are (0.0000, 0.0000, 0.2143)

以上是一个简单的DEA的C++实现及其详细案例和说明，希望对你有所帮助。