prob5.2_3stages.jl

using JuMP, Clp

n = 4
m = 3
# Investment cost
ic = [16, 5, 32, 2]
# Fuel cost
C = [25, 80, 6.5, 160]
# Duration
T = [8760, 7000, 1500] / 8760
# Height
D2 = [diff([0, 3919, 7329, 10315])  diff([0, 7086, 9004, 11169])]


mod = Model(solver=ClpSolver())
@variables(mod, begin
    x1[1:n] >= 0
    x2[1:n] >= 0
    y2[1:n,1:m, 1:2]  >= 0
    y3[1:n,1:m, 1:2, 1:2] >= 0
end)
@constraints(mod, begin
    [i=1:n, s=1:2], sum(y2[i,:, s]) <= x1[i]
    [i=1:n, s1=1:2, s2=1:2], sum(y3[i,:, s1,s2]) <= x2[i]
    [i=1:n], x2[i] >= x1[i]
    [j=1:m, s=1:2], sum(y2[:,j, s]) == D2[j,s]
    [j=1:m, s1=1:2,s2=1:2], sum(y3[:,j, s1, s2])  == D2[j,s2]
end)
obj = dot(ic, x2)
for s1 in 1:2
    append!(obj, 0.5 * dot(C, y2[:,:,s1] * T))
    for s2 in 1:2
        append!(obj, 0.25 * dot(C, y3[:,:,s1,s2] * T))
    end
end
@objective(mod, Min, obj)

solve(mod)

det_obj = getobjectivevalue(mod)
det_x   = getvalue(x1)

println("""
    Deterministic Equivalent Solved.

                  Objective = $(det_obj)
     First stage investment = $(det_x)
    Second stage investment = $(getvalue(x2))

""")

#==
    Deterministic Equivalent Solved.

                  Objective = 437367.7157534246
     First stage investment = [1311.0, 854.0, 9004.0, 0.0]
    Second stage investment = [1311.0, 854.0, 9004.0, 0.0]
==#