Retrieving information each iteration.

[Richardzhang0416]

Hi Oscar,

We are trying to retrieve the following information for a minimisation problem at each iteration:

1. lower bound
2. statistical upper bound
3. Individual sampled scenario costs
4. Number of scenarios sampled

I tried to set the parameter log_frequency to 1 but it does not print the information each iteration. My understanding is that if this parameter works, information 1 and 2 can be retrieved from the log. If there any simple way to retrieve 3 and 4 as well?

Thanks!

Hongyu

[odow]

You're looking for SDDP.train(model; log_every_iteration = true).

You don't need to print the log though.

Just look at model.most_recent_training_results. This can tell you:

SDDP.jl/src/user_interface.jl

Lines 680 to 689 in 014e054

	struct Log
	iteration::Int
	bound::Float64
	simulation_value::Float64
	time::Float64
	pid::Int
	total_solves::Int
	duality_key::String
	serious_numerical_issue::Bool
	end

For (3): we don't store the scenario or the individual stage costs in the log

For (4): it is always 1. We never sample more than one scenario per iteration

What, precisely, do you need, when and why? There is the undocumented forward_pass_callback argument to SDDP.train that might work.

[odow]

julia> using SDDP, HiGHS

julia> model = SDDP.LinearPolicyGraph(;
           stages = 3,
           lower_bound = 0.0,
           optimizer = HiGHS.Optimizer,
       ) do sp, stage
           @variable(
               sp,
               0 <= stored_production <= 100,
               Int,
               SDDP.State,
               initial_value = 0
           )
           @variable(sp, 0 <= production <= 200, Int)
           @variable(sp, overtime >= 0, Int)
           @variable(sp, demand)
           DEMAND = [[100.0], [100.0, 300.0], [100.0, 300.0]]
           SDDP.parameterize(ω -> JuMP.fix(demand, ω), sp, DEMAND[stage])
           @constraint(
               sp,
               stored_production.out ==
               stored_production.in + production + overtime - demand
           )
           @stageobjective(
               sp,
               100 * production + 300 * overtime + 50 * stored_production.out
           )
       end
A policy graph with 3 nodes.
 Node indices: 1, 2, 3


julia> iterations = Any[]
Any[]

julia> function my_callback(result)
           push!(iterations, result)
           return
       end
my_callback (generic function with 1 method)

julia> SDDP.train(model; forward_pass_callback = my_callback)
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-25
-------------------------------------------------------------------
problem
  nodes           : 3
  state variables : 1
  scenarios       : 4.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [6, 6]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  VariableRef in MOI.EqualTo{Float64}     : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [4, 4]
  VariableRef in MOI.Integer              : [3, 3]
  VariableRef in MOI.LessThan{Float64}    : [2, 3]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [1e+00, 3e+02]
  bounds range     [1e+02, 2e+02]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   3.000000e+04  6.250000e+04  1.711211e-01         8   1
        20   6.000000e+04  6.250000e+04  1.145327e+00       172   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 1.145327e+00
total solves   : 172
best bound     :  6.250000e+04
simulation ci  :  6.450000e+04 ± 9.641401e+03
numeric issues : 0
-------------------------------------------------------------------

julia> iterations
20-element Vector{Any}:
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 0.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 30000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 40000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 60000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 60000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 95000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 40000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 40000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 60000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 300.0), (3, 100.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 0.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 55000.0)
 (scenario_path = Tuple{Int64, Any}[(1, 100.0), (2, 100.0), (3, 300.0)], sampled_states = [Dict(:stored_production => 100.0), Dict(:stored_production => 100.0), Dict(:stored_production => 0.0)], objective_states = Tuple{}[], belief_states = Tuple{Int64, Dict{Int64, Float64}}[], cumulative_value = 60000.0)

[Richardzhang0416]

Hi Oscar, thanks for your reply and suggestions. The argument and function from your suggestions are very helpful. Below are some follow up:

For (3): we don't store the scenario or the individual stage costs in the log
I think this information is in the callback function output? I am taking a closer look at the call back output and figure out what are the elements related to this. Just to clarify because I might not be very clear in my previous question: we need the cost of each scenario and do not need the costs at each stage.

For (4): it is always 1. We never sample more than one scenario per iteration.
Thanks. I see. Is there a way to change some parameter to make it sample multiple scenarios? Or it is hard coded that only one scenario is sampled.

What, precisely, do you need, when and why? There is the undocumented forward_pass_callback argument to SDDP.train that might work.
Ok, what we need is: at the end of each SDDP iteration, the lower bound, statistical upper bound, and the number of scenarios sampled in the current iteration. This is of course based on the case where we sample multiple scenarios each iteration.

We need this information to improve the stopping criteria of SDDP.

Thanks!

Best,

Hongyu

[odow]

what we need is: at the end of each SDDP iteration, the lower bound,

This is log.bound

statistical upper bound,

this is log.simulation_value. But note that we do not compute or store the statistical upper bound as you might think about it. We store only the total cost of the single scenario that was iterated on the forward pass.

and the number of scenarios sampled in the current iteration

This is only ever one.

Is there a way to change some parameter to make it sample multiple scenarios? Or it is hard coded that only one scenario is sampled.

It's hard-coded. It is only ever one scenario. We do not implement the "classical" algorithm with N forward paths.

You can write your own stopping rule which performs a simulation of the policy to estimate an upper bound. See:

SDDP.jl/src/plugins/stopping_rules.jl

