In [ ]:

    

using JuMP
using Gadfly


    

In [ ]:

    

eps = 1e-6


    

In [ ]:

    

function solve_spp(num_systems, num_decks, work_rate_budget, difficulty_rate)
    m = Model()
    @defVar(m, u[1:num_systems, 1:num_decks] >= 0)
    @defVar(m, a >= 0)
    @defVar(m, l[1:num_systems, 1:num_decks] >= 0)
    @addConstraint(m, a + sum(u) <= work_rate_budget)
    for i = 1:num_systems        
        l_max = -log(1-i/num_systems+eps)/difficulty_rate
        l_min = -log(1-(i-1)/num_systems+eps)/difficulty_rate
        difficulty = (exp(-l_min * difficulty_rate) * (l_min * difficulty_rate + 1) - 
            exp(-l_max * difficulty_rate) * (l_max * difficulty_rate + 1)) / difficulty_rate
        for j = 1:num_decks
            @addConstraint(m, l[i, j] <= u[i, j])
        end
        @addNLConstraint(m, l[i, 1] == a/num_systems + (1 - (u[i, 1] - l[i, 1]) / 
            (u[i, 1] - l[i, 1] + difficulty)) * l[i, 1] + (1 - (u[i, 2] - l[i, 2]) / 
            (u[i, 2] - l[i, 2] + difficulty / 2)) * l[i, 2])
        for j = 2:num_decks-1
            @addNLConstraint(m, l[i, j] == (u[i, j-1] - l[i, j-1]) / (u[i, j-1] - l[i, j-1] 
                + difficulty / (j-1)) * l[i, j-1] + (1 - (u[i, j+1] - l[i, j+1]) / 
                (u[i, j+1] - l[i, j+1] + difficulty / (j+1))) * l[i, j+1])
        end
        @addNLConstraint(m, l[i, num_decks] == (u[i, num_decks-1] - l[i, num_decks-1]) /
            (u[i, num_decks-1] - l[i, num_decks-1] + difficulty / (num_decks-1)) * l[i, num_decks-1])
    end
    @setNLObjective(m, Max, a)
    TT = STDOUT
    redirect_stdout()
    status = solve(m)
    redirect_stdout(TT)
    return getObjectiveValue(m), getValue(u), getValue(l)
end


    

In [ ]:

    

num_decks = 10
work_rate_budget = 0.19020740740740741


    

In [ ]:

    

num_systems = 1
difficulty_rate = 1 / 0.0076899999999998905


    

In [ ]:

    

solve_spp(num_systems, num_decks, work_rate_budget, difficulty_rate)


    

In [ ]:

    

difficulty_rates = [1 / 0.0076899999999998905]


    

In [ ]:

    

difficulty_rates = 1 ./ collect(0.001:0.001:0.01)


    

In [ ]:

    

num_systems_set = 1:1:5


    

In [ ]:

    

throughputs = [-1. for x in difficulty_rates, y in num_systems_set]
allocations = [[-1. for z in 1:y, w in 1:num_decks] for x in difficulty_rates, y in num_systems_set];


    

In [ ]:

    

for (i, difficulty_rate) in enumerate(difficulty_rates)
    for (j, num_systems) in enumerate(num_systems_set)
        m, u, l = solve_spp(num_systems, num_decks, work_rate_budget, difficulty_rate)
        throughputs[i, j] = m
        allocations[i, j] = u
    end
end


    

In [ ]:

    

colors = ["gray", "gray", "yellow", "yellow", "orange", "orange", "pink", "pink", "red", "red"]


    

In [ ]:

    

d = 1
p1 = plot(x=[z[1] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))], 
y=[z[2] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))], Geom.point, Geom.line,
Guide.xlabel("Number of Systems"), Guide.ylabel("Max Arrival Rate"), Guide.title("Value of Information"))


    

In [ ]:

    

p2 = plot(map(d -> layer(x=[z[1] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))], 
y=[z[2] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))], 
Theme(default_color=color(colors[d])), Geom.line)[1], 1:length(difficulty_rates)), 
Guide.xlabel("Number of Systems"), Guide.ylabel("Max Arrival Rate"), Guide.title("Value of Information"))


    

In [ ]:

    

draw(PDF("figures/lqn/ext1.pdf", 12cm, 8cm), p2)


    

In [ ]:

    

p3 = plot(map(d -> layer(x=[z[1] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))], 
y=([z[2] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))] .- minimum([z[2] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))])) ./ minimum([z[2] for z in collect(filter(a -> a[2]>0, zip(1:num_systems, throughputs[d, :])))]), 
Theme(default_color=color(colors[d])), Geom.line)[1], 1:length(difficulty_rates)), 
Guide.xlabel("Number of Systems"), Guide.ylabel("Relative Throughput Gain"), Guide.title("Value of Information"))


    

In [ ]:

    

draw(PDF("figures/lqn/ext2.pdf", 12cm, 8cm), p3)


    

In [ ]:

    

n = 10
p4 = plot(map(d -> layer(x=1:n, y=[sum(allocations[d, n][x, :]) for x=1:n], 
Theme(default_color=color(colors[d])), Geom.line)[1], 1:length(difficulty_rates)), 
Guide.xlabel("System"), Guide.ylabel("Total Work Rate"), Guide.title("Optimal Work Rate Allocation"))


    

In [ ]:

    

draw(PDF("figures/lqn/ext3.pdf", 12cm, 8cm), p4)


    

In [ ]:

    

colors = ["gray", "yellow", "orange", "pink", "red"]


    

In [ ]:

    

n = 5
p5 = plot(map(d -> layer(x=1:num_decks, y=allocations[1, n][d, :], 
Theme(default_color=color(colors[d])), Geom.line)[1], 1:n), 
Guide.xlabel("Deck"), Guide.ylabel("Work Rate"), Guide.title("Optimal Work Rate Allocation"))


    

In [ ]:

    

draw(PDF("figures/lqn/ext4.pdf", 12cm, 8cm), p5)


    

In [ ]:

    

n = 5
p6 = plot(map(d -> layer(x=1:num_decks, y=(allocations[1, n][d, :] ./ sum(allocations[1, n][d, :])), 
Theme(default_color=color(colors[d])), Geom.line)[1], 1:n), 
Guide.xlabel("Deck"), Guide.ylabel("Work Rate (Normalized)"), Guide.title("Optimal Work Rate Allocation"))


    

In [ ]:

    

draw(PDF("figures/lqn/ext5.pdf", 12cm, 8cm), p6)


    

In [ ]:

    

println(collect(throughputs))


    

In [ ]:

    

for x=1:num_systems
    println(collect(allocations[1, x]))
end


    

In [ ]: