SQP算法 java代码

以下是一个简单的 Java 实现 SQP 算法的代码，供参考：

public class SQPAlgorithm {
    private double epsilon = 0.001; // 精度要求
    private int maxIterations = 100; // 最大迭代次数

    public void setEpsilon(double epsilon) {
        this.epsilon = epsilon;
    }

    public void setMaxIterations(int maxIterations) {
        this.maxIterations = maxIterations;
    }

    public double[] optimize(Function<double[], Double> objectiveFunction, double[] initialGuess) {
        int n = initialGuess.length;
        double[] x = Arrays.copyOf(initialGuess, n); // 当前估计值
        double[] gradient = new double[n]; // 梯度
        double[][] hessian = new double[n][n]; // 黑塞矩阵
        double lambda = 0.1; // 初始步长
        int iteration = 0;
        while (iteration < maxIterations) {
            iteration++;
            double f = objectiveFunction.apply(x);
            gradient = computeGradient(objectiveFunction, x);
            hessian = computeHessian(objectiveFunction, x);
            double[] p = solveQP(hessian, gradient);
            double[] xNew = newtonLineSearch(objectiveFunction, x, p, lambda);
            double fNew = objectiveFunction.apply(xNew);
            double rho = (f - fNew) / (0.5 * p.dot(hessian.dot(p)) + gradient.dot(p)); // 计算鲁棒系数
            if (rho > 0) {
                x = xNew;
                lambda = Math.max(lambda / 3.0, 1e-7); // 缩小步长
            } else {
                lambda = Math.min(lambda * 3.0, 1e7); // 增加步长
            }
            if (p.norm() < epsilon) {
                break; // 满足精度要求
            }
        }
        return x;
    }

    private double[] computeGradient(Function<double[], Double> f, double[] x) {
        int n = x.length;
        double[] gradient = new double[n];
        double f0 = f.apply(x);
        for (int i = 0; i < n; i++) {
            double[] x1 = Arrays.copyOf(x, n);
            x1[i] += epsilon;
            gradient[i] = (f.apply(x1) - f0) / epsilon;
        }
        return gradient;
    }

    private double[][] computeHessian(Function<double[], Double> f, double[] x) {
        int n = x.length;
        double[][] hessian = new double[n][n];
        double f0 = f.apply(x);
        for (int i = 0; i < n; i++) {
            for (int j = i; j < n; j++) {
                double[] x1 = Arrays.copyOf(x, n);
                double[] x2 = Arrays.copyOf(x, n);
                x1[i] += epsilon;
                x1[j] += epsilon;
                x2[i] += epsilon;
                x2[j] -= epsilon;
                double f1 = f.apply(x1);
                double f2 = f.apply(x2);
                hessian[i][j] = hessian[j][i] = (f1 - f2 - (f.apply(x1) - f0)) / (4 * epsilon * epsilon);
            }
        }
        return hessian;
    }

    private double[] solveQP(double[][] hessian, double[] gradient) {
        int n = gradient.length;
        Matrix P = new Matrix(hessian);
        Matrix q = new Matrix(gradient, n);
        double[][] G = new double[2 * n][n];
        double[] h = new double[2 * n];
        for (int i = 0; i < n; i++) {
            G[i][i] = -1;
            h[i] = 0;
        }
        for (int i = n; i < 2 * n; i++) {
            G[i][i - n] = 1;
            h[i] = 1;
        }
        Matrix GMatrix = new Matrix(G);
        Matrix hMatrix = new Matrix(h, 2 * n);
        Solver solver = new Solver(P, q, GMatrix, hMatrix);
        return solver.solve();
    }

    private double[] newtonLineSearch(Function<double[], Double> f, double[] x, double[] p, double lambda) {
        int n = x.length;
        double f0 = f.apply(x);
        double[] gradient = computeGradient(f, x);
        double[] xNew = Arrays.copyOf(x, n);
        double alpha = 1.0;
        while (true) {
            for (int i = 0; i < n; i++) {
                xNew[i] = x[i] + alpha * p[i];
            }
            double fNew = f.apply(xNew);
            if (fNew < f0 + 0.5 * alpha * gradient.dot(p)) {
                break; // Armijo 准则
            }
            alpha *= 0.5; // 缩小步长
        }
        return xNew;
    }
}

其中，`Function<double[], Double>` 表示目标函数，`Solver` 是一个二次规划求解器，`Matrix` 是矩阵类，可以使用 `commons-math3` 库来实现。在 `optimize` 方法中，我们首先初始化一些参数，然后在循环中反复执行以下步骤：

- 计算当前估计值的目标函数值、梯度和黑塞矩阵；
- 解二次规划问题，得到搜索方向；
- 在搜索方向上执行牛顿线搜索，得到新的估计值；
- 根据鲁棒系数更新步长。

当搜索方向的范数小于精度要求时，算法停止迭代并返回当前估计值。

上述实现中的牛顿线搜索采用了 Armijo 准则，而没有采用 Wolfe 准则，这是因为 Armijo 准则更加简单并且通常也能够收敛。如果你对 Wolfe 准则更感兴趣，可以参考其他资料进行实现。