JuMP or Convex: stuck with very specific constraint, not sure if possible as originally intended

Specific Domains Optimization (Mathematical)
,
1

Hi everyone,

Let me start with a “visual” of my issue: imagine a binary matrix with some constraints on how the rows and columns can “interact” with each other. In my case, I need a constraint within each row that guarantees if any values are set to 1, exactly n of them are, and they are all consecutive indices (similar to an SOS2 constraint, though not exactly, more like SOSn in lp_solve).

E.g.
n=3
Matrix of binary decision vars:
[0 0 0 1 0 1 1 0 0 0] (not valid)
[0 0 0 1 1 1 0 0 0 0] (valid)
[0 0 0 0 0 0 0 0 0 0] (valid)
[0 1 1 0 0 0 0 0 0 0] (not valid)

All rows must be valid (other constraints also imposed, but not necessary to understand this sub problem).

I thought I found a really elegant representation of this continuity constraint using a convolution with a box blur kernel vector of length n. Apply this to each row, verify the output is an integer (0 or 1), and we’re done. However, I have encountered a few roadblocks that have been quite vexing, and I’m not sure I can still leverage my original formulation.

  1. In JuMP, I can’t use convolutions from DSP or Convex or the maximum function at all, even though these are all affine transformations according to Convex.

  2. In Convex, the convolution operator doesn’t work between a vector kernel and a matrix variable with the following error: ERROR: convolution only supported between two vectors (though this seems to contradict the documentation with the vector being in R^n and matrix in R^m being permissible). Maybe this could be fixed by using multiple variables, one for each row, but this causes other issues as the number of rows can vary problem to problem.

  3. I need to verify the maximum result of the convolution is an integer (either 0 or 1 specifically), which I intended to do with the modulus operation, or by saying the interval (0,1) is excluded, but I don’t think this is supported by Convex, but maybe it is by JuMP?

It feels like I’m super close, but this constraint is just barely out of my reach currently. Any help would be greatly appreciated. I’d post code, but what I have currently likely isn’t of much help, but I will do so upon request.

1 Like
2

Can you provide a reproducible example of the things you have currently tried? In general, you can’t just use functions from other packages in JuMP constraints, you’d have to formulate it as a MILP.
One approach would be something like this:

using JuMP
num_variables, N = 10, 3
num_patterns = num_variables - N
patterns = zeros(Int, num_variables, num_patterns)
for i in 1:num_patterns
    patterns[i:(i+N), i] .= 1
end
model = Model()
@variable(model, x[1:num_variables], Bin)
@variable(model, z[1:num_patterns], Bin)
@constraint(model, sum(z) <= 1)
@constraint(model, [i=1:num_variables], x[i] == sum(patterns[i, :]' * z))
2 Likes
3

Hi John,

Building off Oscar’s answer, here is an approach that wraps it in a function. I’m using the standard library module SparseArrays for the times where your continuity parameter n (or N) is much smaller than your overall matrix row length. The other adjustment is to change N to N-1 in Oscar’s answer to match up with your example; e.g. if N = 3, you have 3 consecutive ones in a rows.

using SparseArrays
using JuMP

function create_continuity_constraint(model::JuMP.Model, num_variables::T, N::T) where {T<:Integer}
    num_patterns = num_variables - (N - 1)
    patterns = spzeros(Int, num_variables, num_patterns) # allocate a sparse array
    for i in 1:num_patterns
        patterns[i:(i+(N-1)), i] .= 1
    end

    @variable(model, x[1:num_variables], Bin)
    @variable(model, z[1:num_patterns], Bin)
    @constraint(model, sum(z) <= 1)
    @constraint(model, [i=1:num_variables], x[i] == sum(patterns[i, :]' * z))
    return nothing
end

# Test it out:
model = Model()
num_variables, N = 10, 3
create_continuity_constraint(model, num_variables, N)
print(model)
1 Like
4

Thank you to both @odow and @jd-foster ! Both of these answers were a huge help! I did catch that N-1 fix, but I think I found a more compact representation when I was finagling with it.
Your help got me super close to my final set of constraints/variables; this is what I ended up starting with (and I’ll probably adjust and use that function/SparseArray tip):

num_variables, N = 10, 3
num_patterns = num_variables - (N-1)
patterns = zeros(Int, num_variables, num_patterns)
for i in 1:num_patterns
    patterns[i:(i+(N-1)), i] .= 1
end
model = Model()
@variable(model, z[1:num_patterns], Bin)
@constraint(model, [i=1:num_variables], (patterns*z)[i] <= 1)

The reason this last constraint works is that the product of the patterns matrix with our binary decision variable takes the sum of the product of the rows. As long as that sum is <= 1, there’s no overlap for that position.

For my own edification, can I ask why x, the sum(z) <=1 and the last constraint were all formulated as such? I really didn’t understand the intention there.

5

sum(z) <= 1 says “choose at most one pattern” if we choose zero patterns, then that correspond to the [0, 0, .., 0] pattern.

The x recovers your original vector of interest.

@constraint(model, [i=1:num_variables], (patterns*z)[i] <= 1)

This is better written as @constraint(model, patterns * z .<= 1), but you’re much better just doing sum(z) <= 1. If you don’t want to introduce the explicit x variable, you can use

@expression(model, x, patterns * z)
6

Oh I see; in my case, I actually needed to be able to choose multiple patterns, and the patterns can’t overlap; hence sum(z) can actually be >1, but patterns*z yields the product of the possible “shifts” and if that only contains 1s and 0s then I know there are no overlaps.

Thank you for the easier expression with the broadcast notation; I’m still getting used to JuMP so I’m not completely sure what’s allowed and what isn’t :slight_smile: