妖怪幻 发表于 2017-5-5 07:27:46

Simplex单纯性算法的Python实现

  单纯性算法是解决线性规划的经典方法,上世纪50年代就提出了,其基本思想是在可行域内沿着边遍历所有的顶点,找出最优值,即为算法的最优值。
  算法的执行过程如下:
  


[*]求出初始基向量
[*]构建单纯性表格
[*]在所有非基向量对应的c中,找出一个最小的<nobr style=""><span class="math" id="MathJax-Span-71" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-72" style=""><span class="msubsup" id="MathJax-Span-73" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-74" style="">c</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-75" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>,若该<nobr style=""><span class="math" id="MathJax-Span-76" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-77" style=""><span class="msubsup" id="MathJax-Span-78" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-79" style="">c</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-80" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>大于0,则退出,输出obj,否则将<nobr style=""><span class="math" id="MathJax-Span-81" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-82" style=""><span class="msubsup" id="MathJax-Span-83" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-84" style="">a</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-85" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>入基
[*]利用基向量组线性表示<nobr style=""><span class="math" id="MathJax-Span-86" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-87" style=""><span class="msubsup" id="MathJax-Span-88" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-89" style="">a</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-90" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>,得到该线性表示的系数向量<nobr style=""><span class="math" id="MathJax-Span-91" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-92" style=""><span class="mi" id="MathJax-Span-93" style="">Λ</span></span><span style=""></span></span></span><span style=""></span></span></nobr>
[*]对于<nobr style=""><span class="math" id="MathJax-Span-94" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-95" style=""><span class="mi" id="MathJax-Span-96" style="">Λ</span></span><span style=""></span></span></span><span style=""></span></span></nobr>中所有大于0的分量,求出<nobr style=""><span class="math" id="MathJax-Span-97" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-98" style=""><span class="munderover" id="MathJax-Span-99" style=""><span style=""><span style=""><span class="mo" id="MathJax-Span-100" style="">min</span><span style=""></span></span><span style=""><span class="texatom" id="MathJax-Span-101" style=""><span class="mrow" id="MathJax-Span-102" style=""><span class="mi" id="MathJax-Span-103" style="">m</span></span></span><span style=""></span></span><span style=""><span class="texatom" id="MathJax-Span-104" style=""><span class="mrow" id="MathJax-Span-105" style=""><span class="mi" id="MathJax-Span-106" style="">j<span style=""></span></span><span class="mo" id="MathJax-Span-107" style="">=</span><span class="mn" id="MathJax-Span-108" style="">1</span></span></span><span style=""></span></span></span></span><span class="mfrac" id="MathJax-Span-109" style=""><span style=""><span style=""><span class="msubsup" id="MathJax-Span-110" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-111" style="">b</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-112" style="">j<span style=""></span></span><span style=""></span></span></span></span><span style=""></span></span><span style=""><span class="msubsup" id="MathJax-Span-113" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-114" style="">Λ</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-115" style="">j<span style=""></span></span><span style=""></span></span></span></span><span style=""></span></span><span style=""><span style=""></span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>对应的j,然后将<nobr style=""><span class="math" id="MathJax-Span-116" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-117" style=""><span class="msubsup" id="MathJax-Span-118" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-119" style="">a</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-120" style="">j<span style=""></span></span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>出基
[*]问题矩阵旋转,即通过高斯行变换,将<nobr style=""><span class="math" id="MathJax-Span-121" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-122" style=""><span class="msubsup" id="MathJax-Span-123" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-124" style="">a</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-125" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>变换成<nobr style=""><span class="math" id="MathJax-Span-126" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-127" style=""><span class="msubsup" id="MathJax-Span-128" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-129" style="">a</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-130" style="">j<span style=""></span></span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>
[*]如果<nobr style=""><span class="math" id="MathJax-Span-131" style=""><span style=""><span style=""><span class="mrow" id="MathJax-Span-132" style=""><span class="msubsup" id="MathJax-Span-133" style=""><span style=""><span style=""><span class="mi" id="MathJax-Span-134" style="">c</span><span style=""></span></span><span style=""><span class="mi" id="MathJax-Span-135" style="">i</span><span style=""></span></span></span></span></span><span style=""></span></span></span><span style=""></span></span></nobr>全大于0,则退出算法,输出obj,否则,重复第3步

