/****************************/ /* 線形計画法 */ /* coded by Y.Suganuma */ /****************************/ #include <stdio.h> /*************/ /* クラス LP */ /*************/ class LP { int n; // 制約条件の数 int m; // 変数の数 int m_s; // スラック変数の数 int m_a; // 人為変数の数 int mm; // m + m_s + m_a int *a; // 人為変数があるか否か int *cp; // 比較演算子(-1:左辺<右辺, 0:左辺=右辺, 1:左辺>右辺) int *row; // 各行の基底変数の番号 double *z; // 目的関数の係数 double **s; // 単体表 double eps; // 許容誤差 int err; // エラーコード (0:正常終了, 1:解無し) public: /******************************/ /* コンストラクタ */ /* n1 : 制約条件の数 */ /* m1 : 変数の数 */ /* z1 : 目的関数の係数 */ /* eq_l : 制約条件の左辺 */ /* eq_r : 制約条件の右辺 */ /* cp1 : 比較演算子 */ /******************************/ LP(int n1, int m1, double *z1, double **eq_l, double *eq_r, int *cp1) { int i1, i2, k; // 初期設定 eps = 1.0e-10; err = 0; n = n1; m = m1; a = new int [n]; cp = new int [n]; for (i1 = 0; i1 < n; i1++) { a[i1] = 0; cp[i1] = cp1[i1]; } z = new double [m]; for (i1 = 0; i1 < m; i1++) z[i1] = z1[i1]; // スラック変数と人為変数の数を数える m_s = 0; m_a = 0; for (i1 = 0; i1 < n; i1++) { if (cp[i1] == 0) { m_a++; if (eq_r[i1] < 0.0) { eq_r[i1] = -eq_r[i1]; for (i2 = 0; i2 < m; i2++) eq_l[i1][i2] = -eq_l[i1][i2]; } } else { m_s++; if (eq_r[i1] < 0.0) { cp[i1] = -cp[i1]; eq_r[i1] = -eq_r[i1]; for (i2 = 0; i2 < m; i2++) eq_l[i1][i2] = -eq_l[i1][i2]; } if (cp[i1] > 0) m_a++; } } // 単体表の作成 // 初期設定 mm = m + m_s + m_a; row = new int [n]; s = new double * [n+1]; for (i1 = 0; i1 <= n; i1++) { s[i1] = new double [mm+1]; if (i1 < n) { s[i1][0] = eq_r[i1]; for (i2 = 0; i2 < m; i2++) s[i1][i2+1] = eq_l[i1][i2]; for (i2 = m+1; i2 <= mm; i2++) s[i1][i2] = 0.0; } else { for (i2 = 0; i2 <= mm; i2++) s[i1][i2] = 0.0; } } // スラック変数 k = m + 1; for (i1 = 0; i1 < n; i1++) { if (cp[i1] != 0) { if (cp[i1] < 0) { s[i1][k] = 1.0; row[i1] = k - 1; } else s[i1][k] = -1.0; k++; } } // 人為変数 for (i1 = 0; i1 < n; i1++) { if (cp[i1] >= 0) { s[i1][k] = 1.0; row[i1] = k - 1; a[i1] = 1; k++; } } // 目的関数 if (m_a == 0) { for (i1 = 0; i1 < m; i1++) s[n][i1+1] = -z[i1]; } else { for (i1 = 0; i1 <= m+m_s; i1++) { for (i2 = 0; i2 < n; i2++) { if (a[i2] > 0) s[n][i1] -= s[i2][i1]; } } } } /*******************************/ /* 最適化 */ /* sw : =0 : 中間結果無し */ /* =1 : 中間結果あり */ /* return : =0 : 正常終了 */ /* : =1 : 解無し */ /*******************************/ int optimize(int sw) { int i1, i2, k; // フェーズ1 if (sw > 0) { if (m_a > 0) printf("\nphase 1\n"); else printf("\nphase 2\n"); } opt_run(sw); // フェーズ2 if (err == 0 && m_a > 0) { // 目的関数の変更 mm -= m_a; for (i1 = 0; i1 <= mm; i1++) s[n][i1] = 0.0; for (i1 = 0; i1 < n; i1++) { k = row[i1]; if (k < m) s[n][0] += z[k] * s[i1][0]; } for (i1 = 0; i1 < mm; i1++) { for (i2 = 0; i2 < n; i2++) { k = row[i2]; if (k < m) s[n][i1+1] += z[k] * s[i2][i1+1]; } if (i1 < m) s[n][i1+1] -= z[i1]; } // 最適化 if (sw > 0) printf("\nphase 2\n"); opt_run(sw); } return err; } /*******************************/ /* 最適化(単体表の変形) */ /* sw : =0 : 中間結果無し */ /* =1 : 中間結果あり */ /*******************************/ void opt_run(int sw) { int i1, i2, p, q, k; double x, min; err = -1; while (err < 0) { // 中間結果 if (sw > 0) { printf("\n"); output(); } // 列の選択(巡回を防ぐため必ずしも最小値を選択しない,Bland の規則) q = -1; for (i1 = 1; i1 <= mm && q < 0; i1++) { if (s[n][i1] < -eps) q = i1 - 1; } // 終了(最適解) if (q < 0) err = 0; // 行の選択( Bland の規則を採用) else { p = -1; k = -1; min = 0.0; for (i1 = 0; i1 < n; i1++) { if (s[i1][q+1] > eps) { x = s[i1][0] / s[i1][q+1]; if (p < 0 || x < min || x == min && row[i1] < k) { min = x; p = i1; k = row[i1]; } } } // 解無し if (p < 0) err = 1; // 変形 else { x = s[p][q+1]; row[p] = q; for (i1 = 0; i1 <= mm; i1++) s[p][i1] /= x; for (i1 = 0; i1 <= n; i1++) { if (i1 != p) { x = s[i1][q+1]; for (i2 = 0; i2 <= mm; i2++) s[i1][i2] -= x * s[p][i2]; } } } } } } /****************/ /* 単体表の出力 */ /****************/ void output() { int i1, i2; for (i1 = 0; i1 <= n; i1++) { if (i1 < n) printf("x%d", row[i1]+1); else printf(" z"); for (i2 = 0; i2 <= mm; i2++) printf(" %f", s[i1][i2]); printf("\n"); } } /**************/ /* 結果の出力 */ /**************/ void result() { double x; int i1, i2; if (err > 0) printf("\n解が存在しません\n"); else { printf("\n("); for (i1 = 0; i1 < m; i1++) { x = 0.0; for (i2 = 0; i2 < n; i2++) { if (row[i2] == i1) { x = s[i2][0]; break; } } if (i1 == 0) printf("%f", x); else printf(", %f", x); } printf(") のとき,最大値 %f\n", s[n][0]); } } }; /****************/ /* main program */ /****************/ int main() { double **eq_l, *eq_r, *z; int i1, i2, n, m, *cp, sw; char c[2]; // 入力 // 変数の数と式の数 scanf("%d %d", &m, &n); cp = new int [n]; z = new double [m]; eq_l = new double * [n]; eq_r = new double [n]; for (i1 = 0; i1 < n; i1++) eq_l[i1] = new double [m]; // 目的関数の係数 for (i1 = 0; i1 < m; i1++) scanf("%lf", &z[i1]); // 制約条件 for (i1 = 0; i1 < n; i1++) { for (i2 = 0; i2 < m; i2++) scanf("%lf", &eq_l[i1][i2]); scanf("%s", c); if (c[0] == '<') cp[i1] = -1; else if (c[0] == '>') cp[i1] = 1; else cp[i1] = 0; scanf("%lf", &eq_r[i1]); } // 実行 LP lp = LP(n, m, z, eq_l, eq_r, cp); sw = lp.optimize(1); // 結果の出力 lp.result(); return 0; }