Provided the value function is still convex w.r.t. the state variables, and we can obtain valid subgradients, there is nothing to say that the controls have to be convex.
In a previous iteration of SDDP.jl, we implemented this by solving the MIP's on the forward pass, then on the backwards pass, we solve the MIP, fix the integer variables and then solve the LP relaxation. This doesn't necessarily provide valid subgradients but it will create a policy that might be good enough.
This is particularly useful if you want to do things like use SOSII to model things.
Leaving some code that may be potentially useful.
nonzero(x) = abs(x) > 1e-10
function fixvariable!(m::Model, idx, value)
setupperbound(Variable(m, idx), value)
setlowerbound(Variable(m, idx), value)
end
function fixmodel!(m::Model)
!m.internalModelLoaded && error("Unable to fix an unsolved model")
# Follow what Gurobi does for Special Ordered Sets
if length(m.sosconstr) > 0
for sos in m.sosconstr
cols = Int[v.col for v in sos.terms]
if sos.sostype == :SOS1
for c in cols
nonzero(m.colVal[c]) || fixvariable!(m, c, 0.)
end
elseif sos.sostype == :SOS2
colorder = sortperm(sos.weights)
for (i, sosidx) in enumerate(colorder)
c = cols[sosidx]
nonzero(m.colVal[c]) || fixvariable!(m, c, 0.)
end
else
error("Unable to create fixed model with SOS of type $(sos.sostype)")
end
end
m.sosconstr = JuMP.SOSConstraint[]
# Do this for now since no efficient way to update/remove SOS constraints
m.internalModelLoaded = false
end
# Fix all binary and integer variables to their current value
for (i, v) in enumerate(m.colCat)
if v == :Bin || v == :Int
setcategory(Variable(m, i), :Cont)
fixvariable!(m, i, m.colVal[i])
end
end
end
function solvenoncontinuous!(sp::Model, scenario::Int)
lb = copy(sp.colLower)
ub = copy(sp.colUpper)
cat = copy(sp.colCat)
sos = copy(sp.sosconstr)
fixmodel!(sp)
status = solve(sp)
@assert status == :Optimal
storecontinuous!(sp, scenario)
sp.colLower = lb
sp.colUpper = ub
sp.colCat = cat
if length(sos) > 0
sp.sosconstr = sos
sp.internalModelLoaded = false
end
end
function SOSII!(m::JuMP.Model, f::Function, x, lb::Float64, ub::Float64, n::Int64=10)
# This function adds a SOS of Type II to the JuMP model
# Σλ = 1 # Convexity
# a∙λ = z # x value
# f(z) = f(a)∙λ # approximated function value
# where λ is a SOS of Type II
xx = linspace(lb, ub, n) # Set x value
fx = f(xx) # Estimate function value
SOSII!(m, x, xx, fx)
end
function SOSII!(m::JuMP.Model, x, xx::AbstractVector, fx::AbstractVector)
n = length(xx)
@assert length(fx) == n
y = @variable(m) # Estimated function value variable
λ = @variable(m, [1:n], lowerbound=0, upperbound=1) # SOS variable set
@constraints(m, begin
sum(λ) == 1 # Convexity
x == dot(xx, λ)
y == dot(fx, λ)
end)
addSOS2(m, [i*λ[i] for i =1:n])
y # Return estimated function value variable
end
SOSII!{T<:Real}(m::JuMP.Model, x, xx::Vector{NTuple{2, T}}) = SOSII!(m, x, [y[1] for y in xx], [y[2] for y in xx])
function bilinear(x::Variable, y::Variable, n::Int=10)
# This macro is used to replace a bilinear term in a JuMP constraint with the expansion
# xy = 0.25 * [(x+y)^2 - (x-y)^2]
# where (x+y)^2 and (x-y)^2 are approximated with SOS2
m = x.m
# (x + y)^2 SOS2
v1 = SOSII!(m, (a) -> a.^2, x + y, getlowerbound(x) + getlowerbound(y), getupperbound(x) + getupperbound(y), n)
# (x + y)^2 SOS2
v2 = SOSII!(m, (a) -> a.^2, x - y, getlowerbound(x) - getupperbound(y), getupperbound(x) - getlowerbound(y), n)
# Return expression to be placed in constraint
(0.25 * (v1 - v2))
end