__author__ = 'linfuyuan'

class Problem:
def __init__(self):
self.obj = 0
self.coMatrix = []
self.b = []
self.c = []
def pivot(self, basic, l, e):
# the l-th line
self.b /= self.coMatrix
origin = self.coMatrix
for i in range(len(self.coMatrix)):
self.coMatrix /= origin
# the other lines
for i in range(len(self.b)):
if i != l:
self.b = self.b - self.b * self.coMatrix
for i in range(len(self.b)):
if i != l:
origin = self.coMatrix
for j in range(len(self.coMatrix)):
self.coMatrix = self.coMatrix - origin * self.coMatrix
origin = self.c
for i in range(len(self.coMatrix)):
self.c = self.c - self.coMatrix * origin
self.obj = self.obj - self.b * origin
basic = e
def clone(self, another):
self.obj = another.obj
for i in another.b:
self.b.append(i)
for i in another.c:
self.c.append(i)
for v in another.coMatrix:
newv = []
for i in v:
newv.append(i)
self.coMatrix.append(newv)

basic = []
problem = Problem()

def readProblem(filename):
with open(filename)as f:
var_num = int(f.readline())
constrait_num = int(f.readline())
matrix = [( * var_num) for i in range(constrait_num)]
for i in range(constrait_num):
line = f.readline()
tokens = line.split(' ')
for j in range(var_num):
matrix = float(tokens)
problem.b.append(float(tokens[-1]))
for i in range(var_num):
var = []
for j in range(constrait_num):
var.append(matrix)
problem.coMatrix.append(var)
line = f.readline()
tokens = line.split(' ')
for i in range(var_num):
problem.c.append(float(tokens))

def getMinCIndex(c):
min, minIndex = 1, 0
for i in range(len(c)):
if c < 0 and c < min:
min = c
minIndex = i
if min > 0:
return -1
else:
return minIndex

def getLamdaVector(evector):
ld = []
for i in range(len(evector)):
ld.append(evector)
return ld

def simplex(basic, problem):
minCIndex = getMinCIndex(problem.c)
while minCIndex != -1:
ld = getLamdaVector(problem.coMatrix)
# find the l line
l, min = -1, 10000000000
for i in range(len(problem.b)):
if ld > 0:
if problem.b / ld < min:
l = i
min = problem.b / ld
if l == -1:
return False
problem.pivot(basic, l, minCIndex)
minCIndex = getMinCIndex(problem.c)
return True

def initialSimplex(basic, problem):
min, minIndex = 1000000000, -1
for i in range(len(problem.b)):
if problem.b < min:
min = problem.b
minIndex = i
for i in range(len(problem.b)):
basic.append(i+len(problem.b))
if min >= 0:
return True
else:
originC = problem.c
newC = []
for i in range(len(problem.c)):
newC.append(0)
newC.append(1)
problem.c = newC
x0 = []
for i in range(len(problem.b)):
x0.append(-1)
problem.coMatrix.append(x0)
problem.pivot(basic, minIndex, -1)
res = simplex(basic, problem)
if res == False or problem.obj != 0:
return False
else:
problem.c = originC
problem.coMatrix.pop()
# Gaussian row operation
counter = 0
for i in basic:
if problem.c != 0:
origin = problem.c
for j in range(len(problem.c)):
problem.c -= problem.coMatrix * origin
problem.obj -= problem.b * origin
counter += 1
return True

filename = raw_input('please input the problem description file: ')
readProblem(filename)
if initialSimplex(basic, problem):
res = simplex(basic, problem)
if res:
print 'the optimal obj is %.2f' % problem.obj
index = ['0.00'] * len(problem.coMatrix)
counter = 0
for i in basic:
index = '%.2f' % problem.b
counter += 1
print 'the solution is {%s}' % ' '.join(index)
else:
print 'no feasible solution'
else:
print 'no feasible solution'
raw_input('please input any key to quit.')
希望对大家有所帮助。
页: [1]
查看完整版本: Simplex单纯性算法的Python实现