▶ 书中第六章部分程序，加上自己补充的代码，包括单纯形法求解线性规划问题
● 单纯形法求解线性规划问题
1 // 表上作业法，I 为单位阵，y 为对偶变量，z 为目标函数值 2 // n m 1 3 // ┌───────────┬───────┬───┐ 4 // │ │ │ │ 5 // m │ A │ I │ b │ 6 // a[m+1][n+m+1] = │ │ │ │ 7 // ├───────────┼───────┼───┤ 8 // 1 │ c │ y │ z │ 9 // └───────────┴───────┴───┘ 10 11 package package01; 12 13 import edu.princeton.cs.algs4.StdOut; 14 import edu.princeton.cs.algs4.StdRandom; 15 16 public class class01 17 { 18 private static final double EPSILON = 1.0E-10; 19 private double[][] a; // 工作表 20 private int m; // 约束数 21 private int n; // 变量数 22 private int[] basis; // basis[p] = q，第 p 行的基变量是 x[q] 23 24 public class01(double[][] A, double[] b, double[] c) 25 { 26 m = b.length; 27 n = c.length; 28 for (int i = 0; i < m; i++) 29 { 30 if (b[i] < 0) // 要求 b 的分量均不小于 0 31 throw new IllegalArgumentException("RHS must be nonnegative"); 32 } 33 a = new double[m + 1][n + m + 1]; 34 for (int i = 0; i < m; i++) // a[0:m-1][0:n-1] = A（两边包含） 35 { 36 for (int j = 0; j < n; j++) 37 a[i][j] = A[i][j]; 38 } 39 for (int i = 0; i < m; i++) // a[0:m-1][n:n+m-1] = I_m 40 a[i][n + i] = 1.0; 41 for (int i = 0; i < n; i++) // a[m][0:n-1] = c 42 a[m][i] = c[i]; 43 for (int i = 0; i < m; i++) // a[0:m-1][n+m] = b 44 a[i][m + n] = b[i]; 45 basis = new int[m]; 46 for (int i = 0; i < m; i++) 47 basis[i] = n + i; 48 49 for (;;) // 单纯形法求解 50 { 51 int q = -1; 52 for (int i = 0; q == -1 && i < m + n; i++) // 找到可优化的变量对应的列 q，即 c[j] > 0 的最靠前索引 53 q = (a[m][i] > 0) ? i : q; 54 if (q == -1) // 优化已完成 55 break; 56 57 int p = -1; // 依最小比例规则选择离开向量 p 58 for (int i = 0; i < m; i++) // 要求 a[i][q] > 0 且要么 p 选第一个，要么选壁纸最大的那个 59 { 60 if (a[i][q] > EPSILON && (p == -1 || (a[i][m + n] / a[i][q]) < (a[p][m + n] / a[p][q]))) 61 p = i; 62 } 63 if (p == -1) // 无离开向量，无界解 64 throw new ArithmeticException("Linear program is unbounded"); 65 66 for (int i = 0; i <= m; i++) // 对其他行使用a[p][q] 进行高斯消元 67 { 68 for (int j = 0; j <= m + n; j++) 69 { 70 if (i != p && j != q) 71 a[i][j] -= a[p][j] * a[i][q] / a[p][q]; 72 } 73 } 74 for (int i = 0; i <= m; i++) // 消元后其他行清零（消除浮点计算误差） 75 { 76 if (i != p) 77 a[i][q] = 0.0; 78 } 79 for (int j = 0; j <= m + n; j++)// 主元所在行归一化 80 { 81 if (j != q) 82 a[p][j] /= a[p][q]; 83 } 84 a[p][q] = 1.0; 85 86 basis[p] = q; // 更新基向量 87 } 88 assert check(A, b, c); // 检查结果 89 } 90 91 private int dantzig() // 搜索 c 中值大于零的、最靠前的项的索引 92 { 93 int q = 0; 94 for (int j = 1; j < m + n; j++) 95 { 96 if (a[m][j] > a[m][q]) 97 q = j; 98 } 99 return a[m][q] <= 0 ? -1 : q; // c 中各项均小于0，优化已完成 100 } 101 102 public double value() // 最优解的目标函数值 103 { 104 return -a[m][m + n]; 105 } 106 107 public double[] primal() // 最优解的自变量值 108 { 109 double[] x = new double[n]; 110 for (int i = 0; i < m; i++) 111 { 112 if (basis[i] < n) 113 x[basis[i]] = a[i][m + n]; 114 } 115 return x; 116 } 117 118 public double[] dual() // 最优解的对偶变量值 119 { 120 double[] y = new double[m]; 121 for (int i = 0; i < m; i++) 122 y[i] = -a[m][n + i]; 123 return y; 124 } 125 126 private boolean isPrimalFeasible(double[][] A, double[] b) // 解的松弛性 127 { 128 double[] x = primal(); 129 for (int j = 0; j < x.length; j++) // 检查最优解分量是否均非负 130 { 131 if (x[j] < 0.0) 132 { 133 StdOut.println("x[" + j + "] = " + x[j] + " is negative"); 134 return false; 135 } 136 } 137 for (int i = 0; i < m; i++) // 检查约束条件 138 { 139 double sum = 0.0; 140 for (int j = 0; j < n; j++) 141 sum += A[i][j] * x[j]; 142 if (sum > b[i] + EPSILON) 143 { 144 StdOut.println("not primal feasibleb b[" + i + "] = " + b[i] + ", sum = " + sum); 145 return false; 146 } 147 } 148 return true; 149 } 150 151 private boolean isDualFeasible(double[][] A, double[] c) // 对偶解的松弛性 152 { 153 double[] y = dual(); 154 for (int i = 0; i < y.length; i++) // 检查对偶变量不小于 0 155 { 156 if (y[i] < 0.0) 157 { 158 StdOut.println("y[" + i + "] = " + y[i] + " is negative"); 159 return false; 160 } 161 } 162 for (int j = 0; j < n; j++) // 检查对偶约束条件 A*y ≥ c 163 { 164 double sum = 0.0; 165 for (int i = 0; i < m; i++) 166 sum += A[i][j] * y[i]; 167 if (sum < c[j] - EPSILON) 168 { 169 StdOut.println("not dual feasible c[" + j + "] = " + c[j] + ", sum = " + sum); 170 return false; 171 } 172 } 173 return true; 174 } 175 176 private boolean isOptimal(double[] b, double[] c) // 检查解的最优性 value == c*x == y*b 177 { 178 double[] x = primal(), y = dual(); 179 double value = value(), value1 = 0.0; 180 for (int j = 0; j < x.length; j++) 181 value1 += c[j] * x[j]; 182 double value2 = 0.0; 183 for (int i = 0; i < y.length; i++) 184 value2 += y[i] * b[i]; 185 if (Math.abs(value - value1) > EPSILON || Math.abs(value - value2) > EPSILON) 186 { 187 StdOut.println("value = " + value + ", cx = " + value1 + ", yb = " + value2); 188 return false; 189 } 190 return true; 191 } 192 193 private boolean check(double[][]A, double[] b, double[] c) 194 { 195 return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c) && dantzig() == -1; 196 } 197 198 private void show() // 输出工作表 199 { 200 StdOut.println("m = " + m); 201 StdOut.println("n = " + n); 202 for (int i = 0; i <= m; i++) 203 { 204 for (int j = 0; j <= m + n; j++) 205 StdOut.printf("%7.2f ", a[i][j]); 206 StdOut.println(); 207 } 208 StdOut.println("value = " + value()); 209 for (int i = 0; i < m; i++) 210 { 211 if (basis[i] < n) 212 StdOut.println("x_" + basis[i] + " = " + a[i][m + n]); 213 } 214 StdOut.println(); 215 } 216 217 private static void test(double[][] A, double[] b, double[] c) // 测试函数 218 { 219 class01 lp; 220 try // 有解正常输出，无解用异常提前退出 221 { 222 lp = new class01(A, b, c); 223 } 224 catch (ArithmeticException e) 225 { 226 System.out.println(e); 227 return; 228 } 229 StdOut.println("value = " + lp.value()); 230 double[] x = lp.primal(); 231 for (int i = 0; i < x.length; i++) 232 StdOut.println("x[" + i + "] = " + x[i]); 233 double[] y = lp.dual(); 234 for (int j = 0; j < y.length; j++) 235 StdOut.println("y[" + j + "] = " + y[j]); 236 } 237 238 private static void test1() 239 { 240 double[][] A = { { -1, 1, 0 },{ 1, 4, 0 },{ 2, 1, 0 },{ 3, -4, 0 },{ 0, 0, 1 } }; 241 double[] b = { 5, 45, 27, 24, 4 }, c = { 1, 1, 1 }; 242 test(A, b, c); 243 } 244 245 private static void test2() // x[0] = 12, x[1] = 28, value = 800 246 { 247 double[][] A = { { 5.0, 15.0 },{ 4.0, 4.0 },{ 35.0, 20.0 } }; 248 double[] b = { 480.0, 160.0, 1190.0 }, c = { 13.0, 23.0 }; 249 test(A, b, c); 250 } 251 252 private static void test3() // unbounded 253 { 254 double[][] A = { { -2.0, -9.0, 1.0, 9.0 },{ 1.0, 1.0, -1.0, -2.0 } }; 255 double[] b = { 3.0, 2.0 }, c = { 2.0, 3.0, -1.0, -12.0 }; 256 test(A, b, c); 257 } 258 259 260 private static void test4() // x[0] = 1.0, x[1] = 0.0, x[2] = 1.0, x[3] = 0.0, value = 1.0，周期 261 { 262 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 } }; 263 double[] b = { 0.0, 0.0, 1.0 }, c = { 10.0, -57.0, -9.0, -24.0 }; 264 test(A, b, c); 265 } 266 267 public static void main(String[] args) 268 { 269 StdOut.println("----- test 1 --------------------"); 270 test1(); 271 StdOut.println(" ----- test 2 --------------------"); 272 test2(); 273 StdOut.println(" ----- test 3 --------------------"); 274 test3(); 275 StdOut.println(" ----- test 4 --------------------"); 276 test4(); 277 StdOut.println(" ----- test random ---------------"); // 随机测试 278 int m = Integer.parseInt(args[0]), n = Integer.parseInt(args[1]); 279 double[][] A = new double[m][n]; 280 double[] b = new double[m], c = new double[n]; 281 for (int i = 0; i < m; i++) 282 { 283 for (int j = 0; j < n; j++) 284 A[i][j] = StdRandom.uniform(100); 285 } 286 for (int i = 0; i < m; i++) 287 b[i] = StdRandom.uniform(1000); 288 for (int j = 0; j < n; j++) 289 c[j] = StdRandom.uniform(1000); 290 class01 lp = new class01(A, b, c); 291 test(A, b, c); 292 } 293 }