As @ManfredWeis noted, the objective is equivalent to minimizing $\sum_{i=1}^K \frac{x_i^2}{y_i}$. You can reformulate this as a (convex) mixed integer second-order cone programming (MISOCP) problem by introducing a new variable $z_i$ and minimizing $2\sum_i z_i$ subject to
\begin{align}
\sum_i x_i &= X \\
\sum_i y_i &= Y \\
x_i &\le y_i &&\text{for all $i$} \\
2 z_i y_i &\ge x_i^2 &&\text{for all $i$}
\end{align}
The new constraint is a rotated second-order cone.
Note that, even in the original problem, relaxing integrality yields unique optimal solution $(x_i^*,y_i^*)=(X/K,Y/K)$ for all $i$, with objective value $0$. Because the relaxation is convex, every integer optimal solution satisfies $x_i \in \{\lfloor x_i^* \rfloor, \lceil x_i^* \rceil\}$ and $y_i \in \{\lfloor y_i^* \rfloor, \lceil y_i^* \rceil\}$ for all $i$.
An alternative reformulation introduces binary decision variables $u_{ij}$ for $i\in\{1,\dots,K\}$ and $j\in \{0,\dots, X\}$ and $v_{ik}$ for $i\in\{1,\dots,K\}$ and $k\in \{1,\dots, Y\}$. The problem is to minimize the quadratic function
$$\sum_{i=1}^K \sum_{j=0}^X \sum_{k=1}^Y \frac{j^2}{k} u_{ij} v_{ik} \tag1\label1$$
subject to linear constraints
\begin{align}
\sum_i x_i &= X \\
\sum_i y_i &= Y \\
x_i &\le y_i &&\text{for all $i$} \\
\sum_j u_{ij} &= 1 &&\text{for all $i$} \\
\sum_j j u_{ij} &= x_i &&\text{for all $i$} \\
\sum_k v_{ik} &= 1 &&\text{for all $i$} \\
\sum_k k v_{ik} &= y_i &&\text{for all $i$}
\end{align}
You can call an MIQP (or BQP) solver.
Alternatively, you can linearize \eqref{1} by introducing nonnegative variables $w_{ijk}$ to represent the product $u_{ij} v_{ik}$.
You can then minimize $$\sum_{i=1}^K \sum_{j=0}^X \sum_{k=1}^Y \frac{j^2}{k} w_{ijk} \tag{1p}\label{1'}$$
The usual linearization imposes linear constraints
\begin{align}
w_{ijk} &\ge u_{ij} + v_{ik} - 1 \tag2\label2 \\
w_{ijk} &\le u_{ij} \tag3\label3 \\
w_{ijk} &\le v_{ik} \tag4\label4
\end{align}
But you can omit \eqref{3} and \eqref{4} because the objective will drive these constraints to be satisfied naturally.
A more compact linearization instead replaces \eqref{2} with
\begin{align}
\sum_j w_{ijk} &= v_{ik} &&\text{for all $i$ and $k$} \\
\sum_k w_{ijk} &= u_{ij} &&\text{for all $i$ and $j$}
\end{align}
It is also worth noting that without the linking constraints $\sum_i x_i = X$ and $\sum_i y_i = Y$ the problem decomposes into independent problems, one for each $i$. So Dantzig-Wolfe decomposition or Lagrangian relaxation might perform well.
