Central Path Cutting Plane

[FSchmidtDIW]

Hi Oscar,

I was wondering if anyone has implemented the Central Path Cutting Plane algorithm based on Chebyshev centres by Beltran et al. for SDDP.jl. I was going to have a go at a new ForwardPass but wanted to make sure that it is not implemented anywhere yet.

Thanks a lot

Felix

[odow]

I have not implemented it. All code for SDDP.jl is contained in this repo.

[odow]

Hi @FSchmidtDIW, how'd you get on?

Here's a place that might be useful to start

SDDP.jl/src/plugins/forward_passes.jl

Lines 323 to 372 in 3b2570b

	mutable struct RegularizedForwardPass{T<:AbstractForwardPass} <:
	AbstractForwardPass
	forward_pass::T
	trial_centre::Dict{Symbol,Float64}
	ρ::Float64
	function RegularizedForwardPass(;
	rho::Float64 = 0.05,
	forward_pass::AbstractForwardPass = DefaultForwardPass(),
	)
	centre = Dict{Symbol,Float64}()
	return new{typeof(forward_pass)}(forward_pass, centre, rho)
	end
	end
	function forward_pass(
	model::PolicyGraph,
	options::Options,
	fp::RegularizedForwardPass,
	)
	if length(model.root_children) != 1
	error(
	"RegularizedForwardPass cannot be applied because first-stage is " *
	"not deterministic",
	)
	end
	node = model[model.root_children[1].term]
	if length(node.noise_terms) > 1
	error(
	"RegularizedForwardPass cannot be applied because first-stage is " *
	"not deterministic",
	)
	end
	old_bounds = Dict{Symbol,Tuple{Float64,Float64}}()
	for (k, v) in node.states
	if has_lower_bound(v.out) && has_upper_bound(v.out)
	old_bounds[k] = (l, u) = (lower_bound(v.out), upper_bound(v.out))
	x = get(fp.trial_centre, k, model.initial_root_state[k])
	set_lower_bound(v.out, max(l, x - fp.ρ * (u - l)))
	set_upper_bound(v.out, min(u, x + fp.ρ * (u - l)))
	end
	end
	pass = forward_pass(model, options, fp.forward_pass)
	for (k, (l, u)) in old_bounds
	fp.trial_centre[k] = pass.sampled_states[1][k]
	set_lower_bound(node.states[k].out, l)
	set_upper_bound(node.states[k].out, u)
	end
	return pass
	end

[FSchmidtDIW]

Hi Oscar, thanks for checking! I had started implementing a new forward pass but got stuck.
The method requires computing a cut intercept adjustment, cf:

The 'norm' part $\sqrt{(1+||c_t + \beta_t^j||_2^2)}$ remains constant for a given but whereas the $\bar\sigma^k_t$ changes between iterations. So if a new cut has been added in the last iteration, we need to compute the norm part once.
So for every node and every new cut in that node we need to retrieve the cost vector (associated with state variables I assume?) and the cut coefficient vector to compute the 'norm' part. At least for the simple case of a linear graph, this should be doable along the lines of

node = model[node_index] # do I need to loop over all nodes in the forward_pass function?

