Compact form of a set of constraints

rocco_sprmnt21 · May 21, 2025, 1:43pm

I have the following script that works for the case N=4.

using JuMP 
import HiGHS
m = Model(HiGHS.Optimizer);
@variable(m, x[1:4,1:8], Bin)
@variable(m, y[1:4,1:6], Bin)

@constraints(m, begin
    [j in 1:8], sum(x[r,j] for r in 1:4)==1

    [i in 1:3],   y[4,i]-->{x[4,i]+x[4,i+5]==2}
       sum(y[4,j] for j in 1:3)==1
    [i in 1:4],   y[3,i]-->{x[3,i]+x[3,i+4]==2}
       sum(y[3,j] for j in 1:4)==1
    [i in 1:5],   y[2,i]-->{x[2,i]+x[2,i+3]==2}
       sum(y[2,j] for j in 1:5)==1
    [i in 1:6],   y[1,i]-->{x[1,i]+x[1,i+2]==2}
       sum(y[1,j] for j in 1:6)==1
    end)
optimize!(m)
julia> value.(x)
4×8 Matrix{Float64}:
 0.0  1.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0
 0.0  0.0  1.0  0.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0

Can I get a compact form that allows me to do tests for higher values of N, without having to hand write all the constraints, case by case?

Is something like this, with the appropriate and correct syntax possible?


N=4
m = Model(HiGHS.Optimizer);
@variable(m, x[1:N,1:2N], Bin)
@variable(m, y[1:N,1:2(N-1)], Bin)

@constraints(m, begin
    [j in 1:2N], sum(x[r,j] for r in 1:N)==1

    for k in N-1:-1:1
    [i in 1:k],   y[k+1,i]-->{x[k1,i]+x[k+1,i+k+2]==2}
       sum(y[k+1,j] for j in 1:k)==1
    end


    end)
optimize!(m)

value.(x)

rocco_sprmnt21 · May 21, 2025, 2:01pm

maybe like this?

@constraints(m, begin
    [j in 1:2N], sum(x[r,j] for r in 1:N)==1
    [k in N-1:-1:1, i in 1:k],   y[k+1,i]-->{x[k+1,i]+x[k+1,i+k+2]==2}
    [k in N-1:-1:1],    sum(y[k+1,j] for j in 1:k)==1

    end)
optimize!(m)

odow · May 21, 2025, 9:29pm

Your second comment looks like it’s on the right track, but you might want to double check your logic.

I’d probably write your model as:

using JuMP 
import HiGHS
function main(N)
    model = Model(HiGHS.Optimizer)
    @variable(model, x[1:N, 1:2N], Bin)
    @variable(model, y[1:N, 1:2(N-1)], Bin)
    @constraints(model, begin
        [j in 1:2N], sum(x[:,j]) == 1
        [i in 1:N], sum(y[i,1:(2N-1-i)]) == 1
        [k in 1:N, i in 1:(2N-1-k)], y[k,i] --> {x[k,i] + x[k,i+k+1] == 2}
    end)
    optimize!(model)
    assert_is_solved_and_feasible(model)
    return value.(x), value.(y)
end

rocco_sprmnt21 · May 22, 2025, 4:54am

Yes, of course. That was to confirm that the syntax was valid.
I put it that way.
Why does the solution x (declared Bin) have real values and not integers?

m = Model(HiGHS.Optimizer);
@variable(m, x[1:N,1:2N], Bin)
@variable(m, y[1:N,1:2(N-1)], Bin)

@constraints(m, begin
   [k in 1:N, i in 1:2N-k-1],   y[k,i]-->{x[k,i]+x[k,i+k+1]==2}
   [k in 1:N],    sum(y[k,i] for i in 1:2(N-1)-k+1)==1
   [j in 1:2N], sum(x[r,j] for r in 1:N)==1
    end)
optimize!(m)

round.(Int,value.(x))

What role does this expression have?

assert_is_solved_and_feasible(model)

WalterMadelim · May 22, 2025, 5:14am

Because Optimizers work with float64.

A normality check. See the docstring.

help?> JuMP.is_solved_and_feasible
  is_solved_and_feasible(
      model::GenericModel;
      allow_local::Bool = true,
      allow_almost::Bool = false,
      dual::Bool = false,
      result::Int = 1,
  )

.....

    •  the primal_status of the result index result is FEASIBLE_POINT.

  This function is conservative, in that it returns false for situations like the solver terminating with a feasible solution due to a time      
  limit.

  If this function returns false, use termination_status, result_count, primal_status and dual_status to understand what solutions are
  available (if any).

  See also: assert_is_solved_and_feasible.

rocco_sprmnt21 · May 22, 2025, 5:22am

in my system it is not defined (I probably need to update)

ERROR: UndefVarError: `assert_is_solved_and_feasible` not defined

I provide some more information on the problem.
I was trying to understand for which values of N it is possible to find a sequence of 2N numbers, containing all the pairs (1,1) (2,2) …(N,N) positioned at a distance equal to the number. For example (1,x,1), (3,x,x,x,3) …
in the case N=11

[1:11...]'*main(11)
...
 6  8  10  11  1  2  1  6  2  7  8  9  5  10  4  11  3  7  5  4  3  9

odow · May 22, 2025, 5:33am

Take a read of Tolerances and numerical issues · JuMP

And yes, you must have an old version of JuMP. It should work after updating.

rocco_sprmnt21 · September 21, 2025, 1:37pm

