I have the following problem, abstracting away a few details from a real-world application, that I want to solve with linear programming (or any other mathematical optimization with constraints, really).
I have two sets of values $y$ and $\hat{r}$. $y$ are real-world measured values, while $\hat{r}$ are outputs from a simulation with mathematical model. $y$ are integers in range $[0, 1, \dots, 21]$, while $\hat{r}$ are real numbers in range $[0, 21]$ that are rounded to integers $\hat{y}$. So far I used regular rounding, e.g. $1.33 \rightarrow 1$, $2.5 \rightarrow 3$, $21.0 \rightarrow 21$. This can also be written as $\text{round}(x) = \lfloor x+0.5 \rfloor$.
My $y$ values distribution is far from uniform (more like log-normal), which is reflected by the model, which often outputs values too large (for small $y$) or too small (for large $y$). My idea is to "optimize" the rounding, changing the range $[a_i, a_{i+1}]$ for which the number is rounded to $i$. I minimize the sum of absolute deviations, i.e. $\sum_{j=1}^n |y_j - \hat{y}_j|$ (for $n$ observations). For each $\hat{r}_j$, I have $\hat{y}_j = \min_{i} \left|\hat{r}_j - \frac{a_{i+1} - a_i}{2}\right|$, i.e. I find the range $[a_i, a_{i+1}]$ for which the center is the closest to $\hat{r}_j$.
Additionally, for thresholds $a_i$ I have constraints $a_0 = 0$, $a_{21} = 21$, and $a_i < a_{i+1}$ for all $i = 0, 1, \dots, 20$.
My formulation so far is:
$$ \text{minimize } \sum_{j=1}^n \left|y_j - \hat{y}_j\right| \\ \text{for variables } a_i, i = 0, 1, \dots, 21 \\ \text{where } \hat{y}_j = \arg \min_{i} \left|\hat{r}_j - \frac{a_{i+1} - a_i}{2}\right| \\ \text{subject to:} \\ a_0 = 0 \\ a_{21} = 21 \\ a_i < a_{i+1}, i = 0, 1, \dots, 20 $$
The problem lies in formulation of $\hat{y}_j$, which is a nonlinear function. I know that this could be optimized with any blackbox zero-order optimizer, but I am pretty sure this can be reformulated to use linear programming efficiently.
I also had an idea to directly optimize $\sum_{j=1}^n \left|\hat{r}_j - \frac{a_{i+1} - a_i}{2}\right|$, make it an integer and use MIP (Mixed Integer Programming), but I'm not sure if that would be an equivalent formulation.
Example data, for simpler case with $j=4$:
y,r
0,0.1
0,0.3
1,0.45
1,0.49
1,0.65
1,0.67
1,0.78
1,0.99
1,1.1
1,1.27
1,1.7
2,1.9
2,2.1
3,3.3
3,3.4
4,3.5
4,3.99
With regular rounding, i.e. $a_0 = 0$, $a_1 = 0.5$, $a_2 = 1.5$, $a_3 = 2.5$, $a_4 = 3.5$, $a_5 = 4$, I get [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4], with sum of absolute deviations 3.
However, predictions here are heavily skewed towards 1. For thresholds $a_0 = 0$, $a_1 = 0.4$, $a_2 = 1.75$, $a_3 = 2.5$, $a_4 = 3.5$, $a_5 = 4$, I get [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 4] and sum of absolute deviations 0.
You can solve the problem via linear programming by reformulating as a shortest path problem in a layered network defined as follows. Let $R$ be the set of distinct $r$ values. Let $I=\{1,\dots,4\}$. The main nodes are $R \times I$, and the main arcs are from $(r,i)$ to $(r',i+1)$, where $r<r'$. The idea is that arc $(r,i)\to(r',i+1)$ is traversed if all observations with $r \le r_j < r'$ have $\hat{y}_j=i$. The cost for traversing the arc is the sum of absolute deviations $$\sum_{j: r \le r_j < r'} \left|y_j - i\right|.$$ Add a source node $(0,0)$, with directed arcs to all $(r,1)$, and a sink node $(4,5)$, with directed arcs from all $(r,4)$. Now find a shortest path from source to sink. For each arc $(r,i)\to(r',i+1)$ that appears in the shortest path, take $a_i=r$.