Using an integer-valued variable as the number of summands in JuMP?

Specific Domains Optimization (Mathematical)
1

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?

2

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.

1 Like
3

Good point. It may not be an linear program. But can I still model this type of constraint in JuMP?

4

Your application is not clear, but you cannot use a decision variable to specify the number of terms in a constraint!

If your goal from sum(rand() for i=1:N) is to get a random number from 1 to N, you can do the following:

@constraint(m, x <= rand()*N)
5

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)

which gives:

x = [0.557752, -0.31448, 0.721013, -0.141906, -1.37589, -0.343218, 0.369579, 2.0714, 0.674477, -0.953159]
getvalue(y) = -0.8967336932327254
getvalue(z) = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]
3 Likes
6

Thanks a lot!