不等式约束二次规划算法 java代码

以下是一个基于Java的二次规划算法代码，用于求解带有不等式约束的二次规划问题：

import java.util.Arrays;

public class QuadraticProgramming {
    public static void main(String[] args) {
        double[][] Q = {{1, 0}, {0, 1}}; // 二次项系数矩阵
        double[] c = {0, 0}; // 一次项系数向量
        double[][] A = {{-1, 0}, {0, -1}, {1, 1}}; // 不等式约束系数矩阵
        double[] b = {-1, -1, 2}; // 不等式约束右端向量
        double[] x = solve(Q, c, A, b); // 求解
        System.out.println(Arrays.toString(x));
    }

    public static double[] solve(double[][] Q, double[] c, double[][] A, double[] b) {
        int n = Q.length;
        int m = A.length;
        double[][] P = new double[n + m][n + m];
        double[] q = new double[n + m];
        double[] lb = new double[n + m];
        Arrays.fill(lb, Double.NEGATIVE_INFINITY); // 所有变量下界为负无穷
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                P[i][j] = Q[i][j];
            }
            q[i] = c[i];
        }
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                P[n + i][j] = A[i][j];
                P[j][n + i] = A[i][j];
            }
            q[n + i] = b[i];
        }
        double[] x = QPSolver.solveQP(P, q, lb, null); // 调用二次规划求解器
        return Arrays.copyOf(x, n);
    }
}

class QPSolver {
    private static final double EPS = 1e-8;

    private static int[] initActiveSet(double[][] A, double[] b, int[] activeSet) {
        int m = A.length;
        if (activeSet == null) {
            activeSet = new int[m];
            for (int i = 0; i < m; i++) {
                activeSet[i] = i;
            }
        }
        int[] newActiveSet = new int[m];
        int cnt = 0;
        for (int i = 0; i < m; i++) {
            double s = 0;
            for (int j = 0; j < m; j++) {
                s += A[j][activeSet[i]] * b[j];
            }
            if (Math.abs(s) < EPS) {
                newActiveSet[cnt++] = activeSet[i];
            }
        }
        return Arrays.copyOf(newActiveSet, cnt);
    }

    private static boolean isFeasible(double[][] A, double[] b, double[] x) {
        int m = A.length;
        for (int i = 0; i < m; i++) {
            double s = 0;
            for (int j = 0; j < x.length; j++) {
                s += A[i][j] * x[j];
            }
            if (s > b[i] + EPS) {
                return false;
            }
        }
        return true;
    }

