Mathematical Model Formulation#
Based on the explanations in the previous sections, we now describe how to formulate a mathematical model. Decision variables and placeholders are covered in Defining variables, so here we focus on how to set objectives and constraints.
import jijmodeling as jm
Tip
For convenience, we discuss objectives first and then constraints, but in actual code you can update them in any order.
Setting the objective#
When you create a Problem, setting sense to MAXIMIZE makes it a maximization problem, and setting sense to MINIMIZE makes it a minimization problem.
Right after a Problem is created, the objective is initialized to \(0\), and you add terms to it using the += operator on the Problem object.
The Problem object only accepts scalar Expression objects as objective terms.
If you attempt to add array-typed or dictionary-typed expressions, a type error will be raised.
In JijModeling, you can add terms to the objective, but you cannot overwrite or delete the objective once set.
In particular, += adds a new term and does not replace existing terms.
Consider the following example. First, we set the objective with only the term \(x\).
problem = jm.Problem("Sample")
x = problem.BinaryVar("x")
problem += x
problem
Next, define a new decision variable \(y\) and add it to the objective.
y = problem.BinaryVar("y")
problem += y
problem
You can see that the existing term was not replaced; instead, \(y\) was added and the new objective is \(x + y\). If you might need to remove objective terms later, keep a list of terms in Python and set the objective from that list when needed.
As a more practical example, let’s set the objective for the knapsack problem.
import jijmodeling as jm
@jm.Problem.define("Knapsack Problem", sense=jm.ProblemSense.MAXIMIZE)
def knapsack_problem(problem: jm.DecoratedProblem):
N = problem.Length(description="Number of items")
x = problem.BinaryVar(
shape=(N,), description="$x_i = 1$ if item i is put in the knapsack"
)
v = problem.Float(shape=(N,), description="value of each item")
w = problem.Float(shape=(N,), description="weight of each item")
W = problem.Float(description="maximum weight capacity of the knapsack")
# Set the objective by passing an Expression object to the `+=` operator
problem += jm.sum(v[i] * x[i] for i in N)
# Alternatively, using broadcasting, the following is equivalent
# problem += jm.sum(v * x)
knapsack_problem
Setting constraints#
Constraints are also added using the += operator.
However, when adding constraints, you add Constraint objects created by Problem.Constraint().
Problem.Constraint() takes a name and a constraint expression written with ==, <=, or >= as its required arguments.
Important
The only comparison operators available when building constraints are ==, <=, and >=.
As shown below, operators such as > or <, or any logical operators, are not supported.
problem.Constraint("BAD1", 1 < x) # ERROR! `>` cannot be used!
problem.Constraint("BAD2", (x + y) <= 1 or (y + z) >= 2) # ERROR! logical operators not allowed!
problem.Constraint("BAD2", (x + y) <= 1 | (y + z) >= 2) # ERROR! logical operators not allowed!
Now let’s add a constraint to the knapsack model defined above and complete the model.
@knapsack_problem.update
def _(problem: jm.DecoratedProblem):
N = problem.placeholders["N"]
w = problem.placeholders["w"]
W = problem.placeholders["W"]
x = problem.decision_vars["x"]
problem += problem.Constraint("weight", jm.sum(w[i] * x[i] for i in N) <= W)
knapsack_problem
Always use += when adding constraints
When adding a constraint, always use the += operator.
Simply calling Problem.Constraint() does not add the constraint to the model.
Families of constraints#
In JijModeling, you can add not only single constraints, but also a collection of constraints as a “family”. There are several ways to do this:
Indexed constraints using
domain=or comprehensionsComparison expressions between arrays
To see these patterns, consider the quadratic formulation of the traveling salesman problem (TSP). Let \(d_{i,j}\) be the distance matrix between cities \(i\) and \(j\), and let \(x_{t,i}\) be a binary variable that is \(1\) if city \(i\) is visited at time \(t\). Then we can write:
There are two types of constraints here, and each is defined not as a single constraint but as a family indexed by parameters \(t\) and \(i\).
Indexed constraints#
To define indexed constraints with the Decorator API, provide the second argument to
Problem.Constraint() as a list comprehension or generator expression:
@jm.Problem.define("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
def tsp_decorated(problem: jm.DecoratedProblem):
C = problem.CategoryLabel(description="Labels of Cities")
N = C.count()
x = problem.BinaryVar(
dict_keys=(N, C), description="$x_{t,i} = 1$ if City $i$ is visited at time $t$"
)
d = problem.Float(dict_keys=(C, C), description="distance between cities")
problem += jm.sum(
d[i, j] * x[t, i] * x[(t + 1) % N, j] for t in N for i in C for j in C
)
# Definition using a list comprehension
problem += problem.Constraint(
"one time", [jm.sum(x[t, i] for t in N) == 1 for i in C]
)
# Definition using a generator expression
problem += problem.Constraint(
"one city", (jm.sum(x[t, i] for i in C) == 1 for t in N)
)
tsp_decorated
With the Plain API only, provide a lambda that takes the indexing parameters as the second argument,
and specify the domain using the domain= keyword:
tsp_plain = jm.Problem("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
C = tsp_plain.CategoryLabel("C", description="Labels of Cities")
N = C.count()
x = tsp_plain.BinaryVar(
"x",
dict_keys=(N, C),
description="$x_{t,i} = 1$ if City $i$ is visited at time $t$",
)
d = tsp_plain.Float("d", dict_keys=(C, C), description="distance between cities")
tsp_plain += jm.sum(
jm.product(N, C, C), lambda t, i, j: d[i, j] * x[t, i] * x[(t + 1) % N, j]
)
# Definition using the `domain=` keyword argument
tsp_plain += tsp_plain.Constraint(
"one time", lambda i: jm.sum(N, lambda t: x[t, i]) == 1, domain=C
)
tsp_plain += tsp_plain.Constraint(
"one city", lambda t: jm.sum(C, lambda i: x[t, i]) == 1, domain=N
)
tsp_plain
Array-to-array comparisons#
Another way to define a family of constraints is to use comparison expressions between arrays or sets. As mentioned in Constructing Expressions, comparison expressions also support broadcasting. Specifically, the comparison expressions that can be used to construct constraints are those whose left and right sides are one of the following combinations:
Set and scalar
Arrays of the same shape
TotalDictobjects with the same key set
Using this, we can define the TSP constraints as follows:
@jm.Problem.define("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
def tsp_array_comparison(problem: jm.DecoratedProblem):
N = problem.Natural(description="Number of cities")
x = problem.BinaryVar(
shape=(N, N), description="$x_{t,i} = 1$ if City $i$ is visited at time $t$"
)
d = problem.Float(shape=(N, N), description="distance between cities")
problem += jm.sum(
d[i, j] * x[t, i] * x[(t + 1) % N, j] for t in N for i in N for j in N
)
# Definitions using set-scalar comparison
problem += problem.Constraint("one time", x.sum(axis=0) == 1)
problem += problem.Constraint("one city", x.sum(axis=1) == 1)
tsp_array_comparison
Here, giving an axis=i argument to Expression.sum() or jm.sum() works the same way as numpy.sum(): rather than a single total sum, it returns an array of sums along that axis.
You can also pass multiple axes as a list.
In the one city constraint above, x.sum(axis=1) (0-indexed) sums along the second axis, which corresponds to cities, and produces an array representing the number of cities visited at each time.
If you run type inference, you can see that it is a one-dimensional array.
tsp_array_comparison.infer(tsp_array_comparison.decision_vars["x"].sum(axis=1))
We then compare this one-dimensional array (the number of cities per time) with the scalar value \(1\) to define a constraint family.
The same holds for one time.
In this example the comparison is array-to-scalar, but as mentioned earlier, you can also define constraint families by comparing arrays of the same shape.
Constraints with the same name#
When building complex models, you may want to add constraints with the same name in multiple places. In JijModeling, this is only possible if all of the following conditions are satisfied:
The number of indices for the constraints with the same name must match
Constraints without indices are treated as having \(0\) indices
If any indices overlap, the constraint definitions must match exactly, including description and domain
Condition (2) is checked at compile time into an instance, not at model construction, because the result depends on the instance data. This is because, in principle, index overlap can only be determined at compile time. Consider the following example:
@jm.Problem.define("Possibly Overlapping Constraints")
def problem(problem: jm.DecoratedProblem):
N = problem.Natural(ndim=1)
M = problem.Natural(ndim=1)
n = jm.max(N.max(), M.max()) + 1
x = problem.IntegerVar(shape=n, lower_bound=0, upper_bound=10)
problem += problem.Constraint("constr", [x[i] >= 1 for i in N])
problem += problem.Constraint("constr", [x[i] <= 2 for i in M])
problem
In this example, \(N\) and \(M\) are one-dimensional index sets, and you cannot tell whether they overlap until you supply data. For instance, with the following instance data, the problem compiles into an instance without issues:
instance_ok = problem.eval({"N": [0, 2, 4], "M": [1, 3, 5]})
instance_ok.constraints_df
| equality | type | used_ids | name | subscripts | description | |
|---|---|---|---|---|---|---|
| id | ||||||
| 0 | <=0 | Linear | {0} | constr | [0] | <NA> |
| 1 | <=0 | Linear | {2} | constr | [2] | <NA> |
| 2 | <=0 | Linear | {4} | constr | [4] | <NA> |
| 3 | <=0 | Linear | {1} | constr | [1] | <NA> |
| 4 | <=0 | Linear | {3} | constr | [3] | <NA> |
| 5 | <=0 | Linear | {5} | constr | [5] | <NA> |
On the other hand, if \(N\) and \(M\) overlap, a compile-time error is raised:
try:
instance_ng = problem.eval({"N": [0, 2, 4], "M": [1, 2, 5]})
except jm.ModelingError as e:
print(e)
Traceback (most recent last):
while evaluating problem `Problem(name="Possibly Overlapping Constraints", sense=MINIMIZE, objective=0, constraints={constr: [Constraint(name="constr", , lambda i: x[i] >= 1, domain=set(N)),Constraint(name="constr", , lambda i: x[i] <= 2, domain=set(M)),],})',
defined at File "/tmp/ipykernel_719/316836703.py", line 1, col 2-55
File "/tmp/ipykernel_719/316836703.py", line 1, col 2-55:
1 | @jm.Problem.define("Possibly Overlapping Constraints")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Constraint 'constr' already has subscript 2 with conflicting definion!
old: UserInput(Constraint(name="constr", , lambda i: x[i] >= 1, domain=set(N)))
new: UserInput(Constraint(name="constr", , lambda i: x[i] <= 2, domain=set(M)))
Cases where conflicts are immediately obvious (e.g., mismatched index dimensions, or name collisions for non-indexed constraints) are reported at model construction time rather than at compile time.
try:
@jm.Problem.define("Trivially Conflicting Constraints")
def problem(problem: jm.DecoratedProblem):
x = problem.IntegerVar(lower_bound=0, upper_bound=10)
problem += problem.Constraint("constr", x >= 1)
problem += problem.Constraint("constr", x <= 2)
except jm.ModelingError as e:
print(e)
Traceback (most recent last):
while adding constraint 'constr',
defined at File "/tmp/ipykernel_719/2812696639.py", line 7, col 9-56
File "/tmp/ipykernel_719/2812696639.py", line 7, col 9-56:
7 | problem += problem.Constraint("constr", x <= 2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Constraint 'constr' has conflicting definition!
existing: Constraint(name="constr", sense=GREATER_THAN_EQUAL, left=x, right=1, shape=Scalar(Integer))
defined at: File "/tmp/ipykernel_719/2812696639.py", line 6, col 20-56
new: Constraint(name="constr", sense=LESS_THAN_EQUAL, left=x, right=2, shape=Scalar(Integer))
defined at: File "/tmp/ipykernel_719/2812696639.py", line 7, col 20-56
That said, because constraint-name collisions can only be determined at compile time in some cases, we generally recommend avoiding same-name constraints.