# 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?

• Please use MathJax to format math on this site. To begin with, enclose all math expressions (including numbers) in $ signs. For example, $x_1^2$ will give you$x_1^2$. You'll get a much better response if your posts are easy to read. Commented Apr 8, 2021 at 22:38 • I've Mathjax'ed it, but OP could you please check that I didn't change the problem in doing so? Commented Apr 8, 2021 at 22:47 • Thank you @ConMan - I appreciate that and it's correct. Commented Apr 9, 2021 at 2:01 ## 3 Answers 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)  • Thank you, do you know of a Python library that can do the equivalent and ideally support constraints of the for ai*bi/(xi+ai) == aj*bj/(xj+aj) Commented Apr 12, 2021 at 18:21 • I would assume cvxpy should be able to do something equivalent, otherwise you just model the inverse manually using standard SOCP components in cvxpy. aibi/(xi+ai) == ajbj/(xj+aj) is much worse as it is a nonlinear nonconvex inequality. With aibi == ajbj/(xj+aj)(xi+ai) you at least land at a nonconvex quadratic equality. Still theoretically intractable but more support available if you go for global solver. Commented Apr 12, 2021 at 18:27 • I'm kind of familiar with LP and QP, but totally new to SOCP. I saw your code references yalmip.github.io/command/cpower but can't find much more information on what that does and in general what it means to model the inverse problem. Do you know a good reference/tutorial that might be helpful? Commented Apr 13, 2021 at 23:48 • Second-order cone programming is a generalization of linear programming and is implemented in most commercial packages previously only doing linear and quadratic programming. It allows for much more expressive power (such as modelling powers when convex) at almost the same computational complexity en.wikipedia.org/wiki/Second-order_cone_programming. Commented Apr 14, 2021 at 5:13 • The definition seems to imply that the objective function can only be linear as it's f^T.x in the wikipedia article. How can that be used to accommodate an objective function of the form (ab) / (a+x) Commented Apr 14, 2021 at 18:18 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. • How to impose the restrictions$x_i \ge 0$? Commented Apr 8, 2021 at 23:24 • Hmm, good point, I originally thought it would naturally occur by taking the positive square root when solving for the$x_i\$ but I can see that it's not guaranteed. I suspect that you will get the right solution if you try the above, force any negative x_i to 0 and resolve for the remaining values, and you can probably prove that with a convexity argument, but I'm not sure if there is a more efficient method. Commented Apr 9, 2021 at 0:34
• Unfortunately your suspicion is not true, at least not in general. Commented Apr 9, 2021 at 1:30

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.

• Charnes-Cooper does not work because of the coupling constraint (it almost never works when the objective is a sum of fractions). You can solve it as a convex optimization problem though. Commented Apr 8, 2021 at 22:49
• @LinAlg Agreed. I modified my answer. Commented Apr 8, 2021 at 23:14