Proof of NP-completeness via extra information
No.  Your argument is not valid; it doesn't prove NP-completeness.  You have shown that an ILP solver can be used to solve your problem.  But an ILP solver can be used to solve easy problems, too, so it doesn't rule out the possibility that your problem might be easy (solvable in polynomial time).  The reduction needs to go the other way.  I suggest studying standard material on NP-completeness and reductions, e.g., What is the definition of P, NP, NP-complete and NP-hard? and What are common techniques for reducing problems to each other? and How do I construct reductions between problems to prove a problem is NP-complete?.
As it happens, I believe your problem can be solved in polynomial time; in other words, I believe it is not hard.  I believe the following algorithm solves your problem:

Enumerate over all possibilities $m_L,m_R$ for $\text{median}(K(L))$ and $\text{median}(K(R))$.  For each such possibility:

Find the partition $L,R$ that minimizes $E(L)+E(R)$, subject to the requirement that $\text{median}(L)=\mu_L$ and $\text{median}(R)=\mu_R$.


Output the best solution found at any of these stages.


I show below how to do the innermost step in polynomial time.  This gives a polynomial-time algorithm for your problem, as there are only $O(n^2)$ possibilities for $m_L,m_R$ (i.e., each $m_L,m_R$ must be an element of $\cup_i X_i$).
We'll use dynamic programming for this.  Let $\text{imb}(T,m)$ be defined by
$$\text{imb}(T,m) = (\sum_{t \in T} \text{sgn}(t-m), \sum_{t \in T} (t=m))$$
where the sum is &quot;with repetition&quot;, and where
$$\text{sgn}(x) = \begin{cases}
1 &amp;\text{if }x&gt;0\\
0 &amp;\text{if }x=0\\
-1 &amp;\text{if }x&lt;0
\end{cases}$$
$$(x=y) = \begin{cases}
1 &amp;\text{if }x=y\\
0 &amp;\text{otherwise}
\end{cases}$$
Define the imbalance of a partition $L,R$ with respect to $m_L,m_R$ to be the pair $(\vec{b_L},\vec{b_R})$ where $\vec{b_L}=\text{imb}(L,m_L)$ and $\vec{b_R}=\text{imb}(R,m_R)$.  Then we can reformulate the problem as follows:
Given $m_L,m_R$, find $L,R$ that minimizes $E(L)+E(R)$ subject to the requirement that the imbalance of $L,R$ with respect to $m_L,m_R$ is $(0,0)$ (or $+\infty$ if no such partition exists).
Now we can solve that using dynamic programming.  Define $A[k,\vec{b_L},\vec{b_R}]$ to be the smallest possible value of $E(L)+E(R)$ subject to the requirements that $L,R$ are a partition of $X_1,\dots,X_k$ with imbalance $(\vec{b_L},\vec{b_R})$ with respect to $m_L,m_R$ (or $+\infty$ if no such partition exists).  Then $A$ satisfies the recurrence relation
$$A[k,\vec{b_L},\vec{b_R}] = \min(A[k-1,\vec{b_L}-\text{imb}(X_k,m_L),\vec{b_R}] + \sum_{i \in X_k} |i - m_L|,A[k-1,\vec{b_L},\vec{b_R}-\text{imb}(X_k,m_R)] + \sum_{i \in X_k} |i - m_R|).$$
Using this recurrence, you can fill in the array $A$ in $O(Kn^4)$ time, where $n=|\cup_i X_i|$ is the number of elements (counting repetitions): you just initialize every entry to $+\infty$ except that $A[0,0,0]=0$, then fill it in, in order of increasing $k$.  Then, scan  all $A[K,(a,b),(c,d)]$ such that (i) $|a|&lt;b$ and $b&gt;0$, or $a=0$, and (ii) $|c|&lt;d$ and $d&gt;0$, or $c=0$; the smallest such value is the answer to the original problem.