Advanced Usage of JijModeling

Advanced Usage of JijModeling#

🌐 日本語 | English

Overview#

This document explains advanced ways of using JijModeling together with JijZept Tools.

Unlike many mathematical-optimization modelers, JijModeling keeps the structure of the input expressions as-is. This makes it possible to transform the expression structure freely, and JijZept Tools implements techniques that take advantage of this property.

The following sections introduce techniques provided by JijZept Tools.

This document covers advanced usage of JijModeling. For the basics of JijModeling, see here.

Converting a problem with no feasible solution into an always-feasible problem#

In mathematical optimization, when the target problem has no feasible solution, it is generally difficult to identify which constraint is to blame.

In such situations, introducing slack variables (auxiliary variables that absorb constraint violations) converts a problem with no feasible solution into an always-feasible problem — and by looking at the values of the slack variables, you can get a good idea of which constraints are causing the trouble.

Let’s look at a concrete example.

Consider the following problem formulated with JijModeling.

import jijmodeling as jm

# Initialize the problem
problem = jm.Problem("Binary_LP")

# Define the parameters
A = jm.Placeholder("A", ndim=2, description="constraint matrix")
b = jm.Placeholder("b", ndim=1, description="right-hand side of the constraints")
c = jm.Placeholder("c", ndim=1, description="cost coefficient vector")

# Define the indices
m = A.len_at(0)
m.set_latex("M")  # number of rows
i = jm.Element("i", belong_to=(0, m), description="row index of the constraints")
n = A.len_at(1)
n.set_latex("N")  # number of columns
j = jm.Element("j", belong_to=(0, n), description="column index of the variables")

# Define the binary variables
x = jm.BinaryVar("x", shape=(n,), description="binary variable for select/deselect")

# Define the objective
problem += jm.sum(j, c[j] * x[j])  # minimize c^T x

# Define the constraints
problem += jm.Constraint(
    "linear_constraint", jm.sum(j, A[i, j] * x[j]) >= b[i], forall=i
)
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{Binary\_LP} & & \\& & \min \quad \displaystyle \sum_{j = 0}^{N - 1} c_{j} \cdot x_{j} & \\\text{{s.t.}} & & & \\ & \text{linear\_constraint} & \displaystyle \sum_{j = 0}^{N - 1} A_{i, j} \cdot x_{j} \geq b_{i} & \forall i \in \left\{0,\ldots,M - 1\right\} \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}& \text{binary variable for select/deselect}\\\end{array}\end{split}\]

Define the instance data as follows and try to solve the problem.

import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

instance_data = {
    "A": np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
    "b": np.array([10, 1, 20]),
    "c": np.array([7, 8, 9]),
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(problem)

# Solve the problem with OMMXPySCIPOptAdapter
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)

---------------------------------------------------------------------------
InfeasibleDetected                        Traceback (most recent call last)
Cell In[2], line 13
     10 ommx_instance = interpreter.eval_problem(problem)
     12 # Solve the problem with OMMXPySCIPOptAdapter
---> 13 solution = OMMXPySCIPOptAdapter.solve(ommx_instance)

File ~/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:163, in OMMXPySCIPOptAdapter.solve(cls, ommx_instance, use_sos1, initial_state)
    161 model = adapter.solver_input
    162 model.optimize()
--> 163 return adapter.decode(model)

File ~/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:213, in OMMXPySCIPOptAdapter.decode(self, data)
    171 """Convert optimized pyscipopt.Model and ommx.v1.Instance to ommx.v1.Solution.
    172 
    173 This method is intended to be used if the model has been acquired with
   (...)
    208 
    209 """
    211 # TODO: Add the feature to store dual variables in `solution`.
--> 213 state = self.decode_to_state(data)
    214 solution = self.instance.evaluate(state)
    216 if (
    217     data.getStatus() == "optimal"
    218 ):  # pyscipopt does not appear to have an enum or constant for this

File ~/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:259, in OMMXPySCIPOptAdapter.decode_to_state(self, data)
    254     raise OMMXPySCIPOptAdapterError(
    255         "The model may not be optimized. [status: unknown]"
    256     )
    258 if data.getStatus() == "infeasible":
--> 259     raise InfeasibleDetected("Model was infeasible")
    261 if data.getStatus() == "unbounded":
    262     raise UnboundedDetected("Model was unbounded")

InfeasibleDetected: Model was infeasible

We got an infeasibility error. As it stands, we cannot tell which constraint is at fault.

Now, let’s use JijZept Tools to introduce slack variables.

A single line of code is all it takes.

