How can I solve this form of optimization problem? I am looking to solve an optimization problem of $n$ real variables $x_1, \dots, x_n$.
Maximize: $\sum_{i=1}^n \left(b_i - \frac{a_i b_i}{x_i + a_i}\right)$
Such that:
$x_1 + \dots + x_n = c$
and
$x_i \geq 0$
Where $a_i$, $b_i$ and $c$ are positive constants.
I initially thought that I could solve this with linear programming but that's not the case I don't think.
I also tried piece-wise linear approximation but that consistently under-estimates the optimal solution.
Can this be solved some other way like as quadratic programming? Or is there a way to transform this problem and solve it as linear programming?
 A: The problem is convex so basically any nonlinear solver should solve this without issues. If you want to have something closer to linear programming, you can use the fact that $x^{-1}$ is second-order cone representable and thus use a second-order cone programming solver.
Here tested in MATLAB Toolbox YALMIP (requires an SOCP solver such as Mosek, Gurobi, SeDuMi, ECOS etc for the reformulated problem)
b = rand(n,1);
a = rand(n,1);
c = rand(1);
x = sdpvar(n,1);
% Solve as general nonlinear program
obj = sum(a.*b./(x+a));
optimize([x>=0, sum(x)==c],obj)
% Model inverse using socp cone
obj = sum(a.*b.*cpower(x+a,-1));
optimize([x>=0, sum(x)==c],obj)

A: If the variables must be integer-valued, you can solve the problem via dynamic programming as follows.  Let
$$f(k,c) = \max_{\substack{(x_1,\dots,x_k)\in \mathbb{Z}_+^k:\\x_1+\dots+x_k = c}} \sum_{i=1}^k \left(b_i - \frac{a_i b_i}{x_i + a_i}\right)$$
By conditioning on $x_k$, we obtain DP recursion:
$$f(k,c) = \max_{x_k\in\{0,\dots,c\}}\left\{b_k - \frac{a_k b_k}{x_k + a_k} + f(k-1,c-x_k)\right\}$$
You want to compute $f(n,c)$.
Also, the objective has a constant term $\sum_{i=1}^n b_i$, so you might as well instead minimize $$\sum_{i=1}^n \frac{a_i b_i}{x_i + a_i}$$

For the continuous problem, if you had $0 \le x_i \le c$ instead of the equality constraint or if your objective were a single ratio of linear functions (rather than a sum of ratios), you could linearize via a Charnes-Cooper transformation.
A: Assuming you want real-valued answers (and not, say, integers), this kind of question is amenable to a Lagrange multiplier approach:
Let $f(\vec{x}) = \sum_{i=1}^n \left( b_i - \frac{a_i b_i}{a_i + x_i} \right)$ and $g(\vec{x}) = c - \sum_{i=1}^n x_i$. Then define
$$\begin{eqnarray}\mathcal{L}(\vec{x}; \lambda) & = & f(\vec{x}) + \lambda g(\vec{x}) \\
& = & \sum_{i=1}^n \left(b_i - \frac{a_i b_i}{a_i + x_i}\right) + \lambda \left (c - \sum_{i=1}^n x_i \right ) \\
& = & \lambda c + \sum_{i=1}^n \left( b_i - \lambda x_i - \frac{a_i b_i}{a_i + x_i} \right) \end{eqnarray}$$
We then solve for $\nabla \mathcal{L} = 0$:
$$\begin{eqnarray} \frac{\partial \mathcal{L}}{\partial x_i} & = & -\lambda + \frac{a_i b_i}{(a_i + x_i)^2} \\
& = & 0 \mbox{ when } x_i = \sqrt{\frac{\lambda}{a_i b_i}} - a_i \\
\sum_{i=1}^n x_i & = & c \\
\sum_{i=1}^n \left( \sqrt{\frac{\lambda}{a_i b_i}} - a_i \right) & = & c \\
\sqrt\lambda \sum_{i=1}^n \frac{1}{\sqrt{a_i b_i}} & = & c + \sum_{i=1}^n a_i \\
\lambda & = & \left( \frac{c + \sum a_i}{\sum (a_i b_i)^{-\frac{1}{2}}} \right)^2 \\
x_i & = & \frac{c + \sum a_k}{\sqrt{a_i b_i} \sum (a_k b_k)^{-\frac{1}{2}}} - a_i
\end{eqnarray}$$
and you can substitute that into $f$ to get an expression for its extremum.
