initial commit

.gitignore

@@ -1 +1,3 @@
 Manifest.toml
+.DS_Store
+.log

README.md

@@ -1 +1,2 @@
 # DRMSMIP
+Distributionally Robust Multistage Stochastic Mixed Integer Programming

examples/investment.jl

@@ -0,0 +1,209 @@
+using JuMP, Gurobi
+using DRMSMIP
+using DualDecomposition
+using Random
+
+const DD = DualDecomposition
+
+const rng = Random.MersenneTwister(1234)
+"""
+α: interest rate
+π: unit stock price
+ρ: unit dividend price
+
+
+K = 3 number of stages
+L = 2 number of investment vehicles
+2^L scenarios in each stage
+2^L^(K-1)=16 scenarios in total
+ρ = 0.05 * π
+bank: interest rate 0.01
+asset1: 1.03 or 0.97
+asset2: 1.06 or 0.94
+
+In each node, we have Np=10 samples from a log-normal distribution
+"""
+
+function create_tree(K::Int, L::Int, Np::Int)::DRMSMIP.Tree
+    π = ones(L)
+    π_samp = generate_sample(L, π, Np)
+    set = DRMSMIP.WassersteinSet(π_samp, 1)
+    cost = zeros(L+1)
+    tree = DRMSMIP.Tree(π, set, cost)
+    add_nodes!(K, L, tree, 1, 1, Np)
+    return tree
+end
+
+
+function generate_sample(L::Int, π::Array{Float64}, Np::Int)::Array{DRMSMIP.Sample}
+    # generates random samples following a lognormal distribution
+    ret = Array{DRMSMIP.Sample}(undef, Np)
+    for ii in 1:Np
+        ξ = Array{Float64}(undef, L)
+        for l in 1:L
+            sig = sqrt( log( 0.5+sqrt( ( 0.03*l )^2+0.25 ) ) )
+            rnd = sig * randn(rng) .+ log(π[l])
+            ξ[l] = exp(rnd)
+        end
+        ret[ii] = DRMSMIP.Sample(ξ, 1/Np)
+    end
+    return ret
+end
+
+function add_nodes!(K::Int, L::Int, tree::DRMSMIP.Tree, id::Int, k::Int, Np::Int)
+    if k < K-1
+        ls = iterlist(L,tree.nodes[id].ξ)
+        for π in ls
+            π_samp = generate_sample(L, π, Np)
+            set = DRMSMIP.WassersteinSet(π_samp, 1)
+            cost = zeros(L+1)
+            DRMSMIP.addchild!(tree, id, π, set, cost)
+            childid = length(tree.nodes)
+            add_nodes!(K, L, tree, childid, k+1, Np)
+        end
+    elseif k == K-1
+        ls = iterlist(L,tree.nodes[id].scenario)
+        for π in ls
+            cost = vcat(π, [1])
+            DRMSMIP.addchild!(tree, id, k+1, π, nothing, cost)
+        end
+    end
+end
+
+function iterlist(L::Int, π::Array{Float64})::Array{Float64, 2}
+    # generates all combinations of up and down scenarios
+    ls = Array{Float64, 2}(undef, 2^L, L)
+    ii = 1
+
+    function foo(L::Int, l::Int, arr::Vector{Float64})
+        up = (1.0 + 0.03 * l) * π[l]
+        dn = (1.0 - 0.03 * l) * π[l]
+
+        if l < L
+            arr1 = copy(arr)
+            arr1[l] = up
+            foo(L, l+1, arr1)
+
+            arr2 = copy(arr)
+            arr2[l] = dn
+            foo(L, l+1, arr2)
+        else
+            arr1 = copy(arr)
+            arr1[l] = up
+            ls[ii] = arr1
+            ii+=1
+
+            arr2 = copy(arr)
+            arr2[l] = dn
+            ls[ii] = arr2
+            ii+=1
+        end
+    end
+
+    foo(L, 1, Array{Float64}(undef, L))
+    return ls
+end
+
+"""
+    max     B_K+∑_{l=1}^{L}π_{K,l}y_{K,l}
+
+    s.t.    B_1+∑_{l=1}^{L}π_{1,l}x_{1,l} = b_1
+
+            b_k+(1+α)B_{k-1}+∑_{l=1}^{L}ρ_{k,l}y_{k-1,l} = B_k+∑_{l=1}^{L}π_{k,l}x_{k,l}, ∀ k=2,…,K
+    
+            y_{1,l} = x_{1,l}, ∀ l=1,…,L
+    
+            y_{k-1,l}+x_{k,l} = y_{k,l}, ∀ k=2,…,K, l=1,…,L
+    
+            x_{k,l} ∈ ℤ , ∀ k=1,…,K, l=1,…,L
+    
+            y_{k,l} ≥ 0, ∀ k=1,…,K, l=1,…,L
+    
+            B_k ≥ 0, ∀ k=1,…,K.
+"""
+const K = 3
+const L = 2
+const α = 0.01
+const b = [100, 30, 30]
+const Np = 10
+
+function create_scenario_model(K::Int, L::Int, tree::DRMSMIP.Tree, id::Int)
+    hist = DRMSMIP.get_history(tree, id)
+    m = Model(Gurobi.Optimizer) 
+    set_optimizer_attribute(m, "OutputFlag", 0)
+    @variable(m, x[1:K,1:L], integer=true)
+    @variable(m, y[1:K,1:L]>=0)
+    @variable(m, B[1:K]>=0)
+
+    π = tree.nodes[1].ξ
+
+    con = @constraint(m, B[1] + sum( π[l] * x[1,l] for l in 1:L) == b[1])
+    set_name(con, "con[1]")
+
+    for l in 1:L
+        bal = @constraint(m, y[1,l]-x[1,l]==0)
+        set_name(bal, "bal[1,$(l)]")
+    end
+
+    for k = 2:K
+        π = tree.nodes[hist[k]].ξ
+        ρ = tree.nodes[hist[k-1]].ξ * 0.05
+
+        con = @constraint(m, B[k] + sum( π[l] * x[k,l] - ρ[l] * y[k,l] for l in 1:L)
+            - (1+α) * B[k-1] == b[k])
+        set_name(con, "con[$(k)]")
+        for l in 1:L
+            bal = @constraint(m, y[k,l]-x[k,l]-y[k-1,l]==0)
+            set_name(bal, "bal[$(k),$(l)]")
+        end
+    end
+    π = tree.nodes[id].ξ
+    #@objective(m, Max, B[K] + sum( π[l] * x[K,l] for l in 1:L )
+    @objective(m, Max, 0 )
+    return m
+end
+
+
+tree = create_tree(K,L,Np)
+
+# Create DualDecomposition instance.
+algo = DD.LagrangeDual()
+
+# Add Lagrange dual problem for each scenario s.
+nodes = DRMSMIP.get_stage_id(tree)
+models = Dict{Int,JuMP.Model}(id => create_scenario_model(K,L,tree,id) for id in nodes[K])
+for id in nodes[K]
+    DD.add_block_model!(algo, id, models[id])
+end
+
+coupling_variables = Vector{DD.CouplingVariableRef}()
+for k in 1:K-1
+    for root in nodes[k]
+        leaves = DRMSMIP.get_future(tree, root)
+        for id in leaves
+            model = models[id]
+            xref = model[:x]
+            for l in 1:L
+                push!(coupling_variables, DD.CouplingVariableRef(id, [root, l], xref[k, l]))
+            end
+            Bref = model[:B]
+            push!(coupling_variables, DD.CouplingVariableRef(id, [root, L+1], Bref[k]))
+        end
+    end
+end
+# dummy coupling variables
+for id in nodes[K]
+    model = models[id]
+    xref = model[:x]
+    for l in 1:L
+        push!(coupling_variables, DD.CouplingVariableRef(id, [id, l], xref[K, l]))
+    end
+end
+
+# Set nonanticipativity variables as an array of symbols.
+DD.set_coupling_variables!(algo, coupling_variables)
+
+
+
+# Solve the problem with the solver; this solver is for the underlying bundle method.
+DD.run!(algo, optimizer_with_attributes(Gurobi.Optimizer, "print_level" => 0))

