Assignment 4
杨寿国 201718018670071
December 21, 2017
1 Gas Station Placement
假设 n 个加油站的位置分别为 x1:x2; x3; :::; xn。令 M = max
表⽰连续两个加油站之间的最⼤距离则有 LP:
i
8>>>>>><>>>>>>:
化为标准形式为 8>>><>>>:
min(M )
xi =< di + r; i = 1; 2; 3; :::; n
xi >= di r; i = 1; 2; 3; :::; n
xi xi1 M; i = 1; 2; 3; :::; n
xi 0; M > 0
min(M )
xi =< di + r; i = 1; 2; 3; :::; n
xi r di; i = 1; 2; 3; :::; n
xi xi1 M; i = 1; 2; 3; :::; n 0
(xi xi1) ,
(1)
(2)
2 Duality
L(; ; ; fi(u; v)) = 0
fi(u; v))
∑
u;(u;v)2E
fi(u; v))
(3)
(u;v)(c(u; v) k∑
∑
i=1
fi(u; v)
∑
i(
i=1
w;(v;w)2E
i(di
(u;v)2E
∑
k∑
k∑
∑
L = (
i=1
v;(s;v)2E
fi(si; v))
∑
1
g(; ; ) = inf
fi(u;v)
(u;v)2E (u;v)c(u; v) +
k
i=1 idi)
∑
∑
条件是 (u;v) 0; i 0;
u;(u;v)2E
其中 (u;v) 表⽰在节点 u 到节点 v 之间运输的代价
使⽤ lec8 课件中的数据构造原问题的⼀个实例:即 c(s; u) = 7; c(u; v) =
w;(v;w)2E
fi(u; v) =
fi(u; v)
(4)
1; c(s; v) = 8; c(u; t) = 6; c(v; t) = 2; d1 = 5; d2 = 3
得到⼀个线性规划问题为8>>>>>>>>>>>>>>>>>>>>>>><>>>>>>>>>>>>>>>>>>>>>>>:
max 0
f1(s; u) + f2(s; u) 7
f1(s; v) + f2(s; v) 8
f1(u; v) + f2(u; v) 1
f1(u; t) + f2(u; t) 6
f1(v; t) + f2(v; t) 2
f1(s; u) = f1(u; v) + f1(u; t)
f1(s; v) + f1(u; v) = f1(v; t)
f2(s; u) = f2(u; v) + f2(u; t)
f2(s; v) + f2(u; v) = f2(v; t)
f1(s; v) + f1(s; u) = 5
f2(s; v) + f2(s; u) = 3
使⽤ GLPK 进⾏求解,得到如下结果:
3 Dual Simplex Algorithm
⾃⼰代码运⾏结果截图如下
使⽤ glpk ⼯具进⾏求解得到:
2
可以看到结果是⼀样的,最⼩值是-16.5, 最优解是 (0,0,0,1.5,2.5,16.5)
对偶单纯形代码如下:
#encoding=utf 8
__author__ = ’ ysg ’
import numpy as np #python 矩阵操作 l i b
class DualSimplex ( ) :
def __init__( s e l f ) :
s e l f ._A = ”” # 系数矩阵
s e l f ._b = ”” #
s e l f . _c = ’ ’ #约束
s e l f ._B = ’ ’ #基 变 量 的 下 标 集 合
s e l f . row = 0 #约束个数
def
filename ) :
s o l v e ( s e l f ,
#读取⽂件内容 , ⽂ 件 结 构 前 两 ⾏ 分 别 为 变量数 和 约束条件个数
#接 下 来 是 系 数 矩 阵
#然后是 b 数组
#然 后 是 约 束 条 件 c
#假 设 线 性 规 划 形 式 是 标 准 形 式 ( 都是等式 )
A = [ ]
b = [ ]
c = [ ]
with open( filename , ’ r ’ ) as
f :
s e l f . var = int ( f . r e a d l i n e ( ) )
s e l f . row = int ( f . r e a d l i n e ( ) )
for i
in range ( s e l f . row ) :
x =map( int ,
A. append ( x)
f . r e a d l i n e ( ) . s t r i p ( ) . s p l i t ( ’ ␣ ’ ))
b=(map( int ,
l i s t ( f . r e a d l i n e ( ) . s p l i t ( ’ ␣ ’ ) ) ) )
3
c=(map( int ,
l i s t ( f . r e a d l i n e ( ) . s p l i t ( ’ ␣ ’ ) ) ) )
s e l f ._A = np . array (A, dtype=float )
s e l f ._b = np . array (b , dtype=float )
s e l f . _c = np . array ( c , dtype=float )
# s e l f ._A = np . array ([[3 , 1 ,1 , 2 ,0 ,0] ,[2 ,1 ,0 ,1 ,1 ,0] ,[ 1 ,3 ,0 , 3 ,0 ,1]] , dtype=f l o a t )
# s e l f . _b = np . array ([ 3 ,4 ,12] , dtype=f l o a t )
# s e l f . _c = np . array ([7 , 7 , 2, 1, 6, 0] , dtype=f l o a t )
s e l f ._B = l i s t ()
s e l f . row = len ( s e l f ._b)
s e l f . var = len ( s e l f . _c)
(x , obj ) = s e l f . DualSimplex ( s e l f ._A, s e l f ._b, s e l f . _c)
s e l f . pprint (x , obj ,A)
def pprint ( s e l f , x , obj ,A) :
px = [ ’x_%d␣=␣%f ’%( i +1,x [ i ] )
print
’ , ’ . j o i n ( px )
print ( ’ o b j e c t i v e ␣ value ␣ i s ␣ : ␣%f ’% obj )
’’
print
for i
in range ( len (A) ) :
for i
in range ( len (x ) ) ]
print
’%dth␣ l i n e ␣ c o n s t r a i n t ␣ value ␣ i s ␣ : ␣%f ’ % ( i +1, x . dot (A[ i ] ) )
#得 到 ⼀ 个 对 偶 可 ⾏ 初 始 解
def
I n i t i a l i z e S i m p l e x ( s e l f ,A, b , c , obj ) :
#初始基变量 现 在 加 什 么 数 ⽆ 所 谓 , 只要⼤⼩为m⾏就 ok
B = l i s t (np . arange ( c . shape [ 0 ] , c . shape [0]+b . shape [ 0 ] ) )
#变换基 , 得到⼀个基解
l = A. shape [0] 1
while l >= 0:
for j
in range ( s e l f . var 1,1,1):
i f A[ l , j ]
break ;
!= 0:
(B, A, b , c , obj ) = s e l f ._PIVOT(B, A, b , c , obj ,
l =1
l ,
j )
#如 果 不 满 ⾜ 对 偶 可 ⾏ , 添加约束条件
i f np .min( c ) < 0:
newline = np . ones (A. shape [ 1 ] )
newline [B] = 0
A = np . concatenate ((A, newline . reshape (1 , 1)) , axis =0)
4
newcolumn = np . zeros (A. shape [ 0 ] )
newcolumn [ 1] = 1
A = np . concatenate ((A, newcolumn . reshape ( 1 ,1)) , axis =1)
c = np . concatenate (( c , np . array ( [ 0 ] ) ) , axis =1)
b = np . concatenate (( b , np . array ( [ 1 0 0 0 0 0 0 ] ) ) , axis =1)
B. append (A. shape [1] 1)
l = len (B)1
e = np . argmin ( c )
(B, A, b , c , obj ) = s e l f ._PIVOT(B, A, b , c , obj ,
l , e )
return (B,A, b , c , obj )
#算法⼊⼜
def DualSimplex ( s e l f ,A, b , c ) :
# 函数⽬标值
obj = 0
(B, A ,b , c , obj ) = s e l f . I n i t i a l i z e S i m p l e x (A, b , c , obj )
#寻找 b i < 0
b_min = np .min(b)
; b_min_i = np . argmin (b)
while b_min < 0 :
#寻找最⼩ t h e t a
theta = [ ]
for j
in range (A. shape [ 1 ] ) :
i f A[ b_min_i , j ] < 0 :
theta . append ( c [ j ]/A[ b_min_i , j ]*1)
else :
theta . append (10000)
theta_j = np . argmin ( theta )
i f
theta [ theta_j ] == 10000:
print
break
’ no␣ f e a s i b l e ␣ s o l u t i o n ’
(B, A, b , c , obj ) = s e l f ._PIVOT(B,A, b , c , obj , b_min_i , theta_j )
# 寻找 b i < 0
b_min = np .min(b ) ;
b_min_i = np . argmin (b)
5
x = s e l f . _CalculateX (B,A, b , c )
return (x, obj )
#得到完整解
def _CalculateX ( s e l f ,B,A, b , c ) :
x = np . zeros ( s e l f . var , dtype=float )
x [B] = b
return x
# 基变换
def _PIVOT( s e l f ,B,A, b , c , z , l , e ) :
i s a_le
# main element
# l
l e a v i n g b a s i s
repre sents
# e repre sents entering b a s i s
main_elem = A[ l ] [ e ]
#s c a l i n g the lth l i n e
A[ l ] = A[ l ]/ main_elem
b [ l ] = b [ l ]/ main_elem
#change eth column to unit array
for i
i f
in range (A. shape [ 0 ] ) :
i
!= l :
b [ i ] = b [ i ] A[ i ] [ e ] * b [ l ]
A[ i ] = A[ i ] A[ i ] [ e ] * A[ l ]
#update o b j e c t i v e value
z = b [ l ]* c [ e ]
#c = c c [ e ] * A[ l ]
c = c c [ e ] * A[ l ]
# change the b a s i s
B[ l ] = e
return (B, A, b , c , z )
s = DualSimplex ()
s . s o l v e ( ’ pro2 . txt ’ )
pro2.txt ⽂件内容如下:
6
3
3 -1 1 -2 0 0
2 1 0 1 1 0
-1 3 0 -3 0 1
-3 4 12
-7 7 -2 -1 -6 0
6