Cut-sharing property

[Esnilg]

Hello Oscar,

I hope you're well. I am writing to you because I have a question regarding SDDP.LinearPolicyGraph and SDDP.MarkovianPolicyGraph. Suppose I have a problem where there is uncertainty in the coefficients of the objective function. I can use SDDP.LinearPolicyGraph, with 3 realizations in each stage. Here the Cut-sharing property is applied. But when I use SDDP.MarkovianPolicyGraph with a transition matrix of equal probability (that is, independent with 3 markov states per stage). Here the Cut-sharing property does not apply. certain?

[odow]

I strongly dislike "cut-sharing" as a concept. I find it unhelpful, confusing, and misleading.

The linear case is not "sharing" any cuts between the 3 realizations. We're building up a single value function at each stage. In the Markovian case, we build up three separate value functions. Cut sharing arises from the literature which considers a scenario tree. We do not use scenario trees.

This theory tutorial manages to implement an infinite horizon SDDP algorithm that can solve both the linear and Markovian policy graph cases without ever mentioning cut-sharing or scenario trees: https://odow.github.io/SDDP.jl/stable/tutorial/theory/21_theory_intro/

[Esnilg]

Okay, thanks, I'll read what you send me. Greetings

[Esnilg]

Hello Oscar,

I already read the tutorial and now I understand why you don't use the term Cut-sharing. But I want to rephrase the question, when I solve a problem with SDDP.MarkovianPolicyGraph with an independent transition matrix (suppose 3 markov states in each stage with equal probability of transition in the target function costs) the solution SDDP.objective_bound ( model) should be the same as SDDP.LinearPolicyGraph using the same characteristics?

[odow]

Yes, the solution should be the same if you have a linear policy graph or a Markovian policy graph with independent transition probabilities.

If you have an example where that is not the case, that's a bug :)

[Esnilg]

Thanks Oscar

Yesterday I can verify with a simple example using SDDP with linear policy and with markov chain. I also made an example using an AR (1) for the demand. I have tried several ways but I keep getting different solutions. I send you the codes. What am I doing wrong?

Best regards

#####################################
#
#  SDDP with AR(1): when alpha = 1 then AR.out = demand (case independent)
#
#####################################

