不等式约束二次规划——有效集法

预备知识:有效不等式约束是等式约束

这个其实很好理解,通过以下两张图片就可以很清晰的明白这句画的意思:

黑色箭头是约束的区域,蓝色五角星是是全局最优点。对于左边的图,最优点在不等式范围之内,但是这个最优点有没有这个约束都可以求出来,所以这个约束可以看成无效的约束,也就是加不加这个约束,都可以通过无约束的全局解法求出来。相反,右边的不等式约束范围求出来的最优解在绿色点上,可以看到,这个绿色最优点在不等式约束的等式约束上,所以我们完全可以将不等式约束变成等式约束,然后将其看作等式约束QP问题进行求解。所以有效的不等书约束是等式约束。

总体思路

因为不等式中真正起作用的是等式约束,所以我们的思路是:把这些真正起作用的约束找出来,将不等号换成等号,其他的无效不等式直接抛弃掉。之后将这些人为更改出来的等式约束和原问题的等式约束放在一起,将其从不等式约束二次规划问题变成一个等式约束二次规划问题。详细的解法见:等式约束二次规划——变量消除法和KKT法
可以看到,上述思想最重要的就是寻找有效集,也即有效不等式约束到底是哪些?一般情况我们肯定不可能一下子就全部找到,所以会使用一些算法对有效集进行寻找。
下面以
x 1 2 + x 2 2 s . t .   x 1 + x 2 ≥ 1 x 1 − x 2 ≥ − 5 x 1 ≥ 0 x 2 ≥ 0 x 1 ≤ 2 \begin{aligned} x_1^2&+x_2^2\\ s.t.  x_1&+x_2≥1\\ x_1&-x_2≥-5\\ &x_1≥0\\ &x_2≥0\\ &x_1≤2\\ \end{aligned} x12s.t. x1x1+x22+x21x25x10x20x12

为例,做一个不严格的寻找方法图解:

这里红色方框内部是可行集,黑色五角星是理论上的最优点。

如何寻找有效集

我们可以首先随便选择一个点,这个点能让不等式约束集中的一个或多个约束变成等式。在图中表示就是我们一开始可以先在红色的线上任意选择一个点作为迭代的起始点,这个起始点我们记作 x 0 x_0 x0,有了起始点也就相当于选择了将一些不等式约束变成等式约束,那我们就可以使用KKT法对这个问题求解,这里我们将求解出来的最优点记作 x 0 ∗ x_0^* x0.这里总体分为两种情况,第一种是比较幸运,你选起始点正好是这个新等式QP问题的最优点,即 x 0 ∗ = x 0 x_0^*=x_0 x0=x0,

1. x 0 ∗ = x 0 , λ ≥ 0 x_0^*=x_0,λ≥0 x0=x0λ0

这种情况一般就是天选之子,拿阳寿在求解的那种,选的初始点(图中橙色的原点)正好是五角星的那个点,这样一下子就算出来了,具体来说就是只让下图的红色线变成了等式,其他的都不要了,这样进行KKT法求解之后直接就可以算出来最优解:
x 1 2 + x 2 2 s . t .   x 1 + x 2 = 1 \begin{aligned} x_1^2&+x_2^2\\ s.t.  x_1&+x_2=1\\ \end{aligned} x12s.t. x1+x22+x2=1
这种情况基本遇不到,遇到了好好珍惜活着的时光。

2. x 0 ∗ = x 0 , λ j ≤ 0 x_0^*=x_0,λ_j≤0 x0=x0λj0

虽然选择的初始点是当前的最优点,但是拉格朗日乘子中有小于零的,根据拉格朗日乘子的含义,说明之前将那个约束变成等式效果会变差,那说明放宽小于零的这个约束我们会获得更好的值。这里举下图这个例子:

