7
$\begingroup$

So I have an optimization problem of the following form:

\begin{aligned} \max_{x,y} \quad & \sum_i x_i \\ \text{s.t.} \quad & \sum_i x_iy_i \leq a \\ \quad & x_{\min} \leq x \leq x_{\max} \\ \end{aligned}

where $x$ is a vector of real (positive) variables, $y$ is a vector of binary variables ($y_i \in \{0, 1\}$), and $a$ is a positive constant. The objective function is linear, but the constraint is bilinear.

One way to solve this problem is the so-called "mountain-climbing" method where I first fix $y$ and solve for $x$, and then use the solution $x^*$ to to solve for $y$, and continue the process until some convergence condition is achieved on the objective.

My first question is, is there an easier way to solve such a problem? (i.e., convexifying the constraint somehow).

My second question is, if I'm using this "mountain-climbing" method, when I obtain a solution $x^*$ and then solve the next problem for $y$, do I just set $x^*$ for the constraint and keep the objective as optimization function, or do I simply solve for $y$ using the constraint with no objective?

Your help is much appreciated!

Improve this question
$\endgroup$
2
  • $\begingroup$ Indeed, this is not the real problem. There are some hard constraints for the value of $x$. Ill modify the problem so that it's clear. Thanks! $\endgroup$
    – Johnny
    Mar 31, 2021 at 12:32
  • $\begingroup$ Hiya! It looks like you've created two accounts by accident - if you'd like to have all of your content on a single account, please consider visiting the contact page where you can merge your accounts. :) $\endgroup$
    – Catija ♦
    Apr 8, 2021 at 20:37

2 Answers 2

Reset to default
10
$\begingroup$

You can linearize the constraint by introducing a nonnegative variable $z_i$, replacing $x_i y_i$ with $z_i$, and imposing linear big-M constraints \begin{align}x_\min y_i &\le z_i \le x_\max y_i\tag1\\ (0-x_\max)(1-y_i) &\le z_i - x_i \le (0-x_\min)(1-y_i)\tag2 \end{align} Constraint $(1)$ enforces $y_i=0\implies z_i=0$. Constraint $(2)$ enforces $y_i=1\implies z_i=x_i$.

Because @Nikos mentioned automated reformulation, I should mention that the LINEARIZE option in SAS will do this for you:

con C1: sum {i in ISET} x[i] * y[i] <= a;
solve linearize;

Alternatively, you can use indicator constraints:

con C1: sum {i in ISET} z[i] <= a;
con C2 {i in ISET}: y[i] = 0 implies z[i] = 0;
con C3 {i in ISET}: y[i] = 1 implies z[i] = x[i];
solve;

If you have $$x_\min q_i \le x_i \le x_\max q_i,$$ the new constraints are instead \begin{align}0 &\le z_i \le x_\max y_i\tag{$1'$}\\ (0-x_\max)(1-y_i) &\le z_i - x_i \le 0\tag{$2'$} \end{align}

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks for this response! Just a quick question...If my second constraint ($x_{\min} \leq x \leq x_{\max}$) changes to $x_{\min}q \leq x \leq x_{\max}q$ (where $q$ is yet another vector of binary variables with $q_i \in \{0, 1 \}$), how would the linearization work? $\endgroup$
    – Johnny
    Apr 1, 2021 at 9:13
  • $\begingroup$ Also, is the second constraint (2) correct? If you replace $z_i = x_i y_i$ in this constraint, you would get $ - (x_{\min} - x_{\max}) \leq x_i \leq - (x_{\max} - x_{\min})$, which doesn't make sense. Is my logic correct? $\endgroup$
    – Johnny
    Apr 1, 2021 at 9:27
  • $\begingroup$ I corrected the big-M values in (2). The lower end was too aggressive. $\endgroup$
    – RobPratt
    Apr 1, 2021 at 13:34
  • $\begingroup$ Interesting...so the $q_i$ variables don't show up in the new constraints? Why is that? $\endgroup$
    – Johnny
    Apr 1, 2021 at 14:26
  • 1
    $\begingroup$ No, (1') and (2') enforce $z_i = x_i y_i$. You still need to explicitly impose $x_\min q_i \le x_i \le x_\max q_i $. $\endgroup$
    – RobPratt
    Apr 1, 2021 at 18:34
7
$\begingroup$

The best way to solve this in general is to reformulate this to an MILP. Although the reformulation itself is easy, it's also incredibly easy to make a mistake/forget to take something into account.

Since this is a toy problem, I will provide a solution that will work for any problem that has the structure $xy, y\in\{0,1\}$, including your actual problem.

The reformulation

$xy$ can be linearized by replacing $xy$ with an auxiliary variable and adding the following constraints:

\begin{aligned} x_ly-w\le0\\ w-x_uy\le0\\ w-x+x_l(1-y)\le0\\ -w+x-x_u(1-y)\le0\\ xy_l\le w\le xy_u \end{aligned}

where $x_l,x_u$ are parameters specifying the lower and upper bounds of $x$.

Reformulating automatically

The easiest way (and the only one that is free) to do this automatically for any problem where this structure is present, is to use our own Octeract Reformulator, which is free to use. It will perform the linearisation for every bilinear term even if your problem only partially has this structure and the rest is e.g. logarithms or cosines or other complicated functions.

You just need to install Octeract Engine (which is also free), grab the reformulation file from our GitHub repo, and copy-paste the following in a Python file:

from octeract import *
from LinearizeBilinearContBinary import LinearizeBilinearContBinaryMod

m = Model()
m.import_model_file('model.nl')

m.apply_mod(LinearizeBilinearContBinaryMod)

print(m)

