Convex.jl affine constraint not respected in solution

round rounds to the nearest integer, so you’re just rounding rho away to the zero matrix. With SCS (I do not have a Mosek license), I get

julia> evaluate(rho)
4×4 Matrix{ComplexF64}:
 0.25+0.0im   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im  0.25+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im  0.25+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  0.25+0.0im

julia> tr(evaluate(rho))
0.999999966223076 + 0.0im