src/DRMSMIP.jl

@@ -1,5 +1,124 @@
 module DRMSMIP
 
-# Write your package code here.
+using JuMP, Gurobi
+using DualDecomposition, BundleMethod
 
+const DD = DualDecomposition
+const BM = BundleMethod
+
+
+"""
+    Ambiguity Set
+"""
+struct Sample
+    ξ::Vector{Float64}  # sampled scenario
+    p::Float64          # associated probability
+end
+
+abstract type AbstractAmbiguitySet end
+
+struct WassersteinSet <: AbstractAmbiguitySet
+    samples::Array{Sample}  # empirical distribution
+    N::Int                  # number of distinct samples
+    ϵ::Float64              # radius of Wasserstein Ball
+    WassersteinSet(samples::Array{Sample}, ϵ::Float64) = new(samples, length(samples), ϵ)
+end
+
+"""
+    Scenario Tree
+"""
+struct TreeNode
+    parent::Int                                 # index of parent node
+    children::Vector{Int}                       # indices of child nodes
+    k::Int                                      # current stage
+    ξ::Vector{Float64}                          # current scenario
+    set::Union{AbstractAmbiguitySet, Nothing}   # ambiguity set for child nodes
+    cost::Vector{Float64}
+end
+mutable struct Tree
+    nodes::Vector{TreeNode}     # list of nodes
+    K::Int                      # length of tree
+end
+
+Tree(ξ::Vector{Float64}, set::AbstractAmbiguitySet, cost::Vector{Float64}) = Tree([TreeNode(0, Vector{Int}(), 1, ξ, set, cost )], 1)
+
+function addchild!(tree::Tree, id::Int, ξ::Vector{Float64}, set::AbstractAmbiguitySet, cost::Vector{Float64})
+    #   adds child node to tree.nodes[id]
+    1 <= id <= length(tree.nodes) || throw(BoundsError(tree, id))   # check if id is valid
+    k = get_stage(tree, id) + 1                                     # get new stage value
+    push!(tree.nodes, TreeNode(id, Vector{}(), k, ξ, set, cost ))   # push to node list
+    child_id = length(tree.nodes)                                   # get current node ID
+    push!(tree.nodes[id].children, child_id)                        # push child_id to parent node children
+    if k > tree.K
+        tree.K = k  # update length of tree to the maximum value
+    end
+end
+
+get_children(tree, id) = tree.nodes[id].children
+get_parent(tree,id) = tree.nodes[id].parent
+get_stage(tree, id) = tree.nodes[id].k
+get_scenario(tree, id) = tree.nodes[id].ξ
+
+function get_history(tree::Tree, id::Int)::Array{Int}
+    # gets a vector of tree node IDs up until current
+    stage = get_stage(tree, id)
+    hist = Array{Int}(undef, stage)
+
+    current_id = id
+    for k = stage:-1:1
+        hist[k] = current_id
+        current_id = get_parent(tree, current_id)
+    end
+    return hist
 end
+
+function get_future(tree::Tree, root_id::Int)::Array{Int}
+    #   output list of all leaf node IDs branching from root_id
+    arr_leaves = Int[]
+
+    function iterate_children(tree::Tree, id::Int)
+        children = get_children(tree, id)
+        if length(children) == 0
+            #buffer output
+            push!(arr_leaves, id)
+        else
+            for child in children
+                iterate_children(tree, child)
+            end
+        end
+    end
+
+    iterate_children(tree, root_id)
+    return arr_leaves
+end
+
+function get_stage_id(tree::Tree)::Array{Array{Int}}
+    # gets a list of tree node IDs separated by stages
+    K = tree.K
+    nodelist = [ Int[] for _ in 1:K]
+
+    for id in 1:length(tree.nodes)
+        k = get_stage(tree, id)
+        push!(nodelist[k], id)
+    end
+    return nodelist
+end
+
+include("LagrangeDual.jl")
+
+"""
+function apply!(model, set::AbstractAmbiguitySet) end
+
+function apply!(model, set::WassersteinSet)
+    # add variables
+    # add constraints
+end
+
+mutable struct BundleProblem
+    objective_functions # has a pointer to your function
+    constraints
+    params::Parameter
+end
+"""
+
+end