有效集法求解不等式约束二次规划算法的java代码
时间: 2023-06-18 21:02:57 浏览: 64
以下是使用有效集法(Active Set Method)求解不等式约束二次规划问题的Java代码:
```java
public class ActiveSetMethod {
private double[] x; // 最优解
private double f; // 目标函数值
private int[] activeSet; // 有效集
private boolean isOptimal; // 是否为最优解
private boolean isFeasible; // 是否可行
private double[] lambda; // 拉格朗日乘子
private double[] gradient; // 梯度
private double[][] hessian; // 黑塞矩阵
private double[] lowerBound; // 变量下界
private double[] upperBound; // 变量上界
private double eps = 1e-8; // 误差容限
private int maxIter = 1000; // 最大迭代次数
/**
* 构造函数
* @param x0 初始解
* @param hessian 黑塞矩阵
* @param gradient 梯度
* @param lowerBound 变量下界
* @param upperBound 变量上界
*/
public ActiveSetMethod(double[] x0, double[][] hessian, double[] gradient, double[] lowerBound, double[] upperBound) {
this.x = x0;
this.hessian = hessian;
this.gradient = gradient;
this.lowerBound = lowerBound;
this.upperBound = upperBound;
this.activeSet = new int[x.length];
this.lambda = new double[x.length];
this.isOptimal = false;
this.isFeasible = true;
for (int i = 0; i < x.length; i++) {
if (lowerBound[i] > x[i] || upperBound[i] < x[i]) {
isFeasible = false;
break;
}
}
}
/**
* 有效集法求解二次规划问题
*/
public void solve() {
int iter = 0;
while (!isOptimal && iter < maxIter) {
iter++;
// 计算拉格朗日乘子
double[] temp = multiply(hessian, x);
for (int i = 0; i < x.length; i++) {
temp[i] += gradient[i];
}
for (int i = 0; i < x.length; i++) {
double sum = 0.0;
for (int j = 0; j < x.length; j++) {
sum += hessian[i][j] * x[j];
}
if (Math.abs(sum + gradient[i]) < eps && x[i] > lowerBound[i] + eps && x[i] < upperBound[i] - eps) {
lambda[i] = 0.0;
} else {
lambda[i] = Math.max(0.0, -temp[i] / (hessian[i][i] + eps));
}
}
// 判断是否为最优解
double[] temp2 = new double[x.length];
for (int i = 0; i < x.length; i++) {
temp2[i] = temp[i] + lambda[i] * hessian[i][i];
}
double maxLambda = 0.0;
int maxIndex = -1;
for (int i = 0; i < x.length; i++) {
if (lambda[i] > maxLambda) {
maxLambda = lambda[i];
maxIndex = i;
}
}
if (maxLambda < eps) {
isOptimal = true;
f = 0.5 * multiply(x, multiply(hessian, x)) + multiply(x, gradient);
continue;
}
// 计算搜索方向
double[][] subMatrix = new double[activeSet.length][activeSet.length];
double[] subGradient = new double[activeSet.length];
int count = 0;
for (int i = 0; i < x.length; i++) {
if (Math.abs(lambda[i]) < eps) {
activeSet[count] = i;
count++;
}
}
for (int i = 0; i < count; i++) {
for (int j = 0; j < count; j++) {
subMatrix[i][j] = hessian[activeSet[i]][activeSet[j]];
}
subGradient[i] = gradient[activeSet[i]];
}
double[] subX = new double[count];
for (int i = 0; i < count; i++) {
subX[i] = x[activeSet[i]];
}
double[] subLambda = new double[count];
for (int i = 0; i < count; i++) {
subLambda[i] = lambda[activeSet[i]];
}
double[] subDirection = multiply(inverse(subMatrix), subGradient);
// 计算步长
double alpha = 1.0;
if (maxIndex != -1) {
if (temp2[maxIndex] > 0) {
alpha = Math.min(alpha, (upperBound[maxIndex] - x[maxIndex]) / (eps + subDirection[maxIndex]));
} else {
alpha = Math.min(alpha, (lowerBound[maxIndex] - x[maxIndex]) / (eps + subDirection[maxIndex]));
}
}
for (int i = 0; i < count; i++) {
if (subDirection[i] > 0) {
alpha = Math.min(alpha, (upperBound[activeSet[i]] - subX[i]) / (eps + subDirection[i]));
} else {
alpha = Math.min(alpha, (lowerBound[activeSet[i]] - subX[i]) / (eps + subDirection[i]));
}
}
// 更新解
for (int i = 0; i < x.length; i++) {
x[i] += alpha * subDirection[i];
}
// 更新有效集
for (int i = 0; i < x.length; i++) {
if (Math.abs(x[i] - lowerBound[i]) < eps || Math.abs(x[i] - upperBound[i]) < eps) {
activeSet[count] = i;
count++;
}
}
// 判断是否为最优解
double[] temp3 = multiply(hessian, x);
for (int i = 0; i < x.length; i++) {
temp3[i] += gradient[i];
}
for (int i = 0; i < x.length; i++) {
double sum = 0.0;
for (int j = 0; j < x.length; j++) {
sum += hessian[i][j] * x[j];
}
if (Math.abs(sum + gradient[i]) < eps && x[i] > lowerBound[i] + eps && x[i] < upperBound[i] - eps) {
lambda[i] = 0.0;
} else {
lambda[i] = Math.max(0.0, -temp3[i] / (hessian[i][i] + eps));
}
}
}
}
/**
* 获取最优解
* @return 最优解
*/
public double[] getX() {
return x;
}
/**
* 获取目标函数值
* @return 目标函数值
*/
public double getF() {
return f;
}
/**
* 获取有效集
* @return 有效集
*/
public int[] getActiveSet() {
return activeSet;
}
/**
* 获取是否为最优解
* @return 是否为最优解
*/
public boolean isOptimal() {
return isOptimal;
}
/**
* 获取是否可行
* @return 是否可行
*/
public boolean isFeasible() {
return isFeasible;
}
/**
* 获取拉格朗日乘子
* @return 拉格朗日乘子
*/
public double[] getLambda() {
return lambda;
}
/**
* 两个向量相乘
* @param a 向量a
* @param b 向量b
* @return 向量a和b的点积
*/
private double multiply(double[] a, double[] b) {
double sum = 0.0;
for (int i = 0; i < a.length; i++) {
sum += a[i] * b[i];
}
return sum;
}
/**
* 矩阵和向量相乘
* @param a 矩阵a
* @param b 向量b
* @return 矩阵a和向量b的乘积
*/
private double[] multiply(double[][] a, double[] b) {
double[] c = new double[b.length];
for (int i = 0; i < b.length; i++) {
for (int j = 0; j < b.length; j++) {
c[i] += a[i][j] * b[j];
}
}
return c;
}
/**
* 求矩阵的逆
* @param a 矩阵a
* @return 矩阵a的逆
*/
private double[][] inverse(double[][] a) {
int n = a.length;
double[][] b = new double[n][n];
for (int i = 0; i < n; i++) {
b[i][i] = 1.0;
}
for (int k = 0; k < n; k++) {
int maxIndex = k;
double max = Math.abs(a[k][k]);
for (int i = k + 1; i < n; i++) {
if (Math.abs(a[i][k]) > max) {
maxIndex = i;
max = Math.abs(a[i][k]);
}
}
if (maxIndex != k) {
for (int j = 0; j < n; j++) {
double temp = a[k][j];
a[k][j] = a[maxIndex][j];
a[maxIndex][j] = temp;
temp = b[k][j];
b[k][j] = b[maxIndex][j];
b[maxIndex][j] = temp;
}
}
double pivot = a[k][k];
if (Math.abs(pivot) < eps) {
throw new RuntimeException("Matrix is singular.");
}
for (int j = k; j < n; j++) {
a[k][j] /= pivot;
}
for (int j = 0; j < n; j++) {
b[k][j] /= pivot;
}
for (int i = k + 1; i < n; i++) {
double factor = a[i][k];
for (int j = k; j < n; j++) {
a[i][j] -= factor * a[k][j];
}
for (int j = 0; j < n; j++) {
b[i][j] -= factor * b[k][j];
}
}
}
for (int k = n - 1; k > 0; k--) {
for (int i = k - 1; i >= 0; i--) {
double factor = a[i][k];
for (int j = 0; j < n; j++) {
b[i][j] -= factor * b[k][j];
}
}
}
return b;
}
}
```
其中,`x`、`hessian`、`gradient`、`lowerBound`和`upperBound`分别表示初始解、黑塞矩阵、梯度、变量下界和变量上界。`activeSet`表示有效集,`lambda`表示拉格朗日乘子,`isOptimal`表示是否为最优解,`isFeasible`表示是否可行,`eps`表示误差容限,`maxIter`表示最大迭代次数。`solve`方法使用有效集法求解二次规划问题,`getX`方法返回最优解,`getF`方法返回目标函数值,`getActiveSet`方法返回有效集,`isOptimal`方法返回是否为最优解,`isFeasible`方法返回是否可行,`getLambda`方法返回拉格朗日乘子。`multiply`方法用于计算两个向量的点积,`inverse`方法用于计算矩阵的逆。