m.write_problem_to_NL_file('/my/file/name.nl')
#m.write_problem_to_LP_file('/my/file/name.lp')
#m.write_problem_to_MPS_file('/my/file/name.mps')
#m.write_problem_to_GAMS_file('/my/file/name.gms')

This will produce a new file in your desired format which will contain the reformulated problem. You can then solve this with any solver/modelling language you want.

For example, the following Octeract Shell script creates a dummy model to try out the reformulation:

from octeract import *

# ==== Create a toy model ====
m = Model()

for i in range(1,4):
    m.add_variable("y"+str(i),0,1,BIN)

m.minimize("x1")

m.set("x1*y1+x2*y2+x3*y3").to(4)

for i in range(1,4):
    m.set_variable_bounds('x'+str(i),-1000,1000)
print(m)
m.write_problem_to_NL_file('model.nl')
# ============================

The print command will give us:

Structure : MBQCQP
Convexity : nonconvex
------------------------------------------
var y1 binary  >= 0, <= 1;
var y2 binary  >= 0, <= 1;
var y3 binary  >= 0, <= 1;
var x1 >= -1000, <= 1000;
var x3 >= -1000, <= 1000;
var x2 >= -1000, <= 1000;

minimize obj : +1*x1+0;

subject to

constraint : (x3*y3)+(x1*y1)+(x2*y2)+0.0 = 4;
------------------------------------------

and if we then run the very first script on the model.nl file we just created, we get:

Structure : MILP
Convexity : convex
------------------------------------------
var x1 >= -1000, <= 1000;
var x2 >= -1000, <= 1000;
var x3 >= -1000, <= 1000;
var x4 binary  >= 0, <= 1;
var x5 binary  >= 0, <= 1;
var x6 binary  >= 0, <= 1;
var w_xy_0 >= -1000, <= 1000;
var w_xy_1 >= -1000, <= 1000;
var w_xy_2 >= -1000, <= 1000;

minimize obj : +1*x1+0;

subject to

con1 : 1.0*w_xy_2+1.0*w_xy_1+1.0*w_xy_0+0.0 = 4;
constraint_0 : -1.0*w_xy_0+-1000.0*x6+0.0 <= 0;
constraint_1 : 1.0*w_xy_0+-1000.0*x6+0.0 <= 0;
constraint_2 : -1.0*x3+1000.0*x6+1.0*w_xy_0+0.0 <= 0;
constraint_3 : -1.0*w_xy_0+1000.0*x6+1.0*x3+0.0 <= 0;
constraint_4 : -1.0*w_xy_1+-1000.0*x4+0.0 <= 0;
constraint_5 : 1.0*w_xy_1+-1000.0*x4+0.0 <= 0;
constraint_6 : -1.0*x2+1000.0*x4+1.0*w_xy_1+0.0 <= 0;
constraint_7 : -1.0*w_xy_1+1.0*x2+1000.0*x4+0.0 <= 0;
constraint_8 : -1.0*w_xy_2+-1000.0*x5+0.0 <= 0;
constraint_9 : 1.0*w_xy_2+-1000.0*x5+0.0 <= 0;
constraint_10 : -1.0*x1+1.0*w_xy_2+1000.0*x5+0.0 <= 0;
constraint_11 : -1.0*w_xy_2+1000.0*x5+1.0*x1+0.0 <= 0;
------------------------------------------

What the reformulation file does

The file uses the Reformulator Python API to encapsulate the mathematical logic, so that the same logic can be applied to any problem. I have added comments in the mod to describe the steps:

from octeract import *

# Linearize bilinear term x*y where y is binary
# =============================================
# x_l*y-w<=0
# w-x_u*y<=0
# w-x+x_l*(1-y)<=0
# -w+x-x_u*(1-y)<=0
# xy_l<=w<=xy_u
# =============================================

# Define symbolic trigger
my_trigger = Match('V(x)*V(y)') 

# Filter binaries and account for permutations - we bind binary to b and non-binary to nb
xbin_filter = (IsBinary('x') & ~ IsBinary('y') & Bind('x', 'b') & Bind('y', 'nb'))
ybin_filter = (IsBinary('y') & ~ IsBinary('x') & Bind('y', 'b') & Bind('x', 'nb'))

# Combine the filters
my_filter = (xbin_filter | ybin_filter)

# Specify how to change the model
# Add auxiliary variable with the same bounds as x*y
add_auxiliary_var = AddVariableSpan('w_xy','b*nb')
# Add parameters for the bounds of the continuous variable
add_parameters = AddParameter('nb_LB','nb','lb') + AddParameter('nb_UB','nb','ub')
substitute_term = SubWith('w_xy')
add_constraint0 = AddConstraint('nb_LB*b-w_xy <= 0')
add_constraint1 = AddConstraint('w_xy-nb_UB*b <= 0')
add_constraint2 = AddConstraint('w_xy-nb-nb_LB*(1-b) <= 0')
add_constraint3 = AddConstraint('-w_xy+nb-nb_UB*(1-b) <= 0')

add_constraints = add_constraint0 + add_constraint1 + add_constraint2 + add_constraint3

# Add all modifications 
my_actions = add_auxiliary_var + add_parameters + substitute_term + add_constraints

# Create mod
LinearizeBilinearContBinaryMod = (my_trigger & my_filter).then(my_actions)

The best part about this approach is that you only need to figure out the reformulation logic once. You can then keep using the same mods on your problems or mixing and matching rules until the end of time.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the detailed answer! I'll check out this function and see how it works. $\endgroup$
    – Johnny
    Apr 1, 2021 at 18:41

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.

Not the answer you're looking for? Browse other questions tagged or ask your own question.