I’m trying to reformulate a stochastic program (approximately) as a mixed-integer linear program. A particular constraint involves the summation of N i.i.d. random variables, where N is a decision variable itself, something like

m = Model()
@variable(m, N >= 0, Int)
@variable(m, x >= 0)
@constraint(m, x <= sum(rand() for i=1:N))

The above formulation does not work. I also tried to replace N with getvalue(N) in the above constraint, but the solution doesn’t seem to make sense. Can anyone share some thoughts on this?

Since you’re trying to formulate a mixed-integer LP, the key question is: if you remove all the integrality constraints, do you end up with an LP? In this case, no, because the constraint doesn’t even make sense for real-valued N.

I don’t think you can model it the way you’ve written it, but there’s another formulation that should be equivalent. Rather than trying to use N as a decision variable, you can represent your summation by creating a new binary decision variable associated with each of your random variables. Then you can use those binary variables to encode whether each random variable is included in the summation:

using JuMP, Cbc
m = Model(solver=CbcSolver())
# The number of random variables
Nmax = 10
# The random variables themselves
x = randn(Nmax)
# The z variables select which x are included in the sum
@variable m z[1:Nmax] Bin
y = sum(x .* z)
# This constraint means that z must always take the form
# [1, 1, ..., 1, 0, 0, ...., 0]
# which is equivalent to saying y = sum(x[1:N]) for some N
for i in 2:Nmax
@constraint m z[i] <= z[i - 1]
end
@objective m Min y
solve(m)
@show x getvalue(y) getvalue(z)