Question regarding `calculate_bound` with `Bin` variables.

Hello, I'm not sure if this is the right place to ask but I've been looking at this for quite some time with no success.

This is about a very simple problem defined in the code below. What I have trouble understanding is how I am getting 1.2 as a the bound for some objective function that is supposed to be the sum of a binary variable over time (I expected 2, definitely not a rational).

The code contains the context, the deterministic JuMP model and the corresponding deterministic SDDP model. I will also add the two plots that the code produce, which also hints to me that the bound should have been 2, based on the generator curve. I also added the program output at the end.

I'm not very familiar with the library yet, but increasing the number of iterations does still produce a bound of 1.2.

using DataFrames
using Plots
using JuMP
using HiGHS
using SDDP

""" 
From deterministic to stochastic version of MODEL1.

Trying to applying the principles found here:
https://sddp.dev/stable/tutorial/example_reservoir/

This is a simplified version of the original problem to build on it later on.
NOTE: The behaviour is not the same of the original problem. F.e. battery dynamics.

----------------------------------
Consider a network with three devices:
(1) A generator:
    This is a button that when pressed at time t, generates 5 units of energy (UOE)
(2) A demand:
    This is a list of UOE from 1 to 10 at time t. F.e. [2, 3, 10, 1, 6, 7...]
(3) A battery:
    Has a capacity of 100 UOE and starts at 20 UOE.
    If at time t the delta = generated - demand is:
        Positive: then it charges delta UOE
        Negative: then it discharges 2 * delta UOE

We want:
(1) That at no point the charge of the battery goes strictly under 10 UOE
(2) Minimize the number of times we use the generator

----------------------------------
Example:
Demand = [8, 2, 1]

(Up is when the generator is used, down when not)
Battery - 14 ----------- 17 - 21
        - 4 (INVALID)  |    - 15 <- SOLUTION1
                       - 10 - 14 <- SOLUTION2
                            - 8 (INVALID)

The answer is 2 since we need at least to use the generator twice to verify every constraint.
"""

demand = [8, 2, 1]
data = DataFrame(; demand = demand)
T = length(demand)
b_initial = 20
b_min = 10
b_capacity = 100
g_output = 5
M = 4 # we know that max_delta_zero <= (g_output - min_demand) = 5 - 1 = 4
u_g = x_b = max_delta_zero = o_demand = 0 # shutup the linter

"Old model1, for reference"
function run_model()
    model = Model(HiGHS.Optimizer)
    set_silent(model)

    @variable(model, b_min <= x_b[1:T+1] <= b_capacity)
    fix(x_b[1], b_initial; force = true)
    @variable(model, u_g[1:T], Bin)
    @variable(model, max_delta_zero[1:T])
    @variable(model, bin[1:T], Bin)
    for t in 1:T
        delta = g_output * u_g[t] - demand[t]
        @constraint(model, max_delta_zero[t] >= 0)
        @constraint(model, max_delta_zero[t] >= delta)
        @constraint(model, max_delta_zero[t] <= M * bin[t])
        @constraint(model, max_delta_zero[t] <= delta + M * (1 - bin[t]))
        @constraint(model, x_b[t+1] == x_b[t] + 2 * delta - max_delta_zero[t])
    end

    # Finds SOLUTION1
    @objective(model, Min, sum(u_g[t] for t in 1:T))
    optimize!(model)
    ov = objective_value(model)
    @assert ov == 2 "The objective value is $(ov)"

    combined_plot = plot(
        plot(value.(u_g); ylabel = "Generator", legend = false),
        plot(demand; ylabel = "Demand", legend = false),
        plot(value.(max_delta_zero); ylabel = "Max DZ", legend = false),
        plot(g_output * value.(u_g) - demand; ylabel = "Delta", legend = false),
        plot(value.(x_b); ylabel = "Battery", xlabel = "JUMP hour", legend = false);
        layout = (5, 1),
        legend = false,
        size = (400, 600),
    )
    hline!(combined_plot[end-1], [0]; color = :red, linestyle = :dash, label = "Zero")
    savefig(combined_plot, "model1.png")