Lines 42 to 178 in c37a4f5

	"""
	Statistical(;
	num_replications::Int,
	iteration_period::Int = 1,
	z_score::Float64 = 1.96,
	verbose::Bool = true,
	disable_warning::Bool = false,
	)
	Perform an in-sample Monte Carlo simulation of the policy with
	`num_replications` replications every `iteration_period`s and terminate if the
	deterministic bound (lower if minimizing) falls into the confidence interval for
	the mean of the simulated cost.
	If `verbose = true`, print the confidence interval.
	If `disable_warning = true`, disable the warning telling you not to use this
	stopping rule (see below).
	## Why this stopping rule is not good
	This stopping rule is one of the most common stopping rules seen in the
	literature. Don't follow the crowd. It is a poor choice for your model, and
	should be rarely used. Instead, you should use the default stopping rule, or use
	a fixed limit like a time or iteration limit.
	To understand why this stopping rule is a bad idea, assume we have conducted
	`num_replications` simulations and the objectives are in a vector
	`objectives::Vector{Float64}`.
	Our mean is `μ = mean(objectives)` and the half-width of the confidence interval
	is `w = z_score * std(objectives) / sqrt(num_replications)`.
	Many papers suggest terminating the algorithm once the deterministic bound
	(lower if minimizing, upper if maximizing) is contained within the confidence
	interval. That is, if `μ - w <= bound <= μ + w`. Even worse, some papers define
	an optimization gap of `(μ + w) / bound` (if minimizing) or `(μ - w) / bound`
	(if maximizing), and they terminate once the gap is less than a value like 1%.
	Both of these approaches are misleading, and more often than not, they will
	result in terminating with a sub-optimal policy that performs worse than
	expected. There are two main reasons for this:
	1. The half-width depends on the number of replications. To reduce the
	computational cost, users are often tempted to choose a small number of
	replications. This increases the half-width and makes it more likely that
	the algorithm will stop early. But if we choose a large number of
	replications, then the computational cost is high, and we would have been
	better off to run a fixed number of iterations and use that computational
	time to run extra training iterations.
	2. The confidence interval assumes that the simulated values are normally
	distributed. In infinite horizon models, this is almost never the case. The
	distribution is usually closer to exponential or log-normal.
	There is a third, more technical reason which relates to the conditional
	dependence of constructing multiple confidence intervals.
	The default value of `z_score = 1.96` corresponds to a 95% confidence
	interval. You should interpret the interval as "if we re-run this simulation
	100 times, then the true mean will lie in the confidence interval 95 times
	out of 100." But if the bound is within the confidence interval, then we know
	the true mean cannot be better than the bound. Therfore, there is a more than
	95% chance that the mean is within the interval.
	A separate problem arises if we simulate, find that the bound is outside the
	confidence interval, keep training, and then re-simulate to compute a new
	confidence interval. Because we will terminate when the bound enters the
	confidence interval, the repeated construction of a confidence interval means
	that the unconditional probability that we terminate with a false positive is
	larger than 5% (there are now more chances that the sample mean is optimistic
	and that the confidence interval includes the bound but not the true mean). One
	fix is to simulate with a sequentially increasing number of replicates, so that
	the unconditional probability stays at 95%, but this runs into the problem of
	computational cost. For more information on sequential sampling, see, for
	example, Güzin Bayraksan, David P. Morton, (2011) A Sequential Sampling
	Procedure for Stochastic Programming. Operations Research 59(4):898-913.
	"""
	struct Statistical <: AbstractStoppingRule
	num_replications::Int
	iteration_period::Int
	z_score::Float64
	verbose::Bool
	function Statistical(;
	num_replications,
	iteration_period = 1,
	z_score = 1.96,
	verbose = true,
	disable_warning::Bool = false,
	)
	if !disable_warning
	@warn(
	"Are you really sure you want to use this stopping rule? " *
	"Read why we don't recommend it by typing `? " *
	"SDDP.Statistical` into the REPL to read the docstring.\n\n" *
	"If you understand what you are doing, you can disable this " *
	"warning with `SDDP.Statistical(; disable_warning = true)`",
	)
	end
	return new(num_replications, iteration_period, z_score, verbose)
	end
	end
	stopping_rule_status(::Statistical) = :statistical
	function convergence_test(
	graph::PolicyGraph,
	log::Vector{Log},
	rule::Statistical,
	)
	if length(log) % rule.iteration_period != 0
	# Only run this convergence test every rule.iteration_period iterations.
	return false
	end
	results = simulate(graph, rule.num_replications)
	objectives =
	map(simulation -> sum(s[:stage_objective] for s in simulation), results)
	sample_mean = Statistics.mean(objectives)
	sample_ci =
	rule.z_score * Statistics.std(objectives) / sqrt(rule.num_replications)
	if rule.verbose
	println(
	"Simulated policy value: [",
	print_value(sample_mean - sample_ci),
	", ",
	print_value(sample_mean + sample_ci),
	"]",
	)
	end
	current_bound = log[end].bound
	if graph.objective_sense == MOI.MIN_SENSE
	return sample_mean - sample_ci <= current_bound
	else
	@assert graph.objective_sense == MOI.MAX_SENSE
	return current_bound <= sample_mean + sample_ci
	end
	end

Part of the design of SDDP.jl is that the stopping rules are separate from the training iterations. This tutorial is pretty close to the true algorithm that we implement: https://sddp.dev/stable/explanation/theory_intro/. It isn't useful to think in terms of the Pereira or Shapiro type papers.

See perhaps Section 4 of https://github.com/odow/SDDP.jl/blob/master/papers/policy_graph/preprint.pdf, section 3 of https://github.com/odow/SDDP.jl/blob/master/papers/partially_observable/preprint.pdf