With regards to the scale of the problem, let us say N\leq{24}. For larger values, I am only interested in finding “a” solution and not necessarily the optimal answer.
The elements of M are indeed structured. Each row corresponds to a binary representation of the numbers 0\cdots{2}^N-1. Each column is the AND of two bits of the binary representation, see my code below. This code has a few additional details that were omitted in the previous description.
With this code, the solver starts to struggle even with N=10. FYI, I have written a python code to perform the search directly and so far have found some solutions even for N=20 (I don’t know about their optimality though).
I think, a promising direction to improve the speed would be utilizing the fact that for N=10, T_cnt_max \leq 100. I should probably try to calculate the cardinality of the 2^N-element columns in multiple steps instead of doing it in a single shot i.e. break the columns into smaller portions, determine the cardinality on each chunk and then combine the results. Doing it that way the number of variables/conditions for T_eql_m does not scale by 2^{2N} (it would be something like 100\times{2^N} which is still better).
Any other suggestions ? Is there anything else that I can improve in the code below ?
option solver './ampl/cplex'; option cplex_options 'mipdisplay=2 populate=1';
param N > 0 integer;
param M_max = 2^(N-3);
set diag := 1..N-1;
set rows := 0..2^N-1;
set cols := 0..N*(N-1)/2-1;
set rnge := 0..M_max;
param M {rows,cols} >= 0 integer;
# How to define SOS1, does it speed up the solver ?
var E {c in cols,d in diag} binary; subject to E_cond {c in cols} : ( sum {d in diag} E[c,d] ) == 1;
# Partitioning M columns based on E into N-1 subsets, the columns in each subset are added together to create T
var T {r in rows,d in diag} = sum {c in cols} M[r,c]*E[c,d];
# T_mod = T mod 2^(N-2) (largest value in T is 2^(N-1))
var T_div {r in rows,d in diag} integer >= 0, <= 3;
var T_mod {r in rows,d in diag} = T[r,d]-T_div[r,d]*2^(N-2); subject to T_mod_cond {r in rows,d in diag}: 0 <= T_mod[r,d] <= 2^(N-2)-1;
# T_adj = min(T_mod,2^(N-2)-T_mod)
var T_sel {r in rows,d in diag} binary;
var T_adj {r in rows,d in diag} integer; subject to T_adj_cond_1 {r in rows,d in diag} : T_mod[r,d]-2^(N-2)* T_sel[r,d] <= T_adj[r,d];
subject to T_adj_cond_2 {r in rows,d in diag} : T_adj[r,d] <= T_mod[r,d];
subject to T_adj_cond_3 {r in rows,d in diag} : 2^(N-2)-T_mod[r,d]-2^(N-2)*(1-T_sel[r,d]) <= T_adj[r,d];
subject to T_adj_cond_4 {r in rows,d in diag} : T_adj[r,d] <= 2^(N-2)-T_mod[r,d];
# T_eql_m = T_adj >= m && T_adj <= m (this is the part that should be improved)
var T_mte_m {r in rows,d in diag,m in rnge} binary; subject to T_mte_m_cond {r in rows,d in diag,m in rnge}: 0<= T_adj[r,d]-m+(2*M_max+1)*(1-T_mte_m[r,d,m]) <= 2*M_max;
var T_lte_m {r in rows,d in diag,m in rnge} binary; subject to T_lte_m_cond {r in rows,d in diag,m in rnge}: 0<= m-T_adj[r,d]+(2*M_max+1)*(1-T_lte_m[r,d,m]) <= 2*M_max;
var T_eql_m {r in rows,d in diag,m in rnge} binary; subject to T_eql_m_cond {r in rows,d in diag,m in rnge}: 0<= T_mte_m[r,d,m]+T_lte_m[r,d,m]-2*T_eql_m[r,d,m] <= 1;
# T_eql_m_all = T_eql_m(r=0) || T_eql_m(r=1) || ... || T_eql_m(r=2^N-1)
var T_eql_m_all {d in diag,m in rnge} binary; subject to T_eql_m_all_cond {d in diag,m in rnge}: 0<= 2^N*T_eql_m_all[d,m] - ( sum {r in rows} T_eql_m[r,d,m] ) <= 2^N-1;
# T_cnt is the number of different elements in each column (cardinality of each column)
var T_cnt {d in diag} = sum {m in 0..M_max} T_eql_m_all[d,m];
# Different objectives
var T_cnt_max integer >= 0; subject to T_cnt_max_cond {d in diag}: T_cnt[d] <= T_cnt_max;
var T_cnt_sum = sum {d in diag} T_cnt[d];
minimize cost: T_cnt_max*1000+T_cnt_sum;
problem find_sol: E, E_cond, T, T_div, T_mod, T_mod_cond, T_sel, T_adj, T_adj_cond_1, T_adj_cond_2, T_adj_cond_3, T_adj_cond_4,
T_mte_m, T_mte_m_cond, T_lte_m, T_lte_m_cond, T_eql_m, T_eql_m_cond, T_eql_m_all, T_eql_m_all_cond,
T_cnt, T_cnt_max, T_cnt_max_cond, T_cnt_sum, cost;
let N := 8;
check;
param r;
param B {0..N-1};
for {k in 0..2^N-1}
{
let r := 0;
for {b in 0..N-1}
let B[b] := ( k div 2^b ) mod 2;
for {j in 1..N-1}
{
for {i in 1..N-j}
{
let M[k,r] := 2^(N-1-j)*B[N-i]*B[N-i-j];
fix {d in diag: d <= i-1} E[r,d] := 0;
fix {d in diag: d >= i+j} E[r,d] := 0;
let r := r + 1;
}
}
}
fix {c in cols,d in diag: c >= N-1 and (d mod 2) == 1} E[c,d] := 0;
solve find_sol;
#display T;
#display T_div;
#display T_mod;
#display T_adj;
display cost, T_cnt_max, T_cnt_sum, T_cnt;
quit