using SDDP, LinearAlgebra, GLPK, Test

    build_cost = 1e4
    use_cost = 4
    num_units = 5
    capacities = ones(num_units)
    # Cost of unmet demand
    penalty = 5e5
    # Discounting rate
    rho = 0.01
    T = 3

    a = 1
    b = 3
    alpha = 1
    model = SDDP.LinearPolicyGraph(
        stages = T,
        sense = :Min,
        lower_bound = 0.0,
        optimizer = GLPK.Optimizer,
    ) do sp, stage
        @variable(
            sp,
            0 <= invested[1:num_units] <= 5,
            SDDP.State,
            initial_value = 0
        )

        @variable(
            sp,
            AR >= 0,
            SDDP.State,
            initial_value = 2
        )

        @variables(sp, begin
            generation >= 0
            unmet >= 0
            demand
        end)

        @constraint(sp, estado, AR.out == (1-alpha)*AR.in + demand);

        if stage < T
            @constraints(
                sp,
                begin
                    # Can't un-invest
                    investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # # Generation capacity
                    # sum(capacities[i] * invested[i].out for i in 1:num_units) >= generation
                    # # Meet demand or pay a penalty
                    # unmet >= demand - sum(generation)
                    # # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        if stage > 1
            @constraints(
                sp,
                begin
                    # # Can't un-invest
                    # investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # Generation capacity
                    sum(capacities[i] * invested[i].in for i in 1:num_units) >= generation
                    # Meet demand or pay a penalty
                    unmet >= AR.in - sum(generation)
                    # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        # Parameterize the subproblem.
        Ω = [
        (demand_vals = 0, fuel_multiplier = 1.5),
        (demand_vals = 0.5, fuel_multiplier = 1.0),
        (demand_vals = 1, fuel_multiplier = 0.75),
        ]

        @expression(
            sp,
            investment_cost,
            build_cost *
            sum(invested[i].out - invested[i].in for i in 1:num_units)
        )

        if stage == 1
            SDDP.parameterize(sp, Ω, [1 / 3, 1 / 3, 1 / 3]) do ω
                JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
                @stageobjective(sp, investment_cost * rho^(stage - 1)  )
            end
        else
            SDDP.parameterize(sp, Ω, [1 / 3, 1 / 3, 1 / 3]) do ω
                JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
                @stageobjective(sp, ( stage < T ? (investment_cost) * rho^(stage - 1) : 0) + (generation * use_cost * ω.fuel_multiplier ) * rho^(stage - 1) + penalty * unmet)
            end
        end 
    end

   SDDP.train(
        model,
        iteration_limit = 500,
        log_frequency = 10,
    )

simulacion = SDDP.simulate(model, 100, [:invested,:AR]);

SDDP.parameterize(model[1], (demand_vals = 0, fuel_multiplier = 1.5))

SDDP.write_subproblem_to_file(model[1], "subproblem.lp")

read("subproblem.lp") |> String |> print

#################################################################################
#####################################
#
#  SDDP simple
#
#####################################

using SDDP, LinearAlgebra, GLPK, Test

    build_cost = 1e4
    use_cost = 4
    num_units = 5
    capacities = ones(num_units)
    # Cost of unmet demand
    penalty = 5e5
    # Discounting rate
    rho = 0.01
    T = 3
    model = SDDP.LinearPolicyGraph(
        stages = T,
        sense = :Min,
        lower_bound = 0.0,
        optimizer = GLPK.Optimizer,
    ) do sp, stage
        @variable(
            sp,
            0 <= invested[1:num_units] <= 5,
            SDDP.State,
            initial_value = 0
        )
        @variables(sp, begin
            generation >= 0
            unmet >= 0
            demand
        end)

        if stage < T
            @constraints(
                sp,
                begin
                    # Can't un-invest
                    investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # # Generation capacity
                    # sum(capacities[i] * invested[i].out for i in 1:num_units) >= generation
                    # # Meet demand or pay a penalty
                    # unmet >= demand - sum(generation)
                    # # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        if stage > 1
            @constraints(
                sp,
                begin
                    # # Can't un-invest
                    # investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # Generation capacity
                    sum(capacities[i] * invested[i].in for i in 1:num_units) >= generation
                    # Meet demand or pay a penalty
                    unmet >= demand - sum(generation)
                    # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        # Parameterize the subproblem.
        Ω = [
        (demand_vals = 1, fuel_multiplier = 1.5),
        (demand_vals = 2, fuel_multiplier = 1.0),
        (demand_vals = 3, fuel_multiplier = 0.75),
        ]

        @expression(
            sp,
            investment_cost,
            build_cost *
            sum(invested[i].out - invested[i].in for i in 1:num_units)
        )

        if stage == 1
            SDDP.parameterize(sp, Ω, [1 / 3, 1 / 3, 1 / 3]) do ω
                JuMP.fix(demand, ω.demand_vals)
                @stageobjective(sp, investment_cost * rho^(stage - 1)  )
            end
        else
            SDDP.parameterize(sp, Ω, [1 / 3, 1 / 3, 1 / 3]) do ω
                JuMP.fix(demand, ω.demand_vals)
                @stageobjective(sp, ( stage < T ? (investment_cost) * rho^(stage - 1) : 0) + (generation * use_cost * ω.fuel_multiplier ) * rho^(stage - 1) + penalty * unmet)
            end
        end 
    end

   SDDP.train(
        model,
        iteration_limit = 500,
        log_frequency = 10,
    )

simulacion = SDDP.simulate(model, 10, [:invested]);

[odow]

Both of these use linear policy graphs. Do you have the code for the Markovian?

[Esnilg]

Hi, here is the code for Markovian and Markovian with AR process. All models should give the same solution, but with AR.out = demand, it does not happen. Best Regards

############################################################3

using SDDP, LinearAlgebra, GLPK, Test

    build_cost = 1e4
    use_cost = 4
    num_units = 5
    capacities = ones(num_units)
    # Cost of unmet demand
    penalty = 5e5
    # Discounting rate
    rho = 0.01
    T = 3

    datos  = [0.75,1.0,1.5]
    Z_aux = Dict(i => datos for i in 1:5)

    model = SDDP.MarkovianPolicyGraph(
    transition_matrices = Array{Float64,2}[
        [1.0]',
        [0.3 0.3 0.3],
        [0.3 0.3 0.3; 0.3 0.3 0.3;0.3 0.3 0.3],
    ],
    sense = :Min,
    lower_bound = 0.0,
    optimizer = GLPK.Optimizer,
    ) do sp, node
    # Unpack the stage and Markov index.
        stage, markov_state = node 

        @variable(
            sp,
            0 <= invested[1:num_units] <= 4,
            SDDP.State,
            initial_value = 0
        )

        @variables(sp, begin
            generation >= 0
            unmet >= 0
            demand
        end)

        if stage < T
            @constraints(
                sp,
                begin
                    # Can't un-invest
                    investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # # Generation capacity
                    # sum(capacities[i] * invested[i].out for i in 1:num_units) >= generation
                    # # Meet demand or pay a penalty
                    # unmet >= demand - sum(generation)
                    # # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        if stage > 1
            @constraints(
                sp,
                begin
                    # # Can't un-invest
                    # investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # Generation capacity
                    sum(capacities[i] * invested[i].in for i in 1:num_units) >= generation
                    # Meet demand or pay a penalty
                    unmet >= demand - sum(generation)
                    # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        # Parameterize the subproblem.
        Ω = [
        (demand_vals = 1, fuel_multiplier = Z_aux[stage][markov_state]),
        (demand_vals = 2, fuel_multiplier = Z_aux[stage][markov_state]),
        (demand_vals = 3, fuel_multiplier = Z_aux[stage][markov_state]),
        ]

        @expression(
            sp,
            investment_cost,
            build_cost *
            sum(invested[i].out - invested[i].in for i in 1:num_units)
        )

        if stage == 1
            SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
                JuMP.fix(demand, ω.demand_vals)
                @stageobjective(sp, investment_cost * rho^(stage - 1)  )
            end
        else
            SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
                #print("$(ω.demand_vals) $(stage) $(markov_state)\n")
                JuMP.fix(demand, ω.demand_vals)
                @stageobjective(sp, ( stage < T ? (investment_cost) * rho^(stage - 1) : 0) + (generation * use_cost * ω.fuel_multiplier ) * rho^(stage - 1) + penalty * unmet)
            end
        end 
    end

   SDDP.train(
        model,
        iteration_limit = 500,
        log_frequency = 10,
    )

simulacion = SDDP.simulate(model, 1000, [:invested]);

##############################################

using SDDP, LinearAlgebra, GLPK, Test

    build_cost = 1e4
    use_cost = 4
    num_units = 5
    capacities = ones(num_units)
    # Cost of unmet demand
    penalty = 5e5
    # Discounting rate
    rho = 0.01
    T = 3
    a = 1
    b = 3
    alpha = 1

    datos  = [0.75,1.0,1.5]
    Z_aux = Dict(i => datos for i in 1:5)

    model = SDDP.MarkovianPolicyGraph(
    transition_matrices = Array{Float64,2}[
        [1.0]',
        [0.3 0.3 0.3],
        [0.3 0.3 0.3; 0.3 0.3 0.3;0.3 0.3 0.3],
    ],
    sense = :Min,
    lower_bound = 0.0,
    optimizer = GLPK.Optimizer,
    ) do sp, node
    # Unpack the stage and Markov index.
        stage, markov_state = node 

        @variable(
            sp,
            0 <= invested[1:num_units] <= 4,
            SDDP.State,
            initial_value = 0
        )

        @variables(sp, begin
            generation >= 0
            unmet >= 0
            demand
        end)

        @variable(
            sp,
            AR >= 0,
            SDDP.State,
            initial_value = 2
        )

        @constraint(sp, estado, AR.out == (1-alpha)*AR.in + demand);

        if stage < T
            @constraints(
                sp,
                begin
                    # Can't un-invest
                    investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # # Generation capacity
                    # sum(capacities[i] * invested[i].out for i in 1:num_units) >= generation
                    # # Meet demand or pay a penalty
                    # unmet >= demand - sum(generation)
                    # # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        if stage > 1
            @constraints(
                sp,
                begin
                    # # Can't un-invest
                    # investment[i in 1:num_units], invested[i].out >= invested[i].in
                    # Generation capacity
                    sum(capacities[i] * invested[i].in for i in 1:num_units) >= generation
                    # Meet demand or pay a penalty
                    unmet >= AR.in - sum(generation)
                    # For fewer iterations order the units to break symmetry, units are identical (tougher numerically)
                    # [j in 1:(num_units - 1)], invested[j].out <= invested[j + 1].out
                end
            )
        end

        # Parameterize the subproblem.
        Ω = [
        (demand_vals = 0, fuel_multiplier = Z_aux[stage][markov_state]),
        (demand_vals = 0.5, fuel_multiplier = Z_aux[stage][markov_state]),
        (demand_vals = 1, fuel_multiplier = Z_aux[stage][markov_state]),
        ]

        @expression(
            sp,
            investment_cost,
            build_cost *
            sum(invested[i].out - invested[i].in for i in 1:num_units)
        )

        if stage == 1
            SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
                JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
                @stageobjective(sp, investment_cost * rho^(stage - 1)  )
            end
        else
            SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
                #print("$(ω.demand_vals) $(stage) $(markov_state)\n")
                JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
                @stageobjective(sp, ( stage < T ? (investment_cost) * rho^(stage - 1) : 0) + (generation * use_cost * ω.fuel_multiplier ) * rho^(stage - 1) + penalty * unmet)
            end
        end 
    end

   SDDP.train(
        model,
        iteration_limit = 500,
        log_frequency = 10,
    )

simulacion = SDDP.simulate(model, 1000, [:invested]);

[odow]

Your linear policy graph is not quite the same as the Markovian one.

Here's an example of a how to turn a linear policy graph into the equivalent Markovian graph. Note that in the linear case we use parameterize, but in the Markovian case we don't need to.

Your Markovian models use both a Markovian graph and parameterize (which is allowed and useful -- see https://odow.github.io/SDDP.jl/stable/tutorial/basic/04_markov_uncertainty/ -- but I don't think it results in exactly the same model).

I suggest you step back and draw the policy graph on paper, label each node, and describe what information flows into each node on a wavy arc. Your code is a bit confusing, so I don't fully understand what you are trying to do.

using SDDP

Ω = [
    (inflow = 0.0, fuel_multiplier = 1.5),
    (inflow = 50.0, fuel_multiplier = 1.0),
    (inflow = 100.0, fuel_multiplier = 0.75),
]

# Here is the linear model

model = SDDP.LinearPolicyGraph(
    stages = 3,
    lower_bound = 0.0,
) do subproblem, t
    @variable(subproblem, 0 <= volume <= 200, SDDP.State, initial_value = 200)
    @variables(subproblem, begin
        thermal_generation >= 0
        hydro_generation >= 0
        hydro_spill >= 0
        inflow
    end)
    @constraints(
        subproblem,
        begin
            volume.out == volume.in + inflow - hydro_generation - hydro_spill
            thermal_generation + hydro_generation == 150.0
        end
    )
    if t == 1
        # Note that the first stage is deterministic, and uses the first sample.
        JuMP.fix(inflow, Ω[1].inflow; force = true)
        @stageobjective(
            subproblem,
            Ω[1].fuel_multiplier * fuel_cost[t] * thermal_generation
        )
    else
        # Parameterize with two-arguments assumes uniform probability.
        SDDP.parameterize(subproblem, Ω) do ω
            JuMP.fix(inflow, ω.inflow; force = true)
            @stageobjective(
                subproblem,
                ω.fuel_multiplier * fuel_cost[t] * thermal_generation
            )
        end
    end
end

# Here is the Markovian model

model = SDDP.MarkovianPolicyGraph(
    transition_matrices = Array{Float64,2}[
        fill(1.0, 1, 1),
        fill(1/3, 1, 3),
        fill(1/3, 3, 3),
    ],
    lower_bound = 0.0,
) do subproblem, node
    t, markov_state = node
    @variable(subproblem, 0 <= volume <= 200, SDDP.State, initial_value = 200)
    @variables(subproblem, begin
        thermal_generation >= 0
        hydro_generation >= 0
        hydro_spill >= 0
        inflow
    end)
    @constraints(
        subproblem,
        begin
            volume.out == volume.in + inflow - hydro_generation - hydro_spill
            thermal_generation + hydro_generation == 150.0
        end
    )
    # Note we don't need to parameterize!
    JuMP.fix(inflow, Ω[markov_state].inflow)
    @stageobjective(
        subproblem,
        Ω[markov_state].fuel_multiplier * fuel_cost[t] * thermal_generation
    )
end

[Esnilg]

Hello Oscar

Thank you for your prompt response, sorry that I did not explain the problem to you first.

The problem has 3 stages and 2 investment periods. In the first stage I only make the first investment decision in capacity, then in the second stage I make the second decision in new capacity but it is also decided how the system is going to operate given the realization of the operating costs, the AR.in and invested.in . In the last stage, I did not make an investment decision but only decided the operation of the system with the operating costs, the invested.in and AR.in.

The model I sent you tries to describe that the operating costs are managed by 3 markov chain states, and the demand follows an AR (1) process that has three independent error realizations. Therefore what I try is that for each (t, markov_state) there are three realizations of the errors with probability p = 1/3

@constraint(sp, estado, AR.out == (1-alpha)*AR.in + demand); # demand = error

Z_aux Dict{Int64,Array{Float64,1}} with 5 entries: 4 => [0.75, 1.0, 1.5] 2 => [0.75, 1.0, 1.5] 3 => [0.75, 1.0, 1.5] 5 => [0.75, 1.0, 1.5] 1 => [0.75, 1.0, 1.5]

Ω = [ (demand_vals = 0, fuel_multiplier = Z_aux[stage][markov_state]), (demand_vals = 0.5, fuel_multiplier = Z_aux[stage][markov_state]), (demand_vals = 1, fuel_multiplier = Z_aux[stage][markov_state]), ]

for example if t = 1, markov_state = 1

Ω = [ (demand_vals = 0, fuel_multiplier =0.75), (demand_vals = 0.5, fuel_multiplier = 0.75), (demand_vals = 1, fuel_multiplier = 0.75), ]

I parameterize the errors (in this case demand) of the AR process (1)

if stage == 1

in the first stage I have no operating costs

        SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
            JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
            @stageobjective(sp, investment_cost * rho^(stage - 1)  )
        end
    else
                    #### in the last stage I have no investment costs
        SDDP.parameterize(sp, Ω, [1/3, 1/3, 1/3]) do ω
            #print("$(ω.demand_vals) $(stage) $(markov_state)\n")
            JuMP.fix(demand, a - (1-alpha)*a +  alpha*(b-a)*ω.demand_vals)
            @stageobjective(sp, ( stage < T ? (investment_cost) * rho^(stage - 1) : 0) + (generation * use_cost * ω.fuel_multiplier ) * rho^(stage - 1) + penalty * unmet)
        end

end

Please correct me if I am doing something wrong in modeling uncertainty.

[odow]

the operating costs are managed by 3 markov chain states, and the demand follows an AR (1) process that has three independent error realizations

Okay. This is good.

The linear model gets a different solution because it models the fuel cost and demand as dependent random variables. You want to use a multidimensional random variable so that there are 3x3 = 9 realizations of the uncertainty at each stage, not 3: https://odow.github.io/SDDP.jl/stable/guides/add_multidimensional_noise_Terms/

[Esnilg]

Hello Oscar

I am super clear. Very grateful for your valuable help !!!

[odow]

No problem. Modeling the timing and dependence of random variables in SDDP is quite complicated. I found drawing out a policy graph of what happens to be enormously useful. It helps you identify differences between Markovian graph structure and the node-wise independent random variables that arrive via the squiggly arcs.