Well, when you start computing the envelope etc, you are essentially reinventing a global solver. That is what is performed in the first iteration. After that, the feasible space is partitioned, and envelopes (or approximations of) over the new regions are computed. In every region an LP relaxation is solvers, and a local nonlinear solve, in order to get lower and upper bounds. The information obtained from the first iteration is very bad, as the lower bound on the objective is 0.
Your new objective would be modeled with
magnitudes = s1.^2+s2.^2 ;
Constraints = [sqrtm(magnitudes) <= t, diff(xn)>=0,s1 == sum(cos(omega*xn*r),1),s2 == sum(sin(omega*xn*r),1)];
BTW, in case you wonder, I am avoiding the use of the max operator (and thus do it manually through t), since it is would fail a convexity propagation (max of sqrt of trigonometric is not structurally convex). Same thing with norm. sqrtm does not perform convexity propagation