赞 | 7 |
VIP | 0 |
好人卡 | 1 |
积分 | 18 |
经验 | 34644 |
最后登录 | 2024-6-30 |
在线时间 | 951 小时 |
Lv3.寻梦者
- 梦石
- 0
- 星屑
- 1784
- 在线时间
- 951 小时
- 注册时间
- 2012-7-5
- 帖子
- 245
|
加入我们,或者,欢迎回来。
您需要 登录 才可以下载或查看,没有帐号?注册会员
x
线性规划能做很多事啊。
听说叶子姐姐会感兴趣然后就拿出了这个破轮子...
具体用法写在注释里面啦....
Inf这次我可用Float::INFINITY表示啦
代码写得丑,大概凑合看吧....
自觉回沟
simplex.zip
(1.49 KB, 下载次数: 83)
=begin 线性规划的单纯形法 较为详细的资料见 [url]https://github.com/sxysxy/math/tree/master/simplex[/url] 用法: lp = LinearProgramming.new(n, [a1, a2, .. an]) 创建一个线性规划对象,它具有n个变量,目标函数为z = a1*x1 + a2*x2 + ... + an*xn lp.addAssert([a1, a2, ..an], b) 添加一条约束,为 a1*x1 + a2*x2 + ... + an*xn <= b 默认具有约束 x1 >= 0, x2 >= 0, ..., xn >= 0 (非负约束) 对于不满足上述形式的约束请使用适当的方法转变为上述形式 lp.solve 求解 如果无解返回 :no_solution 如果目标函数发散返回 :infinity 否则返回 [z, [x1, x2, x3, ..., xn]],目标函数的最大值和此时变量的取值 例子: lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3 lp.addAssert [1, 1, 3], 30 #x1 + x2 + 3x3 <= 30 lp.addAssert [2, 2, 5], 24 #2x1 + 2x2 + 5x3 <= 24 lp.addAssert [4, 1, 2], 36 #4x1 + x2 + 2x3 <= 36 print lp.solve #[28.0, [8.0, 4.0, 0.0]] 目标函数最大值为28,此时x1 = 8, x2 = 4, x3 = 0 =end class LinearProgramming NO_SOLUTION = :no_solution INFINITY = :infinity attr_reader :xs #变量数 attr_reader :ys #约束条件数 def initialize(xs, f) #f为目标函数 @xs = xs @ys = 0 @a = [] @a << [0.0]+f.map {|e| Float(e)} @base = [] @unbase = [] end def addAssert(x, b) @a << ([b]+x).map{|e| Float(e)} @ys += 1 end def solve r = simplex return NO_SOLUTION if r == 0 return INFINITY if r == -1 @unbase.fill 0 (1..@ys).each {|i| @unbase[@base[i]] = i if @base[i] <= @xs} return [@a[0][0], (1..@xs).map{|i| @unbase[i]!=0? @a[@unbase[i]][0]: 0.0}] end private EPS = 1e-9 def dcmp(x) return -1 if x<-EPS return 1 if x>EPS return 0 end def swapVar(x, y) t = @base[x] @base[x] = @unbase[y] @unbase[y] = t end def pivot(x, y) swapVar(x, y) k = @a[x][y] @a[x][y] = 1.0 @a[x].map! {|e| e/k} (0..@ys).each do |i| next if x == i || dcmp(@a[i][y]) == 0 k = @a[i][y] @a[i][0] += (i==0? 1:-1)*k*@a[x][0] @a[i][y] = 0 (1..@xs).each do |j| @a[i][j] -= k*@a[x][j] end end end def init (1..@xs).each{|e| @unbase[e] = e} (1..@ys).each{|e| @base[e] = e+@xs} x = y = 0 while true (1..@ys).each do |i| x = i if dcmp(@a[i][0]) < 0 end return true if x == 0 (1..@xs).each do |i| y = i if dcmp(@a[x][i]) < 0 end return false if y == 0 pivot(x, y) x = y = 0 end end def simplex return 0 unless init x = y = 0 while true (1..@xs).each do |i| if dcmp(@a[0][i]) > 0 y = i break end end return 1 if y == 0 inf = Float::INFINITY (1..@ys).each do |i| t = @a[i][0]/@a[i][y] if dcmp(@a[i][y])>0 && (x==0 || t<inf) inf = t x = i end end return -1 if x == 0 pivot(x, y) x = y = 0 end end end lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3 lp.addAssert [1, 1, 3], 30 #x1 + x2 + 3x3 <= 30 lp.addAssert [2, 2, 5], 24 #2x1 + 2x2 + 5x3 <= 24 lp.addAssert [4, 1, 2], 36 #4x1 + x2 + 2x3 <= 36 print lp.solve
=begin
线性规划的单纯形法
较为详细的资料见 [url]https://github.com/sxysxy/math/tree/master/simplex[/url]
用法:
lp = LinearProgramming.new(n, [a1, a2, .. an])
创建一个线性规划对象,它具有n个变量,目标函数为z = a1*x1 + a2*x2 + ... + an*xn
lp.addAssert([a1, a2, ..an], b)
添加一条约束,为 a1*x1 + a2*x2 + ... + an*xn <= b
默认具有约束 x1 >= 0, x2 >= 0, ..., xn >= 0 (非负约束)
对于不满足上述形式的约束请使用适当的方法转变为上述形式
lp.solve 求解
如果无解返回 :no_solution
如果目标函数发散返回 :infinity
否则返回 [z, [x1, x2, x3, ..., xn]],目标函数的最大值和此时变量的取值
例子:
lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
lp.addAssert [1, 1, 3], 30 #x1 + x2 + 3x3 <= 30
lp.addAssert [2, 2, 5], 24 #2x1 + 2x2 + 5x3 <= 24
lp.addAssert [4, 1, 2], 36 #4x1 + x2 + 2x3 <= 36
print lp.solve #[28.0, [8.0, 4.0, 0.0]] 目标函数最大值为28,此时x1 = 8, x2 = 4, x3 = 0
=end
class LinearProgramming
NO_SOLUTION = :no_solution
INFINITY = :infinity
attr_reader :xs #变量数
attr_reader :ys #约束条件数
def initialize(xs, f) #f为目标函数
@xs = xs
@ys = 0
@a = []
@a << [0.0]+f.map {|e| Float(e)}
@base = []
@unbase = []
end
def addAssert(x, b)
@a << ([b]+x).map{|e| Float(e)}
@ys += 1
end
def solve
r = simplex
return NO_SOLUTION if r == 0
return INFINITY if r == -1
@unbase.fill 0
(1..@ys).each {|i| @unbase[@base[i]] = i if @base[i] <= @xs}
return [@a[0][0], (1..@xs).map{|i| @unbase[i]!=0? @a[@unbase[i]][0]: 0.0}]
end
private
EPS = 1e-9
def dcmp(x)
return -1 if x<-EPS
return 1 if x>EPS
return 0
end
def swapVar(x, y)
t = @base[x]
@base[x] = @unbase[y]
@unbase[y] = t
end
def pivot(x, y)
swapVar(x, y)
k = @a[x][y]
@a[x][y] = 1.0
@a[x].map! {|e| e/k}
(0..@ys).each do |i|
next if x == i || dcmp(@a[i][y]) == 0
k = @a[i][y]
@a[i][0] += (i==0? 1:-1)*k*@a[x][0]
@a[i][y] = 0
(1..@xs).each do |j|
@a[i][j] -= k*@a[x][j]
end
end
end
def init
(1..@xs).each{|e| @unbase[e] = e}
(1..@ys).each{|e| @base[e] = e+@xs}
x = y = 0
while true
(1..@ys).each do |i|
x = i if dcmp(@a[i][0]) < 0
end
return true if x == 0
(1..@xs).each do |i|
y = i if dcmp(@a[x][i]) < 0
end
return false if y == 0
pivot(x, y)
x = y = 0
end
end
def simplex
return 0 unless init
x = y = 0
while true
(1..@xs).each do |i|
if dcmp(@a[0][i]) > 0
y = i
break
end
end
return 1 if y == 0
inf = Float::INFINITY
(1..@ys).each do |i|
t = @a[i][0]/@a[i][y]
if dcmp(@a[i][y])>0 && (x==0 || t<inf)
inf = t
x = i
end
end
return -1 if x == 0
pivot(x, y)
x = y = 0
end
end
end
lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
lp.addAssert [1, 1, 3], 30 #x1 + x2 + 3x3 <= 30
lp.addAssert [2, 2, 5], 24 #2x1 + 2x2 + 5x3 <= 24
lp.addAssert [4, 1, 2], 36 #4x1 + x2 + 2x3 <= 36
print lp.solve
|
评分
-
查看全部评分
|