拿出破轮子:线性规划的单纯形法

浮云半仙 · 发表于 2016-12-15 20:42:11

加入我们，或者，欢迎回来。

您需要登录才可以下载或查看，没有帐号？注册会员

x

线性规划能做很多事啊。
听说叶子姐姐会感兴趣然后就拿出了这个破轮子...
具体用法写在注释里面啦....
Inf这次我可用Float::INFINITY表示啦
代码写得丑，大概凑合看吧....
自觉回沟

simplex.zip (1.49 KB, 下载次数: 83)

RUBY 代码复制

=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

shitake · 发表于 2016-12-15 21:25:16

嘛代码写好是慢慢练习的事。ACM的话，因为通常是写“一次性”代码【要么过，要么不过换下一种】。所以对很多地方要求并不高【也没那没多时间去在意】。
但是实际开发恰恰相反，强调代码的可维护性。
ruby的话，可以用用RubyCritic或者rubocop做代码检查。坚持做一段时间的话，基本就没问题啦
还可以参考参考这个：https://github.com/bbatsov/ruby-style-guide

晴兰 · 发表于 2016-12-16 08:38:26

提示: 作者被禁止或删除内容自动屏蔽

帐号		自动登录	找回密码
密码			注册会员

[原创发布] 拿出破轮子:线性规划的单纯形法

加入我们，或者，欢迎回来。

点评

评分

点评

开拓者