我们选择的点是 ( 2 , 0 ) (2,0) (2,0),这样能变等式的不等式有两个,也就是我们要求解如下问题:
x 1 2 + x 2 2 s . t .   x 2 = 0 x 1 = 2 \begin{aligned} x_1^2&+x_2^2\\ s.t.  &x_2=0\\ &x_1=2\\ \end{aligned} x12s.t. +x22x2=0x1=2
因为比较简单,这里可以快速推导一下:
x 1 2 + x 2 2 + λ 1 ( x 1 − 2 ) + λ 2 x 2 \begin{aligned} x_1^2&+x_2^2+λ_1(x_1-2)+λ_2x_2 \end{aligned} x12+x22+λ1(x12)+λ2x2
求导的结果为:
x 1 = 2 x 2 = 0 λ 1 = − 4 λ 2 = 0 \begin{aligned} x_1 = 2\\ x_2 = 0\\ λ_1 = -4\\ λ_2 = 0 \end{aligned} x1=2x2=0λ1=4λ2=0
可以看到, λ 1 = − 4 < 0 λ_1 = -4<0 λ1=40,这个时候就需要将第一个等式抛弃:

之后重新求解。

3. x 0 ∗ ≠ x 0 x_0^*≠x_0 x0=x0 x 0 ∗ x_0^* x0在可行域中

第三种情况是一开始选的点不是该问题的最优点,当我们求出最优点的时候,发现这个最优点是在可行域的范围之内的,如图:

橙色1号点是选择的初始点,2号点是使用KKT计算之后的当前最优点,且在原问题的红色框之内。这个时候就可以把 x 0 ∗ x_0^* x0当作初始点,扩大一下有效集,之后再次进行迭代:

4. x 0 ∗ ≠ x 0 x_0^*≠x_0 x0=x0 x 0 ∗ x_0^* x0不在可行域中

如果我们使用KKT法求解出来的最优点不在可行域之内:

如上图的点2,超出了原问题的红色边框。
我们的想法是让他沿着从1到2的方向不变,步长选择不超出可行域的最大步长,也就是下图中的点3:

那这个步长要如何选取呢?这其实和一维搜索的格式一样,我们设方向为 d = x 2 − x 1 d=x_2-x_1 d=x2x1,我们希望寻求
x 3 = x 1 + α d x_3=x_1+αd x3=x1+αd
中的 α α α
从图中我们可以清楚的看到点2不在可行域内,就是点2 ( x 2 , y 2 ) (x_2,y_2) (x2,y2)不满足上图的四个红色线。
讨论一般情况,对于不满足的约束:
a i x + b i ≤ 0 a_ix+b_i≤0 aix+bi0
x 3 x_3 x3带入得到:
a i T ( x 1 + α d ) − b i ≤ 0 α a i T d ≤ b i − a i T x 1 \begin{aligned} a_i^T(x_1+αd)-b_i≤0\\ αa_i^Td≤b_i-a_i^Tx_1\\ \end{aligned} aiT(x1+αd)bi0αaiTdbiaiTx1
a i T d ≠ 0 a_i^Td≠0 aiTd=0的时候,有:
α ≤ b i − a i T x 1 a i T d \begin{aligned} α≤\frac{b_i-a_i^Tx_1}{a_i^Td}\\ \end{aligned} αaiTdbiaiTx1
恒成立,那么就就要小于等于最小的那个:
α = m i n   b i − a i T x 1 a i T d \begin{aligned} α=min  \frac{b_i-a_i^Tx_1}{a_i^Td} \\ \end{aligned} α=min aiTdbiaiTx1

当理解了以上四种问题的处理情况之后,就基本理解了有效集算法。

算法框架

找到一个初始可行解x_0
取I(x_0)为x_0的有效集
k = 0
进入for循环:
step1:	求解原等式约束加有效集约束的子优化问题,获得最优解x_1,乘子为λ
step2:	if x_1=x_0:
			if λ≥0:
				x_opt = x_1
				break
			else:
				计算λ_min := min{λ_i,i∈I(x_0)}
				将计算出来的值重新作为初始变量x_0 = x_1
				将此约束从有效集中删除I(x_0) = I(x_0)\{min}
				goto step1
step3:	if x_1≠x_0:
			if x_1是QP问题的可行点:
				将计算出来的值重新作为初始变量x_0 = x_1
				更新有效集I(x_0)
			else:
				d_k := x_1-x_0
				计算α_k = (b_i-a_i^Tx_1}/a_i^Td   i∈I(x_0)\I(x_1)
				x_0 = x_0 +α_kd_k
				更新有效集I(x_0)
			goto step1
							
  • 26
    点赞
  • 66
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论
以下是使用有效集法(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`方法用于计算矩阵的逆。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

薯一个蜂蜜牛奶味的愿

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值