Let $A, B, C, D$ be $4 \times 4$ matrices over $\mathbb{Z}\_2$.

Suppose they satisfy
\begin{equation}
 \begin{cases}
 A^2+ BC+ BCA+ ABC+ A= I\_4\\
 AB+ BCB+ ABD=0\_{4 \times 4} \\
 CA+ DCA+ CBC= 0\_{4 \times 4} \\
 DCB+ CBD= I\_4
 \end{cases}\,.
\end{equation}

Does anyone have an idea on how to find matrices $A, B, C, D$ either theoretically or numerically? I'm thinking that we can find a solution since there are four equations and four unknowns.

Hints would suffice. Thank you so much.

Edit: I found out that the 148 solutions is not yet the total number of complete solutions. To possibly decrease the number of solutions I added another constraint:

$$\sum\_{k=1}^4\sum\_{l=1}^4 A\_{1,k}A\_{k,l}A\_{l,1}+ \sum\_{k=1}^4\sum\_{l=1}^4 B\_{1,k}C\_{k,l}A\_{l,1}+ \sum\_{k=1}^4\sum\_{l=1}^4 A\_{1,k}B\_{k,l}C\_{l,1}+ \sum\_{k=1}^4\sum\_{l=1}^4 B\_{1,k}D\_{k,l}C\_{l,1}=2Y\_{5,1,1}.$$

I added this constraint to the last part of the code:

