Add ImportanceSamplingForwardPass

[odow]

using SDDP
import HiGHS
import Random

model = SDDP.LinearPolicyGraph(;
    stages = 3,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = HiGHS.Optimizer,
) do sp, t
    @variable(sp, 5 <= x <= 15, SDDP.State, initial_value = 10)
    @variable(sp, g_t >= 0)
    @variable(sp, g_h >= 0)
    @variable(sp, s >= 0)
    @constraint(sp, balance, x.out - x.in + g_h + s == 0)
    @constraint(sp, demand, g_h + g_t == 0)
    @stageobjective(sp, s + t * g_t)
    SDDP.parameterize(sp, [(0.0, 7.5), (3.0, 5.0), (10.0, 2.5)]) do w
        set_normalized_rhs(balance, w[1])
        set_normalized_rhs(demand, w[2])
        return
    end
    return
end

Random.seed!(1234)

SDDP.train(
    model;
    log_every_iteration = true,
    iteration_limit = 10,
    risk_measure = SDDP.Entropic(1e-1),
    forward_pass = SDDP.ImportanceSamplingForwardPass(),
    cut_type = SDDP.MULTI_CUT,
)

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 41 lines in your changes missing coverage. Please review.

Project coverage is 95.31%. Comparing base (979558d) to head (309376a).

Files with missing lines	Patch %	Lines
src/plugins/forward_passes.jl	0.00%	41 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #848      +/-   ##
==========================================
- Coverage   96.39%   95.31%   -1.08%     
==========================================
  Files          26       26              
  Lines        3633     3674      +41     
==========================================
  Hits         3502     3502              
- Misses        131      172      +41     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[joaquimg]

This is it:
https://arxiv.org/abs/2307.13190
?

[odow]

Ah I knew there was something! A student of Dave is trying to do something similar.

[odow]

💯 This was the next step I had in mind:

They're coming at it from a different direction/motivation, but it's the same core technique.

[Danjiang2000]

Hi Oscar,

For now, if my understanding is correct, the code for ImportanceSamplingForwardPass is like a skeleton?
I notice the following code block:

for child in node.children
    for noise in model[child.term]
        push!(nominal_probability, child.probability * noise.probability)
        push!(support, (child.term, noise.term))
    end
end

The block above seems to be a standard forward pass. I think of a linear policy graph, the child is deterministic and we sample the noise using the original pmf.
The missing part I need to add is the importance sampling pmf (for entropic measure) and weight adjustment for noise realization at every stage (or shall I record weight adjustment in the simulation part?)

[odow]

@Danjiang2000 it's a complete implementation. The two for-loops build the nominal probability pmf, then we use SDDP.adjust_probability to get the risk-adjusted probabilities, and then we sample a new node and noise from that distribution.

[Danjiang2000]

Hi Oscar,

Thanks for your last answer, it is really helpful.

I am 'drafting' to modify this forward_pass function so that it can better serve the purpose of an importance sampling experiment.
In the current implementation, we sample noise using adjusted pmf, and we can find the cumulative cost along the path we sample. To do the importance sampling, I also want to know the weight adjustment for each path, which means I need to multiply the weight adjustment on each path together.
I am not sure if there is a better way to suggest code modification I want to make, so I will write the code in this comment:

At line 413,
cumulative_weight = 1.0

At line 464,
sel = (node_index, noise)
idx = findfirst(x -> x == sel, support)
ratio = nominal_probability[idx] / adjusted_probability[idx]
cumulative_weight *= ratio

At line 475,
cumulative_weight = cumulative_weight,

Then, to simulate via importance sampling, I can write another function that repeat this modified forward pass multiple times. If I want to have a point estimate for entropic cost, I can find point estimate for exponential cost by sample averaging the product of cumulative weight and exponential of cumulative cost.
Does this sound like a good plan if I want to implement a complete importance sampling experiment?

[odow]

Does this sound like a good plan if I want to implement a complete importance sampling experiment?

Yes, I think at least for now.

What this is really asking for the ability to sample with the risk-adjusted probabilities in SDDP.simulate instead of SDDP.train. We either need to write something separate for SDDP.simulate, or change SDDP.simulate to use the AbstractForwardPass machinery.

[joaquimg]

This was the next step I had in mind: ...

I haven't update arXiv, but yet another alternative version for single cut is:

1. cache the cuts separately, besides the aggregate (single) cuts
2. after solving each optimization problem with single cut, check (apply optimized states) in the individual cuts which groups is "more active".

It is a bit weird in the sense that if you are storing all cuts individually, why not just do multicut? Maybe single cut is faster for that problem.

[odow]

Yeah I thought of this. I initially considered looping through the cuts evaluating them before I realized I could reach in and find the right JuMP.value(theta).

I guess the point is that there is a lot of information in the disaggregated cuts that we are throwing away just for the sake of adding fewer constraints.