objective_dict = objective_function(node.subproblem).terms
beta = Dict{Int64,Vector{Float64}}()
nterm[j] = Dict{Int64,Float64}()
c = [ haskey(objective_dict,i) ? objective_dict[i] : 0.0 for i in values(node.states)]
for j in new_cut_index_set # not sure how to define this set most efficiently
  beta[j] =  values(node.bellman_function.global_theta.cuts[j].coefficients)
  nterm[j] = sqrt(1+sum((beta[j]+c).^2)
end 

I thought it might be best to define $\bar\sigma_t^k$ as a variable and add it to the cut expression. This way, one can fix it to changing values. The nterm becomes simply a coefficient of that variable. We would not have to add and subtract different values from the intercept as $\sigma$ changes. Does that make sense? Where would be the best place to do that?

For the update of $\bar\sigma_t^k$ the authors suggest using the last lower bound improvement relative to an average of prior improvements, see below ($k_CC$ is a scaling factor). Is there an object where the lower bound trajectory is recorded?

If this approach is not completely off-track, I'd be happy to continue with the implementation.

Best

F

[odow]

I strongly dislike any methods which require modifying existing cut intercepts. It is not something that is easily achievable in SDDP.jl. I don't know if you should spend time working on this.

Why do you even need to modify the cuts? Can't you just perturb the states in the forward pass to be closer to the centre of the feasible region polyhedron?

I didn't read their paper, but modifying the cuts is not what I expected to happen.

[FSchmidtDIW]

I understand! That makes sense! I'll close this issue! I suppose if you pertubed the states towards the centre you'd end up with a stabilisation/gravity centre-like approach where the stability centre is the (Chebyshev) centre of the feasible region. May be worth a shot!

[odow]

If you're thinking about this:

• The forward pass can be computed independently of the current policy (the cuts)
• You don't even need the sequence of states in the forward pass to be consistent. You don't need x_2 to be reachable from x_1.
• I could even sample a forward pass with rand(Uniform(lower_bound(x), upper_bound(x))

One approach would be to add some constraint that restricts state variables to some trust region around the Chebyshev centre, solve the forward pass, and then delete the constraint in preparation for the backwards pass. This is kinda what this does (for the first-stage only):

SDDP.jl/src/plugins/forward_passes.jl

Lines 356 to 370 in 3b2570b

	old_bounds = Dict{Symbol,Tuple{Float64,Float64}}()
	for (k, v) in node.states
	if has_lower_bound(v.out) && has_upper_bound(v.out)
	old_bounds[k] = (l, u) = (lower_bound(v.out), upper_bound(v.out))
	x = get(fp.trial_centre, k, model.initial_root_state[k])
	set_lower_bound(v.out, max(l, x - fp.ρ * (u - l)))
	set_upper_bound(v.out, min(u, x + fp.ρ * (u - l)))
	end
	end
	pass = forward_pass(model, options, fp.forward_pass)
	for (k, (l, u)) in old_bounds
	fp.trial_centre[k] = pass.sampled_states[1][k]
	set_lower_bound(node.states[k].out, l)
	set_upper_bound(node.states[k].out, u)
	end

A second approach would be to solve the forward pass as usual, and then take a convex combination of the forward pass and the Chebyshev center a * x + (1 - a) * chebyshev as the "forward pass" that we computed. (This is what I thought the paper was probably doing).

Then you can make a go to zero as the iterations progress.

Some other things to consider:

• It's not obvious that we even want to do this.
• Regularizing the forward pass to "nice" trajectories means that we add fewer cuts to extreme parts of the state space
• In theory, this improves the in-sample approximation of the policy
• It's not obvious what it does to the out-of-sample approximation of the policy
• The goal of SDDP isn't to improve the lower bound as quickly as possible. It's to find a good enough out-of-sample heuristic policy.
• Finding the "optimal" policy for a highly approximated model isn't necessarily better than an approximate policy for a more detailed model.

[FSchmidtDIW]

It's not obvious that we even want to do this.
Regularizing the forward pass to "nice" trajectories means that we add fewer cuts to extreme parts of the state space
In theory, this improves the in-sample approximation of the policy
It's not obvious what it does to the out-of-sample approximation of the policy

That's actually a very fair point. Working on capacity expansion models, it appears useful for the first stage at least but I suppose the existing RegularizedForwardPass already does that.

Additionally, as the feasible polyhedron in a given node changes between iterations one would have to recompute the centre every time, which is probably not a major computational effort but complicates things. I will probably drop it for my current project and may come back to it after.

[odow]

Working on capacity expansion models, it appears useful for the first stage at least but I suppose the existing RegularizedForwardPass already does that

Yes, now thiiiiiis makes sense. Feel free to open a PR if you have ideas for how we could improve RegularizedForwardPass.

[FSchmidtDIW]

Will do! Thank you! I recall that at some point before adding RegularizedForwardPass you had tried out different bundle methods like level set / quadratic regularisation /trust region, right?
If not, I might give them a go. The real problem, I suppose, is that most hyperparameters in these approaches depend on the upper bound in the classic two-stage setup.

[odow]

That was for the Lagrangian dual problem. I didn't experiment much with RegularizedForwardPass. Any and all ideas accepted 😄

[FSchmidtDIW]

on it :)