I tried to rewrite the set of constraints differently, but equivalently.
The following attempt shows that the equivalence is not true.
But I don’t know how to interpret the second rows of the y and x matrices.
x[2,14]=1 => y[2,14]+y[2,17]==2
An error should be reported or this configuration should be considered outside the region defined by the constraints.
What is the reality in this case?

function main(N)
    model = Model(HiGHS.Optimizer)
    set_silent(model)
    @variable(model, x[1:N, 1:2N], Bin)
    @variable(model, y[1:N, 1:2(N-1)], Bin)
    @constraints(model, begin
        [j in 1:2N], sum(x[:,j]) == 1
        [j in 1:2N-2], sum(y[:,j]) <= 1
        [i in 1:N], sum(x[i,:]) == 2
#[i in 1:N], sum(y[i,1:(2N-1-i)]) == 1
        [i in 1:N], sum(y[i,:]) == 1
        [k in 1:N, i in 1:(2N-1-k)], y[k,i] --> {x[k,i] + x[k,i+k+1] == 2}
    end)
    optimize!(model)
    assert_is_solved_and_feasible(model)
    #return value.(x), value.(y)
    round.(Int,value.(x)), round.(Int,value.(y))
end
# l'insieme di vincoli di questa versione è APPARENTEMENTE equivalente alla versione precedente
# ma il risultato che segue è un controesempio di questa congettura
#
# julia> x
# 8×16 Matrix{Int64}:
#  0  0  0  0  0  0  0  0  0  1  0  1  0  0  0  0
#  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  0
#  0  0  1  0  0  1  0  0  0  0  0  0  0  0  0  0
#  1  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
#  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  0
#  0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0
#  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1
#  0  0  0  1  0  0  0  0  0  0  1  0  0  0  0  0

# julia> y
# 8×14 Matrix{Int64}:
#  0  0  0  0  0  0  0  0  0  1  0  0  0  0
#  0  0  0  0  0  0  0  0  0  0  0  0  0  1
#  0  0  0  0  0  0  0  0  0  0  0  0  1  0
#  0  0  0  0  0  0  0  0  0  0  0  1  0  0
#  0  0  0  0  0  0  0  0  0  0  1  0  0  0
#  0  1  0  0  0  0  0  0  0  0  0  0  0  0
#  0  0  0  0  0  0  0  0  1  0  0  0  0  0
#  0  0  0  0  0  0  0  1  0  0  0  0  0  0

PS

How is the implication y → {x1+x2 ==2} translated into terms of linear constraints?

rocco_sprmnt21 · September 21, 2025, 2:08pm

I’ll answer myself, maybe.
The constraint y[2,14+]–> {x + x} is not actually defined.
So there is no implication

Adding this constraint
[i in 1:N], sum(y[i,2N-i:2N-2]) == 0
seems to work.

           model = Model(HiGHS.Optimizer)
           set_silent(model)
           @variable(model, x[1:N, 1:2N], Bin)
           @variable(model, y[1:N, 1:2(N-1)], Bin)
           @constraints(model, begin
               [j in 1:2N], sum(x[:,j]) == 1
               [j in 1:2N-2], sum(y[:,j]) <= 1
               [i in 1:N], sum(x[i,:]) == 2
       #[i in 1:N], sum(y[i,1:(2N-1-i)]) == 1
               [i in 1:N], sum(y[i,:]) == 1
               [i in 1:N], sum(y[i,2N-i:2N-2]) == 0
               [k in 1:N, i in 1:(2N-1-k)], y[k,i] --> {x[k,i] + x[k,i+k+1] == 2} 
           end)
           optimize!(model)
           assert_is_solved_and_feasible(model)
           #return value.(x), value.(y)
           round.(Int,value.(x)), round.(Int,value.(y))
       end
main (generic function with 1 method)

julia> 

julia> 

julia> x,y = main(N)
([0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 1 0], [0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0])

julia> 

julia> 

julia> 

julia> x
8×16 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0  0  1  0  1  0  0
 0  0  0  1  0  0  1  0  0  0  0  0  0  0  0  0
 1  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  0
 0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  1
 0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0  0  1  0  0  0  0  0
 0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  0

julia> y
8×14 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0  0  1  0  0
 0  0  0  1  0  0  0  0  0  0  0  0  0  0
 1  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  1  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  1  0  0  0  0  0  0  0  0

odow · September 21, 2025, 9:31pm

Is there still a question? Or did you figure it out?

If N = 8 and k = 2, then i = 1:13, so yeah, there is no implication for y[2,14].

How is the implication y → {x1+x2 ==2} translated into terms of linear constraints?

It depends on the solver. But in this specific case:

x1 + x2 + s == 2
s >= 0
s + 2y <= 2

If y = 1, then s = 0. If y = 0, then 0 <= s <= 2, and so the == 2 constraint is satisfiable for any combination of x1 and x2.

Topic		Replies	Views
Help with a Constraint Optimization (Mathematical)	15	3251	November 1, 2018
For syntax in @constraints Optimization (Mathematical) question , jump	1	782	October 30, 2019
How to define a constraint from a string? Optimization (Mathematical)	12	1920	July 10, 2018
How to efficiently generate constraints Optimization (Mathematical)	4	914	January 10, 2019
Trying to solve a puzzle Optimization (Mathematical)	23	1213	May 23, 2021

Compact form of a set of constraints

Related topics