# How to linearize a quadratic constraint to add it then via a callback function

Operations Research Asked by Farouk Hammami on December 9, 2021

Suppose we have a positive continuous variables $$0 le x le UB$$ where $$UB$$ is a known upper bound.

How can we linearize the term $$x^2$$?

Detailled problem:

Suppose that via a callback we compute a factor namely $$A_i in ]0,1]$$. After computing this factor: We need to add the following lazy constraint (using add(modeler,…)):

$$x^2_i le A_i^2 sum_k sigma^2_k y_{ki}$$; ($$x_i ge 0$$, $$y_{ki} in {0,1}$$ are decision variables and $$sigma_k > 0$$ are known parameters).

Adding this lazy constraint results in infeasible status given it is quadratic.

You might want to take a look at two blog posts I wrote earlier this year:

If you approximate $$x^2$$ via tangents, all feasible points will satisfy your lazy constraint, but there will be some infeasible points that satisfy it. If you approximate $$x^2$$ via secants, all points that satisfy the lazy constraint will be feasible, but it will cut off some feasible points. In either case, the more granular the approximation (the more intervals in your piecewise linear function), the closer you get to what you want.

The second post includes some Java code (using CPLEX).

Answered by prubin on December 9, 2021

let me adapt the interpolate example from

to x*x:

float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

int nbSegments=5;

float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
float y2[i in 0..nbSegments]=x2[i]*x2[i];  // y=f(x)

float firstSlope=0;
float lastSlope=0;

tuple breakpoint // y=f(x)
{
key float x;
float y;
}

sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};

float slopesBeforeBreakpoint[b in breakpoints]=
(b.x==first(breakpoints).x)
?firstSlope
:(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);

pwlFunction f=piecewise(b in breakpoints)
{ slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);

assert forall(b in breakpoints) abs(f(b.x)-b.y)<=0.001;

float maxError=max (i in 0..sampleSize) abs(x[i]*x[i]-f(x[i]));
float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(x[i]*x[i]-f(x[i]));

execute
{

// turn an OPL array into a python list
function getPythonListOfArray(_array)
{

var quote=""";
var nextline="\n";

var res="[";
for(var i in _array)
{
var value=_array[i];

if (typeof(value)=="string") res+=quote;
res+=value;
if (typeof(value)=="string") res+=quote;
res+=",";
res+=nextline;
}
res+="]";
return res;
}

// Display a function with points with x and y arrays of x and y
function displayXY(x,y,pythonpath,pythonfile)
{
writeln("displayXY ",x," ",y," ",pythonpath," ",pythonfile);

var python=new IloOplOutputFile(pythonfile);
python.writeln("import matplotlib.pyplot as plt");
python.writeln("x = ",getPythonListOfArray(x))
python.writeln("y = ",getPythonListOfArray(y))
python.writeln("plt.plot(x, y)");
python.writeln("plt.xlabel('x - axis')");
python.writeln("plt.ylabel('y - axis')");
python.writeln("plt.title('xy graph')");
python.writeln("plt.show()");
python.close();
IloOplExec(pythonpath+" "+ pythonfile,true);
}

}

int nbSegments2=10000;

float x3[i in 0..nbSegments2]=(s)+(e-s)*i/nbSegments2;
float y3[i in 0..nbSegments2]=x3[i]*x3[i];  // y=f(x)
float y3pwl[i in 0..nbSegments2]=f(x3[i]);  // y=f(x)

string pythonpath="C:\Python36\python.exe";
string pythonfile="C:\temp\DisplayXY.py";
execute
{

// display x*x function
displayXY(x3,y3,pythonpath,pythonfile);
// display pwl approximation
displayXY(x3,y3pwl,pythonpath,pythonfile);
}


and you will see and later on you may use f as square function:

dvar float xx;
dvar float yy;
subject to
{
xx==2;
yy==f(xx);
}

execute
{
writeln("yy=",yy);
}


gives

yy=4


Answered by Alex Fleischer on December 9, 2021

## Related Questions

### Conditions required for strong duality to hold for SDPs

1  Asked on August 19, 2021

### Relationship between extreme points and optimal solutions of SDPs

1  Asked on August 19, 2021

### How to parallelize metaheuristics algorithms (Island Model)?

1  Asked on August 19, 2021 by antarctica

### Can we get closed form solution for such a problem?

1  Asked on August 19, 2021

### Can we use reinforcement learning and convex optimization to solve an optimization problem?

2  Asked on August 19, 2021 by qinqinxiaoguai

### Nonlinear integer (0/1) programming solver

6  Asked on March 1, 2021 by rajya

### How can I find the shortest path for all nodes in a graph from a source $s$?

1  Asked on March 1, 2021 by windbreeze

### Free solver for MINP problems

1  Asked on February 18, 2021 by dspinfinity

### Where I can study some job shop scheduling by course (video )?

0  Asked on February 18, 2021 by yue-chao

### Linear objective function with power term in constraint

1  Asked on February 15, 2021 by user152503

### Formulating these logical constraint in an ILP

1  Asked on January 18, 2021

### Modeling the multiplication of two binary decision variables in undirected graph in python

0  Asked on January 18, 2021 by amedeo

### Flexible Job Shop with Preemption

0  Asked on January 15, 2021 by robert-hildebrand

### How to handle an equality constraint in metaheuristic algorithms (like GA, PSO)?

3  Asked on January 11, 2021 by stevgates

### What is the difference between min- cut formulation and (bi) partitioning formulation?

1  Asked on January 8, 2021 by fathese

### Logical constraint in ILP

1  Asked on December 22, 2020 by che

### Quasi-convex function must be “partially monotonic”?

1  Asked on December 13, 2020 by high-gpa

### Constraint programming resources

3  Asked on November 28, 2020 by joffrey-l

### Pyomo variable creation dilemma

1  Asked on October 31, 2020 by ethan-deakins

### Convexity of the variance of a mixture distribution

1  Asked on September 25, 2020 by independentvariable