```
proc optmodel;
   num n = 4;
   set ROWS = 1..n;
   set COLS = ROWS;
   set CELLS = ROWS cross COLS;

   var A {CELLS} binary;
   var B {CELLS} binary;
   var C {CELLS} binary;
   var D {CELLS} binary;

   var Y {1..5, CELLS} >= 0 integer;

    /* M[i,j,k] = A[i,k]*A[k,j] */
   var M {1..n, 1..n, 1..n} binary;
   con Mcon1 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[i,k];
   con Mcon2 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[k,j];
   con Mcon3 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] >= A[i,k] + A[k,j] - 1;

    /* X[i,j,k] = B[i,k]*C[k,j] */
   var X {1..n, 1..n, 1..n} binary;
   con Xcon1 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] <= B[i,k];
   con Xcon2 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] <= C[k,j];
   con Xcon3 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] >= B[i,k] + C[k,j] - 1;

     /* O[i,j,k,l] = B[i,k]*C[k,l]*A[l,j] */
   var O {1..n, 1..n, 1..n, 1..n} binary;
   con Ocon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= B[i,k];
   con Ocon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= C[k,l];
   con Ocon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= A[l,j];   
   con Ocon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] >= B[i,k]+C[k,l]+A[l,j] - 2;  

      /* P[i,j,k,l] = A[i,k]*B[k,l]*C[l,j] */
   var P {1..n, 1..n, 1..n, 1..n} binary;
   con Pcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= A[i,k];
   con Pcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= B[k,l];
   con Pcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= C[l,j];   
   con Pcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] >= A[i,k]+B[k,l]+C[l,j] - 2;

   /* A^2 + BC + BCA + ABC + A = I_4 */
   con Con1 {<i,j> in CELLS}:
      sum {k in 1..n} M[i,j,k] 
    + sum {k in 1..n} X[i,j,k] 
    + sum {k in 1..n, l in 1..n} O[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} P[i,j,k,l]
    + A[i,j] 
   = 2*Y[1,i,j] + (i=j);

   /* E[i,j,k] = A[i,k]*B[k,j] */
   var E {1..n, 1..n, 1..n} binary;
   con Econ1 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] <= A[i,k];
   con Econ2 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] <= B[k,j];
   con Econ3 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] >= A[i,k] + B[k,j] - 1;

    /* F[i,j,k,l] = B[i,k]*C[k,l]*B[l,j] */
   var F {1..n, 1..n, 1..n, 1..n} binary;
   con Fcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= B[i,k];
   con Fcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= C[k,l];
   con Fcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= B[l,j];   
   con Fcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] >= B[i,k]+C[k,l]+B[l,j] - 2;

      /* G[i,j,k,l] = A[i,k]*B[k,l]*D[l,j] */
   var G {1..n, 1..n, 1..n, 1..n} binary;
   con Gcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= A[i,k];
   con Gcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= B[k,l];
   con Gcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= D[l,j];   
   con Gcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] >= A[i,k]+B[k,l]+D[l,j] - 2; 

   /* AB + BCB + ABD = 0 */
   con Con2 {<i,j> in CELLS}:
      sum {k in 1..n} E[i,j,k] 
    + sum {k in 1..n, l in 1..n} F[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} G[i,j,k,l]
   = 2*Y[2,i,j];

   /* H[i,j,k] = C[i,k]*A[k,j] */
   var H {1..n, 1..n, 1..n} binary;
   con Hcon1 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] <= C[i,k];
   con Hcon2 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] <= A[k,j];
   con Hcon3 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] >= C[i,k] + A[k,j] - 1;

    /* Q[i,j,k,l] = D[i,k]*C[k,l]*A[l,j] */
   var Q {1..n, 1..n, 1..n, 1..n} binary;
   con Qcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= D[i,k];
   con Qcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= C[k,l];
   con Qcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= A[l,j];   
   con Qcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] >= D[i,k]+C[k,l]+A[l,j] - 2; 

    /* R[i,j,k,l] = C[i,k]*B[k,l]*C[l,j] */
   var R {1..n, 1..n, 1..n, 1..n} binary;
   con Rcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= C[i,k];
   con Rcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= B[k,l];
   con Rcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= C[l,j];   
   con Rcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] >= C[i,k]+B[k,l]+C[l,j] - 2; 

   /* CA + DCA + CBC = 0 */
   con Con3 {<i,j> in CELLS}:
      sum {k in 1..n} H[i,j,k]
    + sum {k in 1..n, l in 1..n} Q[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} R[i,j,k,l]
   = 2*Y[3,i,j];

   /* W[i,j,k,l] = D[i,k]*C[k,l]*B[l,j] */
   var W {1..n, 1..n, 1..n, 1..n} binary;
   con Wcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= D[i,k];
   con Wcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= C[k,l];
   con Wcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= B[l,j];   
   con Wcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] >= D[i,k]+C[k,l]+B[l,j] - 2; 

      /* T[i,j,k,l] = C[i,k]*B[k,l]*D[l,j] */
   var T {1..n, 1..n, 1..n, 1..n} binary;
   con Tcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= C[i,k];
   con Tcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= B[k,l];
   con Tcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= D[l,j];   
   con Tcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] >= C[i,k]+B[k,l]+D[l,j] - 2;

   /* DCB + CBD = I_4 */
   con Con4 {<i,j> in CELLS}:
      sum {k in 1..n, l in 1..n} W[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} T[i,j,k,l]
   = 2*Y[4,i,j] + (i=j);

   /* U[1,1,k,l] = A[1,k]*A[k,l]*A[l,1] */
   var U {1..n, 1..n} binary;
   con Ucon1 {k in 1..n, l in 1..n}:
      U[k,l] <= A[1,k];
   con Ucon2 {k in 1..n, l in 1..n}:
      U[k,l] <= A[k,l];
   con Ucon3 {k in 1..n, l in 1..n}:
      U[k,l] <= A[l,1];   
   con Ucon4 {k in 1..n, l in 1..n}:
      U[k,l] >= A[1,k]+A[k,l]+A[l,1] - 2;

    /* V[1,1,k,l] = B[1,k]*D[k,l]*C[l,1] */
   var V {1..n, 1..n} binary;
   con Vcon1 {k in 1..n, l in 1..n}:
      V[k,l] <= B[1,k];
   con Vcon2 {k in 1..n, l in 1..n}:
      V[k,l] <= D[k,l];
   con Vcon3 {k in 1..n, l in 1..n}:
      V[k,l] <= C[l,1];   
   con Vcon4 {k in 1..n, l in 1..n}:
      V[k,l] >= B[1,k]+D[k,l]+C[l,1] - 2;  

   con 
   sum {k in 1..n, l in 1..n} U[k,l] 
    + sum {k in 1..n, l in 1..n} V[k,l]+ sum {k in 1..n, l in 1..n} O[1,1,k,l] + sum {k in 1..n, l in 1..n} P[1,1,k,l] 
   = 2*Y[5,1,1] ;

   solve noobj with milp / maxpoolsols=100 soltype=best;
   create data out(drop=s) from [s]=(1.._NSOL_) {i in 1..n, j in 1..n} <col('A_'||i||''||j)=A[i,j].sol[s]> {i in 1..n, j in 1..n} <col('B'||i||''||j)=B[i,j].sol[s]> {i in 1..n, j in 1..n} <col('C'||i||''||j)=C[i,j].sol[s]> {i in 1..n, j in 1..n} <col('D'||i||'_'||j)=D[i,j].sol[s]>;  

   print A;
   print B;
   print C;
   print D;
quit;

```

The software found 100 solutions, but when I checked the output data, it is the same solution replicated 100 times. Is there something wrong with the syntax of the constraint that I added to the code?