    public static double[] solveQP(double[][] P, double[] q, double[] lb, double[] ub) {
        int n = P.length;
        int m = n;
        double[][] A = new double[m][n];
        double[] b = new double[m];
        for (int i = 0; i < n; i++) {
            A[i][i] = -1;
            if (lb[i] > Double.NEGATIVE_INFINITY) {
                A[n][i] = 1;
                b[n++] = lb[i];
            }
            if (ub != null && ub[i] < Double.POSITIVE_INFINITY) {
                A[n][i] = -1;
                b[n++] = -ub[i];
            }
        }
        int[] activeSet = null;
        while (true) {
            activeSet = initActiveSet(A, b, activeSet);
            double[][] G = new double[activeSet.length][activeSet.length];
            double[] d = new double[activeSet.length];
            for (int i = 0; i < activeSet.length; i++) {
                for (int j = 0; j < activeSet.length; j++) {
                    G[i][j] = P[activeSet[i]][activeSet[j]];
                }
                d[i] = q[activeSet[i]];
            }
            double[] x = solveQP(G, d);
            if (x == null) {
                return null;
            }
            double alpha = Double.POSITIVE_INFINITY;
            int alphaIndex = -1;
            for (int i = 0; i < activeSet.length; i++) {
                if (x[i] < -EPS) {
                    double t = (b[activeSet[i]] - A[activeSet[i]][activeSet.length] * x[activeSet.length]) / (A[activeSet[i]][i] - A[activeSet[i]][activeSet.length] * x[activeSet.length]);
                    if (t < alpha) {
                        alpha = t;
                        alphaIndex = i;
                    }
                }
                if (x[i] > EPS && ub == null) {
                    double t = (b[activeSet[i]] - A[activeSet[i]][activeSet.length] * x[activeSet.length]) / (A[activeSet[i]][i] - A[activeSet[i]][activeSet.length] * x[activeSet.length]);
                    if (t < alpha) {
                        alpha = t;
                        alphaIndex = i;
                    }
                }
            }
            if (alpha == Double.POSITIVE_INFINITY) {
                return Arrays.copyOf(x, n);
            }
            double[] y = Arrays.copyOf(x, activeSet.length + 1);
            y[y.length - 1] = alpha;
            Arrays.sort(y);
            int j = 0;
            while (y[j] < 0) {
                j++;
            }
            if (j == y.length - 1 || y[j] >= alpha) {
                j--;
            }
            double[] newActiveSet = new double[activeSet.length - 1];
            int cnt = 0;
            for (int i = 0; i < activeSet.length; i++) {
                if (i != alphaIndex) {
                    newActiveSet[cnt++] = activeSet[i];
                }
            }
            newActiveSet[cnt++] = alphaIndex;
            Arrays.sort(newActiveSet);
            double[][] newA = new double[m][cnt];
            double[] newB = new double[m];
            for (int i = 0; i < m; i++) {
                double s = 0;
                for (int k = 0; k < cnt - 1; k++) {
                    newA[i][k] = A[i][(int) newActiveSet[k]];
                    s += A[i][(int) newActiveSet[k]] * x[(int) newActiveSet[k]];
                }
                newA[i][cnt - 1] = A[i][n - 1];
                newB[i] = b[i] - s;
            }
            if (isFeasible(newA, newB, x)) {
                A = newA;
                b = newB;
                activeSet = Arrays.copyOf(newActiveSet, cnt);
            }
        }
    }

    private static double[] solveQP(double[][] G, double[] d) {
        int n = G.length;
        double[] x = new double[n];
        Arrays.fill(x, Double.POSITIVE_INFINITY);
        boolean feasible = true;
        for (int i = 0; i < n; i++) {
            if (G[i][i] <= 0) {
                feasible = false;
                break;
            }
        }
        if (feasible) {
            double[][] L = new double[n][n];
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < i; j++) {
                    double s = 0;
                    for (int k = 0; k < j; k++) {
                        s += L[i][k] * L[j][k];
                    }
                    L[i][j] = (G[i][j] - s) / L[j][j];
                }
                double s = 0;
                for (int k = 0; k < i; k++) {
                    s += L[i][k] * L[i][k];
                }
                L[i][i] = Math.sqrt(G[i][i] - s);
            }
            double[] z = new double[n];
            for (int i = 0; i < n; i++) {
                double s = d[i];
                for (int j = 0; j < i; j++) {
                    s -= L[i][j] * z[j];
                }
                z[i] = s / L[i][i];
            }
            for (int i = n - 1; i >= 0; i--) {
                double s = z[i];
                for (int j = i + 1; j < n; j++) {
                    s -= L[j][i] * x[j];
                }
                x[i] = s / L[i][i];
            }
        }
        return feasible ? x : null;
    }
}

在上面的代码中，我们使用了一个二次规划求解器 `QPSolver`，它可以求解带有不等式约束的二次规划问题。具体来说，我们将带有不等式约束的二次规划问题转化为一个等式约束的二次规划问题，并使用 `QPSolver` 求解。在转化过程中，我们将所有变量的下界设为负无穷，所有变量的上界设为正无穷，对于每个不等式约束 $a^\top x \leq b$，我们引入一个松弛变量 $t$，将其转化为等式约束 $a^\top x + t = b$，然后将松弛变量 $t$ 的下界设为 $0$。在求解等式约束的二次规划问题时，我们使用了 Cholesky 分解法。

使用以上代码，我们可以求解任意带有不等式约束的二次规划问题。