Linearization binary variables / Block-Like Matrix / Indicator Constraint or L2-Norm / Additional SOS1
AnsweredDear all,
I am currently encountering long runtimes in optimization problem would be grateful for advice on how to accelerate the solution.
The following is an MRE with less constraints & objective functions than the original MOO but can reproduce the performance issues.
The MRE describes the decomposition of a 2D object, such as a line in multiple segments. The 2D object is discretized into multiple element \(n\). The elements n are mapped to segments \(i\). To each segment \(i\) a sections \(j\) with a length \(l_j\) is assigned. The objective is to use as less segments as possible.
The problem is formulated with the objective function \(f_1\):
\begin{equation}
Minimize \quad f_1(\phi_{in}) = - || \sum_{n=1}^{n_n}\phi_{in} || _2, \quad \phi_{in} \in \{0,1\},\psi_{ij} \in \{0,1\}
\end{equation}
The available sections are summarized in the vector \(l_j\):
\begin{equation}
l_{j} \in \{l_1,l_2, ..., l_{j_n}\}
\end{equation}
The problem is subjected to the constraints ...
... each segment \(i\) must be assigned to one or no section \(j\)
\begin{equation} \label{eq:g1}
g_1: \quad \sum_{j=1}^{j_n} \psi_{ij} \leq 1 \quad \forall i
\end{equation}
... each element \(n\) must be assigned to a segment \(i\)
\begin{equation} \label{eq:g2}
g_2: \quad \sum_{i=1}^{i_n} \phi_{in} = 1 \quad \forall n
\end{equation}
... the length of each segment \(i\) must be equal to the length of the assigned section \(j\). \(h_e\) represents the length of a discretized element [\(L = n_n \cdot h_e\)].
\begin{equation}
g_3: \quad \psi_{ij} \cdot l_j = \sum_{n=1}^{n_n} \phi_{in} \cdot h_e \quad \forall i
\end{equation}
... each segment \(i\) must consist of neighboring elements \(n\)
\begin{equation}
g_4: \sum_{n=1}^{n_n} \phi_{in} = \sum_{n=1}^{n_n} c_{in} + \sum_{j=1}^{j_n} \psi_{ij} \quad \forall i
\end{equation}
Whereas \(c_{in}\) is the contiguity matrix:
\begin{equation}
c_{in} = \phi_{in} \cdot \phi_{i,n+1} \quad \forall i,n_n-1;
\end{equation}
The bilinear formulation \(\phi_{in} \cdot \phi_{i,n+1}\) is linearized using:
\begin{equation} \label{eq:c_lin_1}
\quad c_{in} \leq c_{in} \quad \& \quad c_{in} \leq c_{i,n+1} \quad \forall i, n_n-1;
\end{equation}
\begin{equation} \label{eq:c_lin_3}
\quad c_{in} \geq c_{in} + c_{i,n+1} - 1 \quad \quad \forall i, n_n-1;
\end{equation}
The MWE is implemented in Python. Typically the size of variables \(i, n, j\) is much larger, which results in runtimes on a HPC Cluster of multiple days. Nevertheless I have the feeling that this can be significantly reduced, which leads me to my questions:
1) As of now, there are multiple equivalent solutions of \(\phi_{in}\). E.g both of the following represent optimal solutions:
\begin{equation} \label{eq:sol_1}
(1) \quad \phi_{in} =
\begin{pmatrix}
1 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}
\end{equation}
\begin{equation} \label{eq:sol_2}
(2) \quad \phi_{in} =
\begin{pmatrix}
0 & 0 & 1 & 1 & 0 & 0 & 0 \\
1 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}
\end{equation}
The desired solution is (1), where the nonzero entries are distributed along the diagonal. Whenever there is a non zero entry at \(i,n\), the entry \(i,n+1\) or \(i+1,n+1\) is also nonzero. Is there a way to implement this as a constraint?
2) \(f_1\) can be also reformulated as an indicator constraint with the auxiliary variable \(b_i\), delivering different objective values but similar results:
\begin{equation}
\sum_{n=1}^{n_n} \phi_{in} \geq 1 \rightarrow b_i = 1
\end{equation}
\begin{equation}
\sum_{n=1}^{n_n} \phi_{in} < 1 \rightarrow b_i = 0
\end{equation}
\begin{equation}
Minimize \quad f_1(b_i) = || b_i || _1
\end{equation}
This would allow to use the \(L_1\)-norm, in my understanding reducing the problem from MIQCP to MILP. Nevertheless, I observed that this increased the runtime. Is there an explanation for that? Which formulation should be preferred?
3) The linearization of \(c_{in}\) increase the number of constraints significantly. Are there other simplifications that could be applied? Is the linearization always better than the bilinear multiplication?
4) The constraints \(g1\) and \(g2\) could be also expressed as SOS1 constraint or an SOS1 constraint could be added. Would this increase performance?
import numpy as np
import gurobipy as gp
from gurobipy import GRB
def doOptimization(obj_f1):
# Main variables
# -------------------------
L_k = 7.00
l_j = np.array([2.0, 1.0])
h_e = 1.00 # Mesh size, use 0.01 for large example
# Discretization
x = np.linspace(0, L_k, int(L_k/h_e) + 1)
n_n = len(x)
n_e = n_n - 1 # Number of elements
i_n = int(L_k / min(l_j)) + 1 # Number of segments
j_n = l_j.shape[0] # Number of sections
# Declare and initialize model
m = gp.Model('RAP')
# Create variables
# ------------
phi = m.addMVar(shape=(i_n, n_e), vtype=GRB.BINARY, name="phi")
psi = m.addMVar(shape=(i_n, j_n), vtype=GRB.BINARY, name="psi")
# Add constraints
# ------------
# Constraint g1
m.addConstr(psi.sum(axis=1) <= np.ones(i_n))
# Constraint g2
m.addConstr((phi.sum(axis=0)) == np.ones(n_e))
# Constraint g3
m.addConstr(psi @ l_j - phi.sum(axis=1) * h_e == np.zeros(i_n))
# Constraint g4
# Auxiliary variable
c = m.addMVar(shape=(i_n, n_e - 1), vtype=GRB.BINARY, name="c")
# Linearization of phi[i, n] * phi[i, n+1]
m.addConstrs((c[i, n] <= phi[i, n] for i in range(i_n)
for n in range(n_e - 1)))
m.addConstrs((c[i, n] <= phi[i, n+1] for i in range(i_n)
for n in range(n_e - 1)))
m.addConstrs((c[i, n] >= phi[i, n] + phi[i, n+1] - 1 for i in range(i_n)
for n in range(n_e - 1)))
# Create constraint g4
m.addConstr(phi[:, :].sum(axis=1) == (c[:, :].sum(axis=1) + psi[:, :].sum(axis=1)))
# Add objective functions
# ------------
# f1 = ||b(phi)||1 -> min
if obj_f1 == "L1":
aux_f1_b = m.addMVar(shape=(i_n), vtype=GRB.BINARY, name="b_f1")
# Constants
eps = 0.0001
M = n_e + eps
# Add big-M constraints
m.addConstr(phi.sum(axis=1)>= 0 + eps - M * (1 - aux_f1_b), name="bigM_f1_1")
m.addConstr(phi.sum(axis=1)<= 0 + M * aux_f1_b, name="bigM_f1_2")
# Add auxiliary variable
aux_f1_norm = m.addVar()
m.addConstr(aux_f1_norm == gp.norm(aux_f1_b.reshape((i_n, )), 1))
m.setObjective(aux_f1_norm, GRB.MINIMIZE)
# f1 = - ||phi.sum()||2 -> min
elif obj_f1 == "L2":
# Add auxiliary variable
aux_f1_phi_i = m.addMVar((i_n), lb=-float("inf"))
m.addConstr(aux_f1_phi_i == phi.sum(axis=1))
# Add auxiliary variable that computes L_n - Norm
aux_f1_norm = m.addVar()
m.addConstr(aux_f1_norm == gp.norm(aux_f1_phi_i.reshape((i_n, )), 2))
m.setObjective(-aux_f1_norm, GRB.MINIMIZE)
# Optimize model
# ------------
m.optimize()
if __name__ == "__main__":
doOptimization(obj_f1 = "L2")
-
-
Gurobi Staff
Which version of Gurobi are you using? I just tried running your code on my local Mac and with Gurobi v11.0.3 the model was solved in a split second.
Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (mac64[arm] - Darwin 23.6.0 23G93)
CPU model: Apple M3 Pro
Thread count: 11 physical cores, 11 logical processors, using up to 11 threads
Optimize a model with 183 rows, 129 columns and 664 nonzeros
Model fingerprint: 0x5488f601
Model has 1 general constraint
Variable types: 9 continuous, 120 integer (120 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Presolve time: 0.00s
Presolved: 192 rows, 138 columns, 689 nonzeros
Presolved model has 1 quadratic constraint(s)
Presolved model has 8 bilinear constraint(s)
Solving non-convex MIQCP
Variable types: 10 continuous, 128 integer (120 binary)
Root relaxation: objective -8.485281e+00, 72 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -8.48528 0 8 - -8.48528 - - 0s
H 0 0 -3.6055513 -8.48528 135% - 0s
0 0 -5.94205 0 4 -3.60555 -5.94205 64.8% - 0s
0 0 -5.94205 0 22 -3.60555 -5.94205 64.8% - 0s
0 0 -5.94205 0 16 -3.60555 -5.94205 64.8% - 0s
0 0 -5.94205 0 11 -3.60555 -5.94205 64.8% - 0s
0 0 -5.94205 0 11 -3.60555 -5.94205 64.8% - 0s
0 2 -5.94205 0 11 -3.60555 -5.94205 64.8% - 0s
* 171 115 15 -3.6055514 -4.58366 27.1% 15.0 0s
Cutting planes:
Implied bound: 1
Clique: 1
Inf proof: 2
Zero half: 3
RLT: 13
Relax-and-lift: 6
BQP: 2
Explored 2099 nodes (22220 simplex iterations) in 0.10 seconds (0.09 work units)
Thread count was 11 (of 11 available processors)
Solution count 2: -3.60555 -3.60555
No other solutions better than -3.60555
Optimal solution found (tolerance 1.00e-04)
Best objective -3.605551415030e+00, best bound -3.605551415030e+00, gap 0.0000%Could you maybe please share the bigger version of the model which may not solve that quickly? Note that uploading files in the Community Forum is not possible but we discuss an alternative in Posting to the Community Forum.
0 -
-
Dear Jaromił,
Thank you for the reply.
I am using the following setup:
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows Server 2019.0 (17763.2))
CPU model: Intel(R) Xeon(R) Gold 5220R CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 24 physical cores, 48 logical processors, using up to 24 threads1) Please find here the MRE with a larger size of variables \(n, j\), which takes much longer to solve. Simply etting \(h_e = 0.001\) will increase the runtime further.
2) Please find here the full model as .mps file. Please note that within this full model, additional constraints / objectives that for the sake of brevity I did not mention are applied. In my opinion, the MRE is able to represent the complexity and issues of the full formulation. The full model applies the indicator constraint as mentioned in 2), the linearization as mentioned in 3) but not the SOS1-constraints mentioned in 4)
0 -
-
-
Thank you for sharing the models, Benedikt.
The models are indeed very complex. You are minimizing the negative of a 2-norm, i.e., you are trying to maximize a convex function, which is known to be hard. Additionally, each variable of the 2-norm is equal to a very dense sum of binaries. This means that actually your original model has tons of dense products of binaries. This makes the problem even harder.
Is there a way to reduce the density of the constraints by doing some manual presolve? For example, you could use some application oriented knowledge to fix some variables.
Is there an alternative objective you could use which would not introduce the (actually hidden) product terms? Would a sum of absolute values be also applicable?
You could also try to provide a MIP Start, see How do I use MIP start?
You might be interested in our webinar about strong MIP formulations.
Unfortunately, the models are just really hard so there is not much more help I can provide at this point.
Best regards,
Jaromił0 -
Please sign in to leave a comment.
Comments
3 comments