基于线性规划寻找含最多公共1的三元矩阵组的技术求助

现有6个二进制矩阵如下：

matrix 1: [0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1]
matrix 2: [0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1]
matrix 3: [0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1]
matrix 4: [1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]
matrix 5: [1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1]
matrix 6: [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]

需求：从上述6个矩阵中选出3个组成一组，找到公共位置上1的数量最多的组。例如矩阵4、5、6的公共1数量为10，对应公共位置的1组成向量：

[1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1]

无法构建适用于Lingo等工具的求解模型，特此求助。

Lingo求解模型 #

变量定义

• 二进制变量 x(i)：x(i)=1 表示选中第i个矩阵（i=1~6），x(i)=0 表示不选
• 二进制变量 y(j)：y(j)=1 表示第j个位置是选中三矩阵的公共1（j=1~21）
• 常数 a(i,j)：第i个矩阵第j个位置的取值（0或1）

目标函数

最大化公共1的数量：

MAX = @SUM(POSITION(j): y(j));

约束条件

1. 恰好选中3个矩阵：

@SUM(MATRIX(i): x(i)) = 3;

2. 仅当选中的所有矩阵在第j位置均为1时，y(j)才能取1：

@FOR(POSITION(j):
    @SUM(MATRIX(i): (1 - a(i,j)) * x(i)) + y(j) <= 1;
);

逻辑说明：若某选中矩阵在j位置为0，则(1-a(i,j))*x(i)=1，此时y(j)必须为0；仅当所有选中矩阵在j位置都是1时，求和项为0，y(j)可取值1。
3. 变量类型约束：

@FOR(MATRIX(i): @BIN(x(i)));
@FOR(POSITION(j): @BIN(y(j)));

完整Lingo代码（代入矩阵数据）

MODEL:
SETS:
    MATRIX /1..6/;
    POSITION /1..21/;
    LINK(MATRIX, POSITION): a;
ENDSETS

DATA:
    a = 0 1 0 1 1 0 1 0 1 0 1 1 1 0 1 0 0 0 1 0 1
        0 1 0 1 1 0 0 1 0 1 0 1 0 0 0 1 1 1 0 1 1
        0 0 0 1 1 0 0 1 0 0 0 1 0 0 1 0 1 1 0 0 1
        1 1 0 1 1 1 0 1 0 0 0 1 1 1 0 0 0 0 0 1 1
        1 1 0 1 1 1 0 0 0 0 0 1 1 1 0 0 0 1 0 1 1
        1 1 1 1 1 1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1;
ENDDATA

MAX = @SUM(POSITION(j): y(j));

@SUM(MATRIX(i): x(i)) = 3;

@FOR(POSITION(j):
    @SUM(MATRIX(i): (1 - a(i,j)) * x(i)) + y(j) <= 1;
);

@FOR(MATRIX(i): @BIN(x(i)));
@FOR(POSITION(j): @BIN(y(j)));

END

内容的提问来源于stack exchange，提问作者Hrruuska