# -*- coding: UTF-8 -*- from math import * import numpy as np ############################################ # クラス LP(線形計画法) # coded by Y.Suganuma ############################################ class LP : ################################ # コンストラクタ # n1 : 制約条件の数 # m1 : 変数の数 # z1 : 目的関数の係数 # eq_l : 制約条件の左辺 # eq_r : 制約条件の右辺 # cp1 : 比較演算子 ################################ def __init__(self, n1, m1, z1, eq_l, eq_r, cp1) : # 初期設定 self.eps = 1.0e-10 # 許容誤差 self.err = 0 # エラーコード (0:正常終了, 1:解無し) self.n = n1 # 制約条件の数 self.m = m1 # 変数の数 self.a = np.zeros(self.n) # 人為変数があるか否か self.cp = cp1 # 比較演算子(-1:左辺<右辺, 0:左辺=右辺, 1:左辺>右辺) self.z = z1 # 目的関数の係数 # スラック変数と人為変数の数を数える self.m_s = 0 # スラック変数の数 self.m_a = 0 # 人為変数の数 for i1 in range(0, self.n) : if self.cp[i1] == 0 : self.m_a += 1 if eq_r[i1] < 0.0 : eq_r[i1] = -eq_r[i1] for i2 in rang(0, self.m) : eq_l[i1][i2] = -eq_l[i1][i2] else : self.m_s += 1 if eq_r[i1] < 0.0 : self.cp[i1] = -self.cp[i1] eq_r[i1] = -eq_r[i1] for i2 in rang(0, self.m) : eq_l[i1][i2] = -eq_l[i1][i2] if self.cp[i1] > 0 : self.m_a += 1 # 単体表の作成 # 初期設定 self.mm = self.m + self.m_s + self.m_a self.row = np.empty(self.n, np.int) # 各行の基底変数の番号 self.s = np.empty((self.n+1, self.mm+1), np.float) # 単体表 for i1 in range(0, self.n+1) : if i1 < self.n : self.s[i1][0] = eq_r[i1] for i2 in range(0, self.m) : self.s[i1][i2+1] = eq_l[i1][i2] for i2 in range(self.m+1, self.mm+1) : self.s[i1][i2] = 0.0 else : for i2 in range(0, self.mm+1) : self.s[i1][i2] = 0.0 # スラック変数 k = self.m + 1 for i1 in range(0, self.n) : if self.cp[i1] != 0 : if self.cp[i1] < 0 : self.s[i1][k] = 1.0 self.row[i1] = k - 1 else : self.s[i1][k] = -1.0 k += 1 # 人為変数 for i1 in range(0, self.n) : if self.cp[i1] >= 0 : self.s[i1][k] = 1.0 self.row[i1] = k - 1 self.a[i1] = 1 k += 1 # 目的関数 if self.m_a == 0 : for i1 in range(0, self.m) : self.s[self.n][i1+1] = -self.z[i1] else : for i1 in range(0 , self.m+self.m_s+1) : for i2 in range(0 , self.n) : if self.a[i2] > 0 : self.s[self.n][i1] -= self.s[i2][i1] ################################ # 最適化 # sw : =0 : 中間結果無し # =1 : 中間結果あり # return : =0 : 正常終了 # : =1 : 解無し ################################ def optimize(self, sw) : # フェーズ1 if sw > 0 : if self.m_a > 0 : print("\nphase 1") else : print("\nphase 2") self.opt_run(sw) # フェーズ2 if self.err == 0 and self.m_a > 0 : # 目的関数の変更 self.mm -= self.m_a for i1 in range(0, self.mm+1) : self.s[self.n][i1] = 0.0 for i1 in range(0 , self.n) : k = self.row[i1] if k < self.m : self.s[self.n][0] += self.z[k] * self.s[i1][0] for i1 in range(0, self.mm) : for i2 in range(0, self.n) : k = self.row[i2] if k < self.m : self.s[self.n][i1+1] += self.z[k] * self.s[i2][i1+1] if i1 < self.m : self.s[self.n][i1+1] -= self.z[i1] # 最適化 if sw > 0 : print("\nphase 2") self.opt_run(sw) return self.err ################################ # 最適化(単体表の変形) # sw : =0 : 中間結果無し # =1 : 中間結果あり ################################ def opt_run(self, sw) : self.err = -1 while self.err < 0 : # 中間結果 if sw > 0 : print() self.output() # 列の選択(巡回を防ぐため必ずしも最小値を選択しない,Bland の規則) q = -1 for i1 in range(1, self.mm+1) : if self.s[self.n][i1] < -self.eps : q = i1 - 1 break # 終了(最適解) if q < 0 : self.err = 0 # 行の選択( Bland の規則を採用) else : p = -1 k = -1 min = 0.0 for i1 in range(0, self.n) : if self.s[i1][q+1] > self.eps : x = self.s[i1][0] / self.s[i1][q+1] if (p < 0) or (x < min) or ((x == min) and (self.row[i1] < k)) : min = x p = i1 k = self.row[i1] # 解無し if p < 0 : self.err = 1 # 変形 else : x = self.s[p][q+1] self.row[p] = q for i1 in range(0, self.mm+1) : self.s[p][i1] /= x for i1 in range(0, self.n+1) : if i1 != p : x = self.s[i1][q+1] for i2 in range(0, self.mm+1) : self.s[i1][i2] -= (x * self.s[p][i2]) ################################ # 単体表の出力 ################################ def output(self) : for i1 in range(0, self.n+1) : if i1 < self.n : print("x" + str(self.row[i1]+1), end="") else : print(" z", end="") for i2 in range(0, self.mm+1) : print(" " + str(self.s[i1][i2]), end="") print() ################################ # 結果の出力 ################################ def result(self) : if self.err > 0 : print("\n解が存在しません") else : print("\n(", end="") for i1 in range(0, self.m) : x = 0.0 for i2 in range(0, self.n) : if self.row[i2] == i1 : x = self.s[i2][0] break if i1 == 0 : print(x, end="") else : print(", " + str(x), end="") print(") のとき,最大値 " + str(self.s[self.n][0])) ---------------------------------- # -*- coding: UTF-8 -*- import numpy as np import sys from math import * from function import LP ############################################ # 線形計画法 # coded by Y.Suganuma ############################################ # 入力 # 変数の数と式の数 line = sys.stdin.readline() x = line.split() m = int(x[0]) n = int(x[1]) cp = np.empty(n, np.int) z = np.empty(m, np.float) eq_r = np.empty(n, np.float) eq_l = np.empty((n, m), np.float) # 目的関数の係数 line = sys.stdin.readline() x = line.split() for i1 in range(0, m) : z[i1] = float(x[i1]) # 制約条件 for i1 in range(0, n) : line = sys.stdin.readline() x = line.split() for i2 in range(0, m) : eq_l[i1][i2] = float(x[i2]) if x[m] == '<' : cp[i1] = -1 elif x[m] == '>' : cp[i1] = 1 else : cp[i1] = 0 eq_r[i1] = float(x[m+1]) # 実行 lp = LP(n, m, z, eq_l, eq_r, cp) lp.optimize(1) # 結果の出力 lp.result()