▶ 书中第六章部分程序，加上自己补充的代码，包括单纯形法求解线性规划问题
● 单纯形法求解线性规划问题
// 表上作业法，I 为单位阵，y 为对偶变量，z 为目标函数值
// n m 1
// ┌───────────┬───────┬───┐
// │ │ │ │
// m │ A │ I │ b │
// a[m+1][n+m+1] = │ │ │ │
// ├───────────┼───────┼───┤
// 1 │ c │ y │ z │
// └───────────┴───────┴───┘ package package01; import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom; public class class01
{
private static final double EPSILON = 1.0E-10;
private double[][] a; // 工作表
private int m; // 约束数
private int n; // 变量数
private int[] basis; // basis[p] = q，第 p 行的基变量是 x[q] public class01(double[][] A, double[] b, double[] c)
{
m = b.length;
n = c.length;
for (int i = 0; i < m; i++)
{
if (b[i] < 0) // 要求 b 的分量均不小于 0
throw new IllegalArgumentException("RHS must be nonnegative");
}
a = new double[m + 1][n + m + 1];
for (int i = 0; i < m; i++) // a[0:m-1][0:n-1] = A（两边包含）
{
for (int j = 0; j < n; j++)
a[i][j] = A[i][j];
}
for (int i = 0; i < m; i++) // a[0:m-1][n:n+m-1] = I_m
a[i][n + i] = 1.0;
for (int i = 0; i < n; i++) // a[m][0:n-1] = c
a[m][i] = c[i];
for (int i = 0; i < m; i++) // a[0:m-1][n+m] = b
a[i][m + n] = b[i];
basis = new int[m];
for (int i = 0; i < m; i++)
basis[i] = n + i; for (;;) // 单纯形法求解
{
int q = -1;
for (int i = 0; q == -1 && i < m + n; i++) // 找到可优化的变量对应的列 q，即 c[j] > 0 的最靠前索引
q = (a[m][i] > 0) ? i : q;
if (q == -1) // 优化已完成
break; int p = -1; // 依最小比例规则选择离开向量 p
for (int i = 0; i < m; i++) // 要求 a[i][q] > 0 且要么 p 选第一个，要么选壁纸最大的那个
{
if (a[i][q] > EPSILON && (p == -1 || (a[i][m + n] / a[i][q]) < (a[p][m + n] / a[p][q])))
p = i;
}
if (p == -1) // 无离开向量，无界解
throw new ArithmeticException("Linear program is unbounded"); for (int i = 0; i <= m; i++) // 对其他行使用a[p][q] 进行高斯消元
{
for (int j = 0; j <= m + n; j++)
{
if (i != p && j != q)
a[i][j] -= a[p][j] * a[i][q] / a[p][q];
}
}
for (int i = 0; i <= m; i++) // 消元后其他行清零（消除浮点计算误差）
{
if (i != p)
a[i][q] = 0.0;
}
for (int j = 0; j <= m + n; j++)// 主元所在行归一化
{
if (j != q)
a[p][j] /= a[p][q];
}
a[p][q] = 1.0; basis[p] = q; // 更新基向量
}
assert check(A, b, c); // 检查结果
} private int dantzig() // 搜索 c 中值大于零的、最靠前的项的索引
{
int q = 0;
for (int j = 1; j < m + n; j++)
{
if (a[m][j] > a[m][q])
q = j;
}
return a[m][q] <= 0 ? -1 : q; // c 中各项均小于0，优化已完成
} public double value() // 最优解的目标函数值
{
return -a[m][m + n];
} public double[] primal() // 最优解的自变量值
{
double[] x = new double[n];
for (int i = 0; i < m; i++)
{
if (basis[i] < n)
x[basis[i]] = a[i][m + n];
}
return x;
} public double[] dual() // 最优解的对偶变量值
{
double[] y = new double[m];
for (int i = 0; i < m; i++)
y[i] = -a[m][n + i];
return y;
} private boolean isPrimalFeasible(double[][] A, double[] b) // 解的松弛性
{
double[] x = primal();
for (int j = 0; j < x.length; j++) // 检查最优解分量是否均非负
{
if (x[j] < 0.0)
{
StdOut.println("x[" + j + "] = " + x[j] + " is negative");
return false;
}
}
for (int i = 0; i < m; i++) // 检查约束条件
{
double sum = 0.0;
for (int j = 0; j < n; j++)
sum += A[i][j] * x[j];
if (sum > b[i] + EPSILON)
{
StdOut.println("not primal feasibleb\nb[" + i + "] = " + b[i] + ", sum = " + sum);
return false;
}
}
return true;
} private boolean isDualFeasible(double[][] A, double[] c) // 对偶解的松弛性
{
double[] y = dual();
for (int i = 0; i < y.length; i++) // 检查对偶变量不小于 0
{
if (y[i] < 0.0)
{
StdOut.println("y[" + i + "] = " + y[i] + " is negative");
return false;
}
}
for (int j = 0; j < n; j++) // 检查对偶约束条件 A*y ≥ c
{
double sum = 0.0;
for (int i = 0; i < m; i++)
sum += A[i][j] * y[i];
if (sum < c[j] - EPSILON)
{
StdOut.println("not dual feasible\nc[" + j + "] = " + c[j] + ", sum = " + sum);
return false;
}
}
return true;
} private boolean isOptimal(double[] b, double[] c) // 检查解的最优性 value == c*x == y*b
{
double[] x = primal(), y = dual();
double value = value(), value1 = 0.0;
for (int j = 0; j < x.length; j++)
value1 += c[j] * x[j];
double value2 = 0.0;
for (int i = 0; i < y.length; i++)
value2 += y[i] * b[i];
if (Math.abs(value - value1) > EPSILON || Math.abs(value - value2) > EPSILON)
{
StdOut.println("value = " + value + ", cx = " + value1 + ", yb = " + value2);
return false;
}
return true;
} private boolean check(double[][]A, double[] b, double[] c)
{
return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c) && dantzig() == -1;
} private void show() // 输出工作表
{
StdOut.println("m = " + m);
StdOut.println("n = " + n);
for (int i = 0; i <= m; i++)
{
for (int j = 0; j <= m + n; j++)
StdOut.printf("%7.2f ", a[i][j]);
StdOut.println();
}
StdOut.println("value = " + value());
for (int i = 0; i < m; i++)
{
if (basis[i] < n)
StdOut.println("x_" + basis[i] + " = " + a[i][m + n]);
}
StdOut.println();
} private static void test(double[][] A, double[] b, double[] c) // 测试函数
{
class01 lp;
try // 有解正常输出，无解用异常提前退出
{
lp = new class01(A, b, c);
}
catch (ArithmeticException e)
{
System.out.println(e);
return;
}
StdOut.println("value = " + lp.value());
double[] x = lp.primal();
for (int i = 0; i < x.length; i++)
StdOut.println("x[" + i + "] = " + x[i]);
double[] y = lp.dual();
for (int j = 0; j < y.length; j++)
StdOut.println("y[" + j + "] = " + y[j]);
} private static void test1()
{
double[][] A = { { -1, 1, 0 },{ 1, 4, 0 },{ 2, 1, 0 },{ 3, -4, 0 },{ 0, 0, 1 } };
double[] b = { 5, 45, 27, 24, 4 }, c = { 1, 1, 1 };
test(A, b, c);
} private static void test2() // x[0] = 12, x[1] = 28, value = 800
{
double[][] A = { { 5.0, 15.0 },{ 4.0, 4.0 },{ 35.0, 20.0 } };
double[] b = { 480.0, 160.0, 1190.0 }, c = { 13.0, 23.0 };
test(A, b, c);
} private static void test3() // unbounded
{
double[][] A = { { -2.0, -9.0, 1.0, 9.0 },{ 1.0, 1.0, -1.0, -2.0 } };
double[] b = { 3.0, 2.0 }, c = { 2.0, 3.0, -1.0, -12.0 };
test(A, b, c);
} private static void test4() // x[0] = 1.0, x[1] = 0.0, x[2] = 1.0, x[3] = 0.0, value = 1.0，周期
{
double[][] A = { { 0.5, -5.5, -2.5, 9.0 },{ 0.5, -1.5, -0.5, 1.0 },{ 1.0, 0.0, 0.0, 0.0 } };
double[] b = { 0.0, 0.0, 1.0 }, c = { 10.0, -57.0, -9.0, -24.0 };
test(A, b, c);
} public static void main(String[] args)
{
StdOut.println("----- test 1 --------------------");
test1();
StdOut.println("\n----- test 2 --------------------");
test2();
StdOut.println("\n----- test 3 --------------------");
test3();
StdOut.println("\n----- test 4 --------------------");
test4();
StdOut.println("\n----- test random ---------------"); // 随机测试
int m = Integer.parseInt(args[0]), n = Integer.parseInt(args[1]);
double[][] A = new double[m][n];
double[] b = new double[m], c = new double[n];
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
A[i][j] = StdRandom.uniform(100);
}
for (int i = 0; i < m; i++)
b[i] = StdRandom.uniform(1000);
for (int j = 0; j < n; j++)
c[j] = StdRandom.uniform(1000);
class01 lp = new class01(A, b, c);
test(A, b, c);
}
}