Your implementations of constraints (1), (6) and (7) do not match the algebraic formulation. (You sum over different variables than what the formulation indicates.) Does fixing that solve the problem? Addendum: After modification, the new solution uses only variables $x_{i,j,k}$ where $i=j$ (and consequently $c_{i,j}=0.$ Without knowing the full details of the model, I suspect that you want either to constrain $x_{i,i,k}=0$ for all $i$ or change the index of summation in (1) (and possibly elsewhere) to skip the case $i=j.$