We have two decision variables \begin{align} & 0<x\le X,\\ & 0<y\le Y, \end{align} where both $X$ and $Y$ are two sensible upper bounds on our decision variables.

We also have a constraint


We discretize the interval $(0,X]$ and denote each piece by $r_i$ where $i=1,\ldots,n$ and $n$ is a finite number.

By defining a binary decision variable say $z_i\in\{0,1\}$ where $z_i=1$ if and only if $x\in(r_{i},r_{i+1}]$ and $z_i=0$, otherwise, we linearize this constraint as


Now the question is how to determine $n$?

I am also aware that there are other techniques of linearization we would like to focus on this method for the moment.

I'm also sure this method has been used before but I could not find any reference on it, mainly because I don't know if it has a name or where to even look!

Could someone help please?

  • $\begingroup$ Why not to set $n$ equal to $X$? $\endgroup$
    – Mostafa
    Jan 17, 2022 at 16:31
  • $\begingroup$ $n$ is a positive integer number and $X$ is a continuous number. $\endgroup$
    – user9659
    Jan 17, 2022 at 16:37
  • $\begingroup$ This does not address your approach, but in this particular case, since $x^2/(1-x)$ is strictly increasing when $x,y>0$, the constraints can be summarised by $0<x<\min\left\{X,\frac{\sqrt{Y^2+4Y}-Y}2\right\}$. $\endgroup$
    – TheSimpliFire
    Jan 17, 2022 at 16:51
  • $\begingroup$ As @prubin points out, it must be that $1-x > 0$. So f $y$ does not appear anywhere in your problem other than the constraint $0 < y \le Y$, then $y$ need never be formed, and that constraint can be replaced by the two constraints, $x \le (1-x)Y,-x \le (1-x)Y$ $\endgroup$ Jan 17, 2022 at 20:13

2 Answers 2


I assume $X\le1$, since the formulation is infeasible for $x > 1$.

What you are doing is a form of piecewise linear approximation of the function $f(x) = x^2/(1-x)$. Choosing $n$ is a matter of trial and error. Obviously, the larger $n$ is the more accurate the approximation is. Also obvious is that larger $n$ means longer solution time (and more memory usage).

There are a couple of things worth noting. One is that you might (or might not) get a more accurate approximation using $(r_i + r_{i+1})/2$ rather than $r_i$ as the point in interval $i$ where you evaluate $f(x)$. Another is that it is often a good idea to concentrate evaluation points where the function has greatest curvature. So you would shrink the spacing between points as $x$ progressed from 0 to 1.

Addendum: If solution times are fairly brief, one other thing that might be worth trying is to discretize, solve, add break points near the optimal solution (to get a better approximation there) and maybe remove a few breakpoints away from the optimum (to conserve model size), solve again, and repeat until satisfied (or bored).

  • $\begingroup$ Thank you. What does this sentence mean "... concentrate evaluation points where the function has greatest curvature" ? $\endgroup$
    – user9659
    Jan 17, 2022 at 16:41
  • 1
    $\begingroup$ Rather than using an evenly spaced grid, put the grid points ($r_i$) closer together where $f(x)$ is changing rapidly (toward $x=1$) and put the grid points points further apart where $f(x)$ is changing slowly (toward $x=0$). $\endgroup$
    – prubin
    Jan 17, 2022 at 17:34
  • $\begingroup$ Thank you for the clear explanation. So let's say $\ell_i$ is the length of the $i$-th interval. Should we set a coefficient $\alpha_i$ and define $\ell_i = \alpha_i \frac{X}{n}$? I'm confused on how to reformulate the constraint above using this idea. $\endgroup$
    – user9659
    Jan 17, 2022 at 17:46
  • 1
    $\begingroup$ There is no single correct way to do this. If I had to prescribe an approach, it might be to choose interval lengths so that the increase in $f(x)$ from left to right endpoint is approximately equal on each interval. $\endgroup$
    – prubin
    Jan 17, 2022 at 18:59

Have you tried solving the problem directly? This is simply a quadratic non-convex constraint:

$$ y = \frac{x^2}{1-x}\implies y - yx = x^2 $$

So e.g. using Gurobi and Python, you can write:

import gurobipy as gp

m = gp.Model()
x = m.addVar(0, x_ub)
y = m.addVar(0, y_ub)
m.addConstr(y - y*x == x*x)
  • $\begingroup$ Thank you. This is helpful when I code the problem in python but $y$ appears in other constraints as well as in the objective function. $\endgroup$
    – user9659
    Jan 18, 2022 at 10:28
  • $\begingroup$ How is the fact that $y$ enters the objective or other constraints in any way change the fact that you can model the relation as a quadratic constraint. If the objective and all other constraints are at most quadratic in $x$, $y$ and any other variables, you still have a nonconvex quadratic model. $\endgroup$ Jan 18, 2022 at 19:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.