logo资料库

卜东波算法第四次作业线性规划.pdf

第1页 / 共6页
第2页 / 共6页
第3页 / 共6页
第4页 / 共6页
第5页 / 共6页
第6页 / 共6页
资料共6页,全文预览结束
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
分享到:
收藏