end

function run_sddp_deterministic_model()
    """ https://sddp.dev/stable/tutorial/example_reservoir/ """
    model = SDDP.LinearPolicyGraph(;
        stages = T,
        sense = :Min,
        lower_bound = 0.0,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        @variable(sp, b_min <= x_b <= b_capacity, SDDP.State, initial_value = b_initial)
        @variable(sp, u_g, Bin)
        @variable(sp, o_demand)
        fix(o_demand, data[t, :demand])
        @variable(sp, max_delta_zero)
        @variable(sp, bin, Bin)

        delta = g_output * u_g - o_demand
        @constraint(sp, max_delta_zero >= 0)
        @constraint(sp, max_delta_zero >= delta)
        @constraint(sp, max_delta_zero <= M * bin)
        @constraint(sp, max_delta_zero <= delta + M * (1 - bin))
        @constraint(sp, x_b.out == x_b.in + 2 * delta - max_delta_zero)
        @stageobjective(sp, u_g)
    end

    # Finds SOLUTION2
    SDDP.train(model; iteration_limit = 10)
    bound = SDDP.calculate_bound(model)
    # Shouldn't this be 2?
    println("bound = $(bound)")

    simulations = SDDP.simulate(
        model,
        1,  # Number of replications
        [:x_b, :u_g, :max_delta_zero],
    )

    g_sim = [sim[:u_g] for sim in simulations[1]]
    mdz_sim = [sim[:max_delta_zero] for sim in simulations[1]]
    b_sim = vcat([b_initial], [sim[:x_b].out for sim in simulations[1]])

    combined_plot = plot(
        plot(g_sim; ylabel = "Generator", legend = false),
        plot(demand; ylabel = "Demand", legend = false),
        plot(mdz_sim; ylabel = "Max DZ", legend = false),
        plot(g_output * g_sim - demand; ylabel = "Delta", legend = false),
        plot(b_sim; ylabel = "Battery", xlabel = "SDDP hour", legend = false);
        layout = (5, 1),
        legend = false,
        size = (400, 600),
    )
    hline!(combined_plot[end-1], [0]; color = :red, linestyle = :dash, label = "Zero")
    savefig(combined_plot, "sddp1.png")
end

function main()
    run_model()
    run_sddp_deterministic_model()
end

main()

And the two plots

model1 sddp1

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 3
  state variables : 1
  scenarios       : 1.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 7]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  AffExpr in MOI.GreaterThan{Float64}     : [2, 2]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.EqualTo{Float64}     : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [2, 2]
  VariableRef in MOI.LessThan{Float64}    : [1, 2]
  VariableRef in MOI.ZeroOne              : [2, 2]
