I want to find $\delta_j$ in the following optimization problem. My variables are $\gamma_i$ and $\delta_j$ (all other symbols are known parameters). Assume $i\in\{1,\ldots,9\}$ and $j=\{1,\ldots,3\}$. Additionally, let $ASN = [asn_i] = [1, 1, 1, 2, 2, 2, 3, 3, 3]$. I think there must be an analytical solution for $\delta_j$ using lagrangian multiplier method, however I cannot find it, since I cannot get rid of the absolute functions within the summation operator.
Some facts I know: $C_j$ is the arithmetic mean of assigned $x_i$, i.e., $C_j = \frac{\sum\limits_{asn_i=j}^{}{x_i}}{3}$. Moreover, I know $\sum\limits_{i=1}^9{x_i} = \sum\limits_{j=1}^3 C_j = 0$.
$f(\delta,\gamma)=minimize \sum\limits_{i=1}^{9}{\gamma_i^2}+3\sum\limits_{j=1}^3{\delta_j^2}$
s.t.
$\sum\limits_{j=1}^{3}{\delta_j} = 0$
$\sum\limits_{j=1}^3{(C_j+\delta_j)}^2 = 3$
$|x_i-C_{asn_i}-\delta_{asn_i}|=1-\gamma_i\ \ \ \forall i\in\{1,\ldots,9\}$
EDIT:
For my application, I also can change the last set of constraints to the following set:
$(x_i-C_{asn_i}-\delta_{asn_i})=1-\gamma_i\ \ \ \forall i\in\{1,\ldots,9\}$
Well, this is not a finished solution. Your function is problematic, due to the absolute value, and I'm not sure whether there is something to be done regarding it. If the function is not differentiable, you probably won't be able to make it so. However, I do believe that the problem can be reduced from the 12 variables you have now, to a single one. From there, you could at least see if there are manageable domains where the expression inside the absolute value preserves it's sign, and differentiate in each region... Or something else. In any case, 1 variable is preferable to 12.
As a warning, I may make an error in precise numbers, but the geometrical idea behind what follows is still valid.
Notation-wise, whenever I omit an index, I mean the corresponding vector. So $\delta = (\delta_1,\delta_2,\delta_3)$.
First of all, you last set of constraints defines $\gamma_i = 1 - |x_i - C_{{asn_i}} - \delta_{asn_i}|$, so there go 9 variables already.
We are left with $\delta_{1,2,3}$ subject to $\sum \delta_i = 0$ and $\sum (\delta_i + C_i)^2 = 3$. Fortunately, the first is an equation of a plane, and the second - of a sphere with center $-C$. Moreover, $\sum C_i = 0$, so the center of the sphere lies on the plane, which means that the intersection is not empty, and in fact is a circle of radius $\sqrt 3$. This is great news, as a circle can be parametrized with a single variable. If we could parametrize it, we would be left with a much simpler expression. However, the way they are presented now, finding the parametrization would be hard.
Fortunately, since both the origin point and the center of the sphere lie on our plane, we can shift the sphere to the origin, without perturbing the plane itself. Thus, define $$\delta = -C + v$$ Now, $v_i$ are located on a $\sqrt 3$ sphere around the origin, and the set we minimize over becomes the intersection of that sphere with the plane $\sum v_i = 0$. Still, parametrization is hard. But! The sphere is around the origin, so rotations of the whole space will preserve the sphere. We can rotate the space such that the plane we are interested in becomes the $X-Y$ plane. Define: $$z = Av$$ Where $A$ is an orthogonal matrix that does precisely the said rotation. If I am not mistaken, $A$ is the matrix that rotates the polar angle by $\theta = \arcsin(\frac1{\sqrt 3}) \approx 35^o$ counterclockwise, and preserves the azimuth angle.
In terms of $z$ the resulting intersection is just a $\sqrt 3$ circle in the $X-Y$ plane, and therefore has the parametrization of $$ z(t) = (\sqrt 3 \cos(t), \sqrt 3 \sin(t), 0) $$
All the transformations we did along the way are nice and invertible, so we can write, finally the parametrization for $\delta$: $$\delta(t) = -C + A^{-1} \begin{pmatrix} \sqrt 3 cos(t)\\ \sqrt 3 sin(t)\\ 0\\ \end{pmatrix} \text{, for t $\in [0,2\pi) $}$$
This, together with expressing $\gamma$ via $\delta$ should reduce you to a problem of optimizing a single variable function over an interval... But the function still remains not very nice. I'm not sure what can be done with it further, other than exploring the sign of the expression inside the absolute value, which should be possible, but tedious. Or, at the very least you could plot the function now, to get some intuition as to how it behaves.
Following the discussion in the comments to this answer, and swapping the absolute value for square in the constraints, here is a sketch of the continuation to the above solution. I manage to get to a trigonometric equation, but do not seem to be able to solve it.
From above, we can write $$\delta_i = a_i \cos(t) + b_i \sin(t) + c_i$$ $$\gamma_i = 1 - (x_i - C_{asn_i} - \delta_{asn_i})^2$$
$$f(t) = \sum_1^9 1 - (x_i - C_{asn_i} - \delta_{asn_i})^2 + 3 \sum_1^3 \delta_i^2 = $$ $$- \sum_1^3 (p_i - \delta_1)^2 - \sum_1^3 (q_i - \delta_2)^2 - \sum_1^3 (s_i - \delta_3)^2 + 3 \sum_1^3 \delta_i^3 + 9$$ (with the constants inside the brackets defined accordingly) $$\dot f(x) = -\sum_1^3 (2 \dot \delta_1 \delta_1 - 2p_i \dot \delta_1) -\sum_1^3 (2 \dot \delta_2 \delta_2 - 2q_i \dot \delta_2) -\sum_1^3 (2 \dot \delta_3 \delta_3 - 2s_i \dot \delta_3) + \sum_1^3 2 \dot \delta_i \delta_i$$
For any $i$:
$$\dot \delta_i \delta_i = (-a_i \sin t + b_i \cos t)(a_i \cos t + b_i\sin t + c_i) =$$ $$a_ib_i(\cos^2 t - \sin^2 t) + (b_i^2 - a_i^2)\sin t\cos t - a_ic_i \sin t + a_ic_i \cos t =$$ $$a_ib_i\cos 2t + \frac{(b_i^2 - a_i^2)}2 \sin 2t - a_ic_i \sin t + a_ic_i \cos t$$
Now, you substitute that and $\dot \delta_i = -a_i \sin t + b_i \cos t$ into the expression, and from linearity all the coefficients before the sines and cosines will bundle up, which will leave you with an expression of the form:
$$A \cos 2t + B \sin 2t + C \cos t + D \sin t = 0$$
And I just do not know what to do with that... It is possible to find the precise $A,B,C,D$ in terms of $x_i, C_i$, but barring a case where one of them vanishes, or $C = D$ of some such thing, I don't see how it would help.