# -*- 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()