from jijzepttools.modeling.replace.problem import (
    generate_always_feasible_problem_forall,
)

replaced_problem = generate_always_feasible_problem_forall(problem)
replaced_problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{Binary\_LP} & & \\& & \min \quad \displaystyle \sum_{j = 0}^{N - 1} c_{j} \cdot x_{j} + bigM \cdot \sum_{i = 0}^{M - 1} \left(v0_{i} + w0_{i}\right) & \\\text{{s.t.}} & & & \\ & \text{linear\_constraint} & \displaystyle \sum_{j = 0}^{N - 1} A_{i, j} \cdot x_{j} + v0_{i} - w0_{i} \geq b_{i} & \forall i \in \left\{0,\ldots,M - 1\right\} \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}& \text{binary variable for select/deselect}\\& v0 & 1\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }100000000 & \\& w0 & 1\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }100000000 & \\\end{array}\end{split}\]

In addition to the original problem, the slack variables v0, w0 have been introduced. The subscripts attached to the slack variables are determined automatically by analyzing the subscripts of the original problem.

The objective function now also contains v0 + w0, so the optimization drives the slack variables toward small values. Let’s solve this problem in the same way as before.

The bigM added to the objective must be defined separately in the instance data.

import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

instance_data = {
    "A": np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
    "b": np.array([10, 1, 20]),
    "c": np.array([7, 8, 9]),
    "bigM": 1000,
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(replaced_problem)

# Solve the problem with OMMXPySCIPOptAdapter
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)

solution.decision_variables_df
kind lower upper name subscripts description substituted_value value
id
0 Binary 0.0 1.0 x [0] <NA> <NA> 1.0
1 Binary 0.0 1.0 x [1] <NA> <NA> 1.0
2 Binary 0.0 1.0 x [2] <NA> <NA> 1.0
3 Continuous 0.0 100000000.0 v0 [0] <NA> <NA> 4.0
4 Continuous 0.0 100000000.0 v0 [1] <NA> <NA> 0.0
5 Continuous 0.0 100000000.0 v0 [2] <NA> <NA> 0.0
6 Continuous 0.0 100000000.0 w0 [0] <NA> <NA> 0.0
7 Continuous 0.0 100000000.0 w0 [1] <NA> <NA> 0.0
8 Continuous 0.0 100000000.0 w0 [2] <NA> <NA> 0.0

Looking at this, we can see that element 0 of the slack variable v0 has a value greater than 0.

Even though the optimization tries to keep the slack variables as small as possible, element 0 of v0 is still positive. From this we can infer that, in the original problem, the constraint corresponding to i=0 of linear_constraint is the cause of the infeasibility.

As a test, let’s decrease the value of b[0] to relax constraint 0.

import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

instance_data = {
    "A": np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
    "b": np.array([1, 1, 20]),
    "c": np.array([7, 8, 9]),
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(problem)

# Solve the problem with OMMXPySCIPOptAdapter
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)

solution.decision_variables_df
kind lower upper name subscripts description substituted_value value
id
0 Binary 0.0 1.0 x [0] <NA> <NA> 1.0
1 Binary 0.0 1.0 x [1] <NA> <NA> 1.0
2 Binary 0.0 1.0 x [2] <NA> <NA> 1.0

This time we obtained a feasible solution. If we also solve the slack-variable version of the problem in the same way, we can see that the slack variables are now 0.

import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

instance_data = {
    "A": np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
    "b": np.array([1, 1, 20]),
    "c": np.array([7, 8, 9]),
    "bigM": 1000,
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(replaced_problem)

# Solve the problem with OMMXPySCIPOptAdapter
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
solution.decision_variables_df
kind lower upper name subscripts description substituted_value value
id
0 Binary 0.0 1.0 x [0] <NA> <NA> 1.0
1 Binary 0.0 1.0 x [1] <NA> <NA> 1.0
2 Binary 0.0 1.0 x [2] <NA> <NA> 1.0
3 Continuous 0.0 100000000.0 v0 [0] <NA> <NA> 0.0
4 Continuous 0.0 100000000.0 v0 [1] <NA> <NA> 0.0
5 Continuous 0.0 100000000.0 v0 [2] <NA> <NA> 0.0
6 Continuous 0.0 100000000.0 w0 [0] <NA> <NA> 0.0
7 Continuous 0.0 100000000.0 w0 [1] <NA> <NA> 0.0
8 Continuous 0.0 100000000.0 w0 [2] <NA> <NA> 0.0

In this way, introducing slack variables converts a problem with no feasible solution into an always-feasible problem, and inspecting the slack-variable values lets you identify which constraints are causing the trouble.

Contents