reset;
param n := 5; # square size
param N := n^2; # numbers
param C := n*(n^2+1)/2; # target sum
var x {i in 1..n,j in 1..n,k in 1..N} binary; # binary variables
subject to NB1{j in 1..n}: sum{i in 1..n,k in 1..N} k*x[i,j,k] = C;
subject to NB2{i in 1..n}: sum{j in 1..n,k in 1..N} k*x[i,j,k] = C;
subject to NB3{i in 1..n,j in 1..n}: sum{k in 1..N} x[i,j,k] = 1;
subject to NB4{k in 1..N}: sum{i in 1..n,j in 1..n} x[i,j,k] = 1;
solve; # solution
var y {i in 1..n,j in 1..n}=sum{k in 1..N} k*x[i,j,k]; # transformed
display y; # output