One approach you could do is to create a series of projected quadratic equality constraints:
model = Model()
@variable(model, A[1:2, 1:2])
@variable(model, B[1:2, 1:2])
@variable(model, C[1:4, 1:4])
@variable(model, D[1:4, 1:4])
@constraint(model, D .== kron(A, B))
@variable(model, E[1:16, 1:16])
@constraint(model, E .== kron(D, C)) # kron(A, B, C)
@objective(model, Min, tr(E))
You could also write the kron expression as an explicit outer product with sums, similar to what you’re doing.
I’ll have a look at the WIP branch.
To clarify, you shouldn’t use this. It’s still under development.