Building Bigger Optimization Model¶

In this document, we will introduce how to decompose the process of building a large optimization model into smaller parts and how to compose them together. This is a common practice in optimization modeling to make the code more readable and maintainable.

Generally speaking, we need two important parts to build an optimization model:

A model object that represents the optimization problem and communicates with the underlying optimizer. Due to the lightweight design philosophy of PyOptInterface, the model object is just a thin wrapper around the optimizer API, and it does not store the optimization problem itself. The model object is responsible for creating and managing variables, constraints, and the objective function, as well as solving the optimization problem.
A dict-like container to store the variables and constraints. We recommend the types.SimpleNamespace object as the container to store the variables and constraints. It is a simple way to store and manipulate the optimization problem in Python. You can use attribute access to get and set the variables and constraints.

import types

container = types.SimpleNamespace()
container.x = 1
container.y = 2

print(container.x, container.y)

1 2

Thus, we can define a class to represent the optimization model and use the container to store the variables and constraints.

import numpy as np
import pyoptinterface as poi
from pyoptinterface import highs

class OptModel:
    def __init__(self):
        self.container = types.SimpleNamespace()
        self.model = highs.Model()

Then we can define small functions to declare the variables and constraints and add them to the model. We take the N-queens problem as an example.

def declare_queen_variables(m:OptModel, N):
    model = m.model
    x = np.empty((N, N), dtype=object)
    for i in range(N):
        for j in range(N):
            x[i, j] = model.add_variable(domain=poi.VariableDomain.Binary)
    container = m.container
    container.x = x
    container.N = N

def add_row_column_constraints(m:OptModel):
    container = m.container
    N = container.N
    x = container.x
    model = m.model

    for i in range(N):
        # Row and column
        model.add_linear_constraint(poi.quicksum(x[i, :]), poi.Eq, 1.0)
        model.add_linear_constraint(poi.quicksum(x[:, i]), poi.Eq, 1.0)

def add_diagonal_constraints(m:OptModel):
    container = m.container
    N = container.N
    x = container.x
    model = m.model

    for i in range(-N+1, N):
        # Diagonal
        model.add_linear_constraint(poi.quicksum(x.diagonal(i)), poi.Leq, 1.0)
        # Anti-diagonal
        model.add_linear_constraint(poi.quicksum(np.fliplr(x).diagonal(i)), poi.Leq, 1.0)

Finally, we can compose these functions to build the optimization model and solve it.

m = OptModel()

N = 8
declare_queen_variables(m, N)
add_row_column_constraints(m)
add_diagonal_constraints(m)

m.model.optimize()

print("Termination status:", m.model.get_model_attribute(poi.ModelAttribute.TerminationStatus))

Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [0e+00, 0e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
42 rows, 64 cols, 252 nonzeros  0s
42 rows, 64 cols, 270 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   42 rows
   64 cols (64 binary, 0 integer, 0 implied int., 0 continuous)
   270 nonzeros
MIP-Timing:      0.0012 - starting analytic centre calculation

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   0               inf                  inf        0      0      6        29     0.0s
 L       0       0         0 100.00%   0               0                  0.00%     1559     28     84       213     0.1s
         1       0         1 100.00%   0               0                  0.00%     1559     28     84       267     0.1s

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.07 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 1
  Nodes             1
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     267 (total)
                    0 (strong br.)
                    184 (separation)
                    54 (heuristics)
Termination status: TerminationStatusCode.OPTIMAL