numerical stability report
  matrix range     [1e+00, 1e+01]
  objective range  [1e+00, 1e+00]
  bounds range     [1e+01, 1e+02]
  rhs range        [4e+00, 4e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   2.000000e+00  1.200000e+00  1.912148e+00         6   1
        10   2.000000e+00  1.200000e+00  2.004904e+00        60   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 2.004904e+00
total solves   : 60
best bound     :  1.200000e+00
simulation ci  :  2.000000e+00 ± 0.000000e+00
numeric issues : 0
-------------------------------------------------------------------

bound = 1.2

See https://sddp.dev/stable/guides/add_integrality/

We do not necessarily converge to an optimal solution if there are integer variables.

You could try:

SDDP.train(model; iteration_limit = 10, duality_handler = SDDP.LagrangianDuality())
# or
SDDP.train(model; iteration_limit = 10, duality_handler = SDDP.StrengthenedConicDuality())

See https://sddp.dev/stable/apireference/#Duality-handlers

Thank you for your answer.

The Lagrangian returns 1.28 and the Conic 1.2 so no luck in that regard.

I will try some other approach and see if I can get rid of the integer variables to model this problem.

The Lagrangian returns 1.28 and the Conic 1.2 so no luck in that regard.

This makes sense because of the big-M formulation (even if the M is quite tight).

It has been a while since I wrote https://sddp.dev/stable/guides/add_integrality, and I think it stands up well.

For small problems like this, you're right. Just use JuMP and get a better answer. No debate there. But:

SDDP.jl cannot guarantee that it will find a globally optimal policy when some of the variables are discrete. However, in most cases we find that it can still find an integer feasible policy that performs well in simulation.

Moreover, when the number of nodes in the graph is large, or there is uncertainty, we are not aware of another algorithm that can claim to find a globally optimal policy.

In general, we recommend that you introduce integer variables into your model without fear of the consequences, and that you treat the resulting policy as a good heuristic, rather than an attempt to find a globally optimal policy.

if your problem is too big to solve as a monolithic JuMP problem, give SDDP a try. If it works. Great. If it doesn't. Try a different heuristic. You can't hope to solve the problem to global optimality no matter what algorithm you use.

Closing because this seems resolved. It is not a bug in SDDP.jl, and I think the docs summarize the issue well enough.

Please comment if you have more questions and I can re-open.

I'm sorry I could not find the time to answer before.

It was indeed not a bug. Maybe I should have written it in some other place.

I have not implemented yet the fully fletched model in SDDP, mostly scared about performance since the equivalent JuMP model can already take ~20mins for some scenarios and I read this reply of yours in another user-question issue #775 :

As a rule of thumb, you should expect the SDDP model with uncertainty to take 10^3 - 10^5 times longer than it to solve. If it solves in less than a second, the SDDP will still likely take hours to solve. If it takes longer. Stop.

But now I'm not sure how to understand this part too:

if your problem is too big to solve as a monolithic JuMP problem, give SDDP a try.

Do you mean by monolithic JuMP problem, the model when feeded the (in my case up to ~8 or so, weekly) scenarios to solve simultaneously?

And last

For small problems like this...

I'm not sure how to assess the effective size of the model for our conversation about performance. To be concrete I will post here the current complete version in JuMP. Note that T can be either 168 (a week), 168 x 2 or 168 x 4, but these last two are quite slow at the moment (a scenario is passed here through ipt).

Do you mean by monolithic feeding multiple (lets say ~5) ipts instances to solve simultaneously?

Thank you again for your time.

model = Model(HiGHS.Optimizer)
set_silent(model)

T = length(ipt.demand)

solar = ipt.solar
demand = ipt.demand
added_demand = 0.1
demand .+= added_demand
grid = ipt.outage
diesel_power = ipt.info.diesel_power
ccoe = ipt.info.ccoe
dcoe = ipt.info.dcoe
b_initial = ipt.info.b_initial
b_min = ipt.info.b_min
capacity = ipt.info.capacity
grid_power = ipt.info.grid_power

coef = div(60, granularity_data)
w1 = ipt.info.w1
w2 = ipt.info.w2 * coef
w3 = ipt.info.w3 * coef
w4 = ipt.info.w4 * coef

# M parameter for the BigM technique
max_possible_produced = maximum([grid[t] + solar[t] - demand[t] for t in 1:T])
pessimist = diesel_power + max_possible_produced
M = pessimist # ~ 30
# Take the smallest M bigger than any possible b_t
M2 = 2

@variable(model, b_min <= b[1:T+1] <= 1)
fix(b[1], b_initial; force = true)

@variable(model, u_diesel[1:T], Bin)
@variable(model, u_grid[1:T], Bin)
@variable(model, big_m_helper[1:T], Bin)
@variable(model, max_delta_zero[1:T])
@variable(model, big_m_helper2[1:T], Bin)
@variable(model, min_trans_one[1:T])

# Track the amount of times we started the diesel generator
@variable(model, diesel_starts[1:T], Bin)
@constraint(model, [t in 2:T], diesel_starts[t] >= u_diesel[t] - u_diesel[t-1])
@constraint(model, diesel_starts[1] >= u_diesel[1])

# Track the maximum amount of time we used the diesel generator
# Uses McCormick envelope constraints
@variable(model, run_time[1:T] >= 0)
@variable(model, max_run_time[1:T] >= 0)
@constraint(model, run_time[1] == 0)
@constraint(model, max_run_time[1] == run_time[1])
@variable(model, bilin_prod[1:T] >= 0) # This is run_time * u_diesel

M3 = 25
xL = 0
xU = 1
yL = 0
yU = M3 # <= max possible of run_time
for t in 1:T
    @constraints model begin
        bilin_prod[t] >= xL * run_time[t] + yL * u_diesel[t] - xL * yL
        bilin_prod[t] >= xU * run_time[t] + yU * u_diesel[t] - xU * yU
        bilin_prod[t] <= xU * run_time[t] + yL * u_diesel[t] - xU * yL
        bilin_prod[t] <= xL * run_time[t] + yU * u_diesel[t] - xL * yU
    end
end
for t in 2:T
    @constraints(model, begin
        run_time[t] == bilin_prod[t-1] + u_diesel[t-1]
        max_run_time[t] >= run_time[t]
        max_run_time[t] >= max_run_time[t-1]
    end)
end

# -------------- Heuristics
# NOTE: There is always a pick in battery at the start, so to increase b_min
@constraint(model, [t in 15:T+1], b[t] >= 0.3)
# NOTE: In solutions without this constraint it seems to be verified anyway
for t in 1:T
    # Never use the grid and the diesel at the same time
    @constraint(model, u_grid[t] + u_diesel[t] <= 1)

    # Always use grid if possible
    fix(u_grid[t], 1; force = true)
end

# -------------- Heuristics for u_diesel?

for t in 1:T
    if grid[t] == 0
        fix(u_grid[t], 0; force = true) # force to not use it if outage
    end
    delta = grid_power * u_grid[t] + diesel_power * u_diesel[t] + solar[t] - demand[t]

    # Big M technique to get max(delta, 0)
    @constraints(model, begin
        max_delta_zero[t] >= 0
        max_delta_zero[t] >= delta
        max_delta_zero[t] <= M * big_m_helper[t]
        max_delta_zero[t] <= delta + M * (1 - big_m_helper[t])
    end)
    # coe = delta > 0 ? ccoe * delta : dcoe * delta
    coe = dcoe * delta - (dcoe - ccoe) * max_delta_zero[t]

    # NOTE: The battery can not go over 100% => take the min(b_(t+1), 1)
    trans = b[t] + coe / (granularity_data * capacity) # transition

    # Big M technique to get min(b_t, 1)
    @constraints(model, begin
        min_trans_one[t] <= 1
        min_trans_one[t] <= trans
        min_trans_one[t] >= 1 - M2 * big_m_helper2[t]
        min_trans_one[t] >= trans - M2 * (1 - big_m_helper2[t])
    end)

    @constraint(model, b[t+1] == min_trans_one[t])
end

@objective(
    model,
    Min,
    w1 * sum(diesel_starts[t] for t in 1:T) +
    w2 * sum(u_diesel[t] for t in 1:T) +
    w3 * max_run_time[T] +
    w4 * sum(u_grid[t] for t in 1:T)
)

optimize!(model)

Do you mean by monolithic JuMP problem,

I mean a single scenario.

Well I guess the only thing left is for me to actually write the code and see.

I think you can close this now, this time for good, and thank you again!

Well I guess the only thing left is for me to actually write the code and see.

Yip :smile:

thank you again

No problem

odow / SDDP.jl

Question regarding `calculate_bound` with `Bin` variables. #776