Skip to main content

Introduction of JijModeling

This section introduces the first time users of JijModeling to mathematical modeling with JijModeling.

Mathematical Modeling

Mathematical Model

A mathematical model is a mathematical expression of a problem. For example, let's look at a mathematical model for the problem of "selecting the smallest element from the elements of a one-dimensional list d=[d0,d1,,dN1]\mathbf{d} = [d_0, d_1, \dots, d_{N-1}].

minn=0N1dnxns.t.n=0N1xn=1xn{0,1},n\begin{align} &\min \qquad \sum_{n=0}^{N-1} d_n x_n \\ &\mathrm{s.t.} \qquad \sum_{n=0}^{N-1} x_n = 1 \\ &x_n \in \{ 0, 1\}, \quad \forall n \end{align}

We aim to minimize equation (1) with NN binary variables x=(x0,x1,,xN1)\mathbf{x} = (x_0, x_1, \dots, x_{N-1}) expressed in equation (3). If we only minimize equation (1), then x=0\mathbf{x} = \mathbf{0}. However, the restriction in equation (2) requires that only one of x0,x1,,xN1x_0, x_1, \dots, x_{N-1} takes one and the rest are 0. Therefore, by minimizing equation (1) while satisfying equation (2), only the binary variable corresponding to the smallest element of d\mathbf{d} will be 1. From the above, we can see that this mathematical model can express the problem of "selecting the smallest element from a one-dimensional list.

Mathematical model and Instances

In the above, we have not specified a specific value for d\mathbf{d}. It could be [10, -3, 6.6, 0.1] or some other value. When solving the problem, we must determine the specific value of d\mathbf{d}.

Such a specific value for solving the problem is called an instance. It is helpful to know the word instance because JijModeling separates mathematical models and instances.

In the above example, the mathematical model and instances are related as shown in the table below.

Mathematical modelminndnxn, s.t.nxn=1,xn{0,1}\min \sum_n d_n x_n,~\mathrm{s.t.} \sum_n x_n =1, x_n \in \{0, 1\}
Instanced=[0.5,1,1,1,3]\mathbf{d} = [0.5, 1, -1, 1, 3]

You can read instances as the specific data of the problem to be solved.

Mathematical modeling by JijModeling

Let's start modeling with JijModeling.

As mentioned above, JijModeling separately represents "mathematical models" and "instances. So, first, we build a mathematical model that does not depend on instances. After that, we assign instances to a mathematical model.

First, we create the problem we want to solve as follows

import jijmodeling as jm

# create problem instance
problem = jm.Problem('my_first_problem')

After this, we will add objective functions and constraints to this problem object.

Variables for instances

Let's define variables for the instances needed in the optimization problem (e.g., distances between cities in a traveling salesman problem).
A variable AA can be defined as follows.

import jijmodeling as jm

# Define variable A
A = jm.Placeholder('A')

If you want to define a variable like BijB_{ij}, which is a two-dimensional array, use the following.

# Define a 2-dimensional variable sequence B
B = jm.Placeholder('B', dim=2)

We can declare only that BijB_{ij} is a 2-dimensional list, without specifying the number of elements in BijB_{ij}. By defining it in this way, we can have the flexibility to change the number of elements in BijB_{ij} on the instance side without having to modify this part of the definition.
If you want to pull information on the number of elements from this 2-dimensional array, use shape as follows.

# Define a 2-dimensional variable sequence B
B = jm.Placeholder('B', dim=2)
# Extract the number of rows N in B
N = B.shape[0].
# Extract the number of columns M of B
M = B.shape[1]

You can specify the dimensions of the variables by dim and get information on the number of elements in each dimensional variable from shape.

If you want to specify the number of elements in advance, define the following

# Define the number of rows C, columns D
C = jm.Placeholder('C')
D = jm.Placeholder('D')
# Define 2-dimensional variable column E with rows C and columns D
E = jm.Placeholder('E', shape=(C, D))

Describe the number of elements in a tuple in the shape argument of jm.Placeholder.

Decision variables

The decision variable used for combinatorial optimization is the binary variable x={0,1}x = \{ 0, 1\} can be defined as follows.

import jijmodeling as jm

# Define a binary variable x
x = jm.Binary('x')

This allows us to define a single binary variable, x.
By also using the Placeholder from earlier, we can also define a three-dimensional binary variable sequence yy.

# Define variables F, G, H
F = jm.Placeholder('F')
G = jm.Placeholder('G')
H = jm.Placeholder('H')
# Define a 3-dimensional binary variable sequence y for F x G x H
y = jm.Binary('y', shape=(F, G, H))

Just as you define instance variables like a multidimensional tensor with jm.Placeholder, you can define a multidimensional decision variable tensor with jm.Binary by specifying the number of elements in a tuple for shape.

We can also use an integer variable. You can define an integer variable that satisfies LaUL \leq a \leq U as follows.

# Define lower L and upper U
L = jm.Placeholder('L')
U = jm.Placeholder('U')
# Define an integer variable a
a = jm.LogEncInteger('a', lower=L, upper=U)

Integer variables for multidimensional lists are also supported. To define an integer variable that satisfies LijaijUij(i,j)L_{ij} \leq a_{ij} \leq U_{ij} \leq (\forall i, j), do the following.

# Define the number of elements in the integer variable
I = jm.Placeholder('I')
J = jm.Placeholder('J')
# Define lower L and upper U
L = jm.Placeholder('L', shape=(I, J))
U = jm.Placeholder('U', shape=(I, J))
# Define an integer variable a
a = jm.LogEncInteger('a', lower=L, upper=U, shape=(I, J))

For such a definition, Lij,aij,UijL_{ij}, a_{ij}, U_{ij} must be multidimensional array of the same dimension and same shape.

If all integer variables in a multidimensional array have the same upper uu and lower \ell, such as biju\ell \leq b_{ij} \leq u, you can define the following.

ell = jm.Placeholder('ell')
u = jm.Placeholder('u')
I = jm.Placeholder('I')
J = jm.Placeholder('J')
b = jm.LogEncInteger('b', lower=ell, upper=u, shape=(I, J))

Variables for indexing

When they are represented in the form of (multidimensional) arrays, such as the variable dnd_n for an instance or the decision variable xnx_n for optimization, we need an index that specifies how many elements they are. Here is how this is defined.
The index nn used for a one-dimensional binary variable xnx_n or an instance variable dnd_n is defined as follows using jm.Element.

# Define variable N
N = jm.Placeholder('N')
# define subscript n with 0≤n<N
n = jm.Element('n', (0, N))

where the (0, N) tuple represents the possible range of the index n, as in 0n<N0 \leq n < N.

caution

Note that NN is not included when written as (0, N). (0, N) means {0,1,...,N1}\{0, 1, ..., N-1\}.

Any set can be used for the part of jm.Element that specifies the range. For example

# Define a variable M
M = jm.Placeholder('M', dim=1)
# define an index m satisfying m ∈ M
m = jm.Element('m', M)

and then write M={2,3,5,7,11}M = \{ 2, 3, 5, 7, 11\} as an instance, you have defined m{2,3,5,7,11}m \in \{ 2, 3, 5, 7, 11\}.
In addition, a two-dimensional list or similar can be used for the part that specifies the range. For example

# Define a 2-dimensional variable E
E = jm.Placeholder('E', dim=2)
# define an index e satisfying e[0] ∈ E[0], e[1] ∈ E[1].
e = jm.Element('e', E)

and give E=[[1,2],[3,4],[5,6]]E = [[1, 2], [3, 4], [5, 6]] as instances, so that e[0] is 1, 3, 5 and e[1] will act as an index that takes values only for 2, 4, and 6. In this way, for example, for a set of edges EE in a graph, (uv)E(uv) \in E, we can define the vertex (uv)(uv) that forms that edge.

Basic Summation operation

Here is how to implement sums of sums appearing in a mathematical model using JijModeling.

Let's start with the basic

n=0N1xn\sum_{n=0}^{N-1} x_n

is an implementation of an operation such as It is written as follows

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N))
n = jm.Element('i', (0, N))
# Calculate the sum
sum_x = jm.Sum(n, x[n])

The desired sum operation can be implemented by entering the subscripts to be summed as the first argument of jm.Sum and the formula to be summed as the second argument.

jm.Sum(nsubscript,x[n]operand)\mathrm{jm.Sum(\overbrace{n}^{subscript}, \underbrace{x[n]}_{operand})}

Substituting the index defined in jm.Element for the index to be summed, the sum of all possible ranges for that index is calculated. In the present case, the operation is to find the sum of n=0,1,,N1n = 0, 1, \dots, N-1.
Such a simple sum can be written in abbreviated form as follows

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N))
# Simplified notation for sum calculation
sum_x = x[:]

The sum operation in JijModeling can also be performed on sets. In the combinatorial optimization problem.

vVvxv\sum_{v \in V} v x_v

mathematical models often appear that take the sum over some set VV, such as $$ \sum_{v \in V} v x_v $$. This implementation is as follows.

# Define variables to be used
V = jm.Placeholder('V', dim=1)
num_V = V.shape[0].
x = jm.Binary('x', shape=(num_V))
v = jm.Element('v', V)
# Compute the sum over a set V
sum_v = jm.Sum(v, v*x[v])

By defining v as an element of a set in jm.Element and writing jm.Sum with v as its first argument, it is easy to write an implementation of summing over a set V. All that remains is to prepare concrete values of that set V as instances.

Multiple Sum Operations

To compute multiple sums, one may write them in a formula as follows.

i,jQijxij\sum_{i, j} Q_{ij} x_{ij}

This can be implemented in JijModeling as follows.

# Define variables to be used
Q = jm.Placeholder('Q', dim=2)
I = Q.shape[0].
J = Q.shape[1]
x = jm.Binary('x', shape=(I, J))
i = jm.Element('i', (0, I))
j = jm.Element('j', (0, J))
# Calculate the sum over the two indices i, j
sum_ij = jm.Sum([i, j], Q[i, j]*x[i, j])

If there are multiple sums, multiple jm.Sum can be omitted by making the subscripts a list [subscript1, subscript2, ...]. Of course, this can be written as i,j=ij\sum_{i, j} = \sum_{i} \sum_{j}, as can be seen from

sum_ij = jm.Sum(i, jm.Sum(j, Q[i, j]*x[i, j]))

Sum Operation with condition

In addition to summing over the entire range of possible indexes, you may want to sum only if the indexes satisfy certain conditions, as in the following.

i<Uxi\sum_{i<U} x_i

This can be implemented using JijModeling as follows

# Define variables to be used
I = jm.Placeholder('I')
x = jm.Binary('x', shape=(I))
i = jm.Element('i', (0, I))
U = jm.Placeholder('U')
# Calculate sum only for parts satisfying i<U
sum_i = jm.Sum((i, i<U), x[i])

In the part of jm.Sum specifying the index, use a tuple to specify (index, condition) like this.

jm.Sum((iindex,i<Ucondition),x[i]operand)\mathrm{jm.Sum((\underbrace{i}_{index}, \overbrace{i<U}^{condition}), \underbrace{x[i]}_{operand})}

The conditional expressions satisfied by the subscripts can be <, <=, >=, >, ==, ! = can be used for the conditional expression that the subscript satisfies.

If you want to impose multiple conditions on a subscript, for example.

i<Ui!=Ndixi\sum_{\substack{i < U \\ i!=N}} d_i x_i

Suppose you want to implement a formula such as This can be implemented as follows

# Define variables to be used
d = jm.Placeholder('d', dim=1)
I = d.shape[0].
x = jm.Binary('x', shape=(I))
U = jm.Placeholder('U')
N = jm.Placeholder('N')
i = jm.Element('i', (0, I))
# Calculate sum only for the part satisfying i<U and i≠N
sum_i = jm.Sum((i, (i<U)&(i!=N)), d[i]*x[i])

When multiple conditions are imposed, the operator & for logical AND and | for logical OR can be used to chain the conditions. Of course, more than three conditions can be imposed.

jm.Sum((isubscript,(i<U)condition1&logicaloperator(i!=N)condition2),d[i]x[i]operand)\mathrm{jm.Sum((\underbrace{i}_{subscript}, \overbrace{(i<U)}^{condition 1} \underbrace{\&}_{logical operator} \overbrace{(i!=N)}^{condition 2}), \underbrace{d[i]*x[i ]}_{operand})}

If there is a condition on the subscripts in a sum operation on multiple subscripts, e.g..

i>Li!=Nj<iRijxij\sum_{\substack{i>L \\ i!=N}} \sum_{j<i} R_{ij} x_{ij}

A formula such as the following can be implemented as follows

# Define variables to be used
R = jm.Placeholder('R', dim=2)
I = R.shape[0]
J = R.shape[1]
x = jm.Binary('x', shape=(I, J))
i = jm.Element('i', (0, I))
j = jm.Element('j', (0, J))
N = jm.Placeholder('N')
L = jm.Placeholder('L')
# Calculate sum only for the part satisfying i>L and i!=N and j<i
sum_ij = jm.Sum([(i, (i>L)&(i!=N)), (j, j<i)], R[i, j]*x[i, j])

The part of jm.Sum specifying the indexes should be replaced with [[(index 1, condition of index 1), (index 2, condition of index 2), ...]] By making it look like [(index1, condition of index1), (index2, condition of index2), ...]', multiple sum operations can be written with conditions attached to each index.

caution

[(j, j<i), (i, (i>L)&(i!=N))] will occur an error because ii is not yet defined in the j<ij<i part. This can be written in a formula like i>Li!=N(j<i)\sum_{\substack{i>L \\ i!=N}} \left( \sum_{j<i} \cdots \right), but j<i(i>Li!=N)\sum_{j<i} \left( \sum_{\substack{i>L \\ i!=N}} \cdots \right) corresponds to the fact that it cannot be written as Please note the order in which conditions are imposed on subscripts in multiple sums.

Notes on Replacing Sums by Variables

For example, consider the following implementation of the formula.

H=i{(j=0j>iN1xij)(j=0j>iN1xij1)}H = \sum_i \left\{ \left( \sum_{\substack{j=0 \\ j>i}}^{N-1} x_{ij} \right) \left( \sum_{\substack{j=0 \\ j>i}}^{N-1} x_{ij} -1 \right) \right\}

Since j>ixij\sum_{j>i} x_{ij} appears twice, the following can be implemented by temporarily replacing it with an appropriate variable.

# Define the variable to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N, N))
i = jm.Element('i', (0, N))
j = jm.Element('j', (0, N))
# Calculate sum only for j satisfying j>i
sum_j = jm.Sum((j, j>i), x[i, j])
# Compute sum for i
H = jm.Sum(i, sum_j*(sum_j-1))

Note that if jm.Sum is replaced by a variable, it will not be written as sum_j[i] in the following.

Adding constraints

Equality Constraints

Let us insert the following constraint condition into the problem.

ixi=1\sum_i x_i = 1

By using an abbreviated description of the sum, we can easily implement the following.

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N))
# Create a problem
problem = jm.Problem('constraint_example_01')
# Add constraint: Σ x = 1
problem += jm.Constraint('equality', x[:]==1)

Set a constraint with jm.Constraint. The first argument is a string that is the name of that constraint, and the second argument is an equation with == representing an equality constraint.
The constraint so created is then added to problem using the += operator.

Inequality Constraints

Using jm.Constraint and <= from earlier, you can also easily implement inequality constraints. For example, consider implementing the following equation.

i=0N1xi<i=0N1yi+U\sum_{i=0}^{N-1} x_i < \sum_{i=0}^{N-1} y_i + U

The implementation using JijModeling is as follows.

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N))
y = jm.Binary('y', shape=(N))
U = jm.Placeholder('U')
# Create problem
problem = jm.Problem('constraint_example_02')
# Add constraint: Σ x < Σ y + U
problem += jm.Constraint('inequality', x[:]-y[:]<=U-1)
note

We can only use >=, <= as a constraint in JijModeling.

Using forall to define constraints in a batch

There are many combinatorial optimization problems where there are NN constraints for subscripts as follows.

jxij=1(i{0,1,,N1})\sum_j x_{ij} = 1 \quad (\forall i \in \{0, 1, \dots, N-1\})

Such constraints can also be implemented cleanly in JijModeling.

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N, N))
i = jm.Element('i', (0, N))
# Create problem
problem = jm.Problem('constraint_example_03')
# Add constraint: Σ_j x_{ij} = 1 (∀ i)
problem += jm.Constraint('forall', x[i, :]==1, forall=i)

By entering a subscript defined in jm.Element as the forall argument of jm.Constraint, you can set constraints on all possible ranges for that subscript at once.

Multiple forall

As the number of subscripts increases, constraint conditions may be imposed on multiple \forall symbols as follows.

jxi,j,k=1(i{0,1,,N1},k{0,1,,N1})\sum_j x_{i, j, k} = 1 \quad (\forall i \in \{0, 1, \dots, N-1\}, k \in \{ 0, 1, \dots, N-1\})

Similar to the way jm.Sum sums over multiple subscripts, [subscript 1, subscript 2, ...]], this can be implemented by listing them as follows.

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N, N, N))
i = jm.Element('i', (0, N))
k = jm.Element('k', (0, N))
# Create problem
problem = jm.Problem('constraint_example_04')
# Add constraint: Σ_j x_{ijk} = 1 (∀ i, k)
problem += jm.Constraint('multi_forall', x[i, :, k]==1, forall=[i, k])

Attaching conditions to forall

Unlike before, we may wish to impose a constraint on all subscripts satisfying a certain condition as follows.

ixi,j,k=Ak(j,kj>k)\sum_i x_{i, j, k} = A_k \quad (\forall j, k \vert j > k)

As with jm.Sum, this can also be implemented by entering a tuple like (subscript, condition).

# Define variables to be used
A = jm.Placeholder('A', dim=1)
N = A.shape[0].
x = jm.Binary('x', shape=(N, N, N))
i = jm.Element('i', (0, N))
j = jm.Element('j', (0, N))
k = jm.Element('k', (0, N))
# Create problem
problem = jm.Problem('constraint_example_05')
# Add constraint: Σ_i x_{ijk} = A_k (∀ j, k | j > k)
problem += jm.Constraint('multi_forall_condition', x[:, j, k]==A[k], forall=[k, (j, j>k)])

Also, if you want to set conditions on forall for multiple indexes, you can use [[(index1, condition for index1), (index2, condition for index2), ...]]. It is also possible to impose a condition on each index by doing. Multiple conditions on a single index can be written using the AND operator & and the OR operator |.

note

A statement like [(j, j>k), k] will occur an error because kk is not yet defined in the j>kj>k part. Please be careful about the order of conditions in forall for multiple subscripts.

Variable Fixation by Obvious Constraints

When mathematically modeling a combinatorial optimization problem, it is sometimes possible to fix variables based on trivial constraints. In JijModeling, this can be easily implemented by using jm.Constraint. For example, to fix x0,0=0x_{0, 0}=0 in a binary variable sequence xi,jx_{i, j}, use the following

# Define the variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N, N))
# Create problem
problem = jm.Problem('constraint_example_06')
# Add an obvious constraint: x_{0, 0} = 1
problem += jm.Constraint('single_fix', x[0, 0]==1)

When such an obvious constraint is added, x0,0x_{0, 0} is not passed to the solver as a decision variable, but is translated into an instance fixed at 0.
However, the obvious constraints must satisfy the following requirements.

  • It is an equational constraint
  • The left-hand side must be a decision variable only and must not contain other operations such as jm.Sum.
  • The right-hand side must not contain a decision variable

It is also possible to set them all at once using forall. For example,

x0,j=0(j)x_{0, j} = 0 \quad (\forall j)

would be implemented as follows.

# Define variables to be used
N = jm.Placeholder('N')
x = jm.Binary('x', shape=(N, N))
j = jm.Element('j', (0, N))
# Create problem
problem = jm.Problem('constraint_example_07')
# Add an obvious constraint: x_{0, j} = 0 (∀ j)
problem += jm.Constraint('forall_fix', x[0, j]==1, forall=j)

Check the implementation

If you are using Jupyter Notebook, you can check the mathematical model implemented in JijModeling while implementing. The following figure shows the actual implementation of equations (1)~(3) on Jupyter Notebook.

By using .set_latex in the variable definition, you can set the formula display on Jupyter Notebook. If not set using .set_latex, the variable N defined in d.shape[0] will be displayed as dmathrmshape(0)d_mathrm{shape(0)} by default, as follows.

Solving with JijZept

JijZept subscribers can solve optimization problems by sending the mathematical model implemented in JijModeling to the JijZept server. In doing so, the numerical data of the actual instance to be solved is also passed at the same time.

For example, the solution to the previous implementation is as follows

import jijzept as jz

# Prepare the instance
instance_data = {'d': [1, 10, -5, 0.45, -0.1]}
# Set up sampler
sampler = jz.JijSASampler(token='*** your API key ***', url='https://api.jijzept.com')
# Calculate using the sampler
response = sampler.sample_model(problem, instance_data, {'onehot': 5.0})

If a constraint is inserted into problem using jm.Constraint, the third argument to sampler.sample_model is the undetermined multiplier of that constraint, of type {'constraint name': undetermined multiplier number} in the dictionary. If there is no constraint, specify an empty dictionary.

JijZept uses the .sample_model of the Sampler class to submit the problem to be solved. See the JijZept documentation for more details.

Solving with OpenJij in the local environment

You can also solve optimization problems implemented in JijModeling in your local environment by converting the JijModeling implementation to PyQUBO, which can then be used in OpenJij PyQUBO and other solvers.

The conversion to PyQUBO uses the information of a previously prepared instance.

import openjij as oj

# Convert to PyQUBO
pyq_obj = problem.to_pyqubo(instance_data, {})
# Convert for OpenJij using PyQUBO functionality
bqm = pyq_obj.compile().to_bqm(feed_dict={'onehot': 5.0})
# Prepare SA sampler
sampler = oj.SASampler()
# Solve using SA
response = sampler.sample(bqm)

Decoding the computation, viewing the results

Both JijZept and OpenJij return the result in the form of a response. However, this response has a structure that is difficult to see when checking the obtained solution. Therefore, JijModeling decodes the response to analyze and visualize whether the obtained results satisfy the constraint conditions.

The decoding is done as follows

decoded = problem.decode(response, instance_data, {})

Decoding is done by passing the result, response, and the instance, instance_data, to .decode of the Problem class.

Also, to check the state of the solution (decision variable sequence), the value of the objective function, and the constraint breaking from the decode obtained in this way, do the following.

print(decoded.solutions) # <- check solution state
print(decoded.objectives) # <- check the value of the objective function
print(decoded.constraint_violations) # <- check for constraint violation

Running this, for example, yields the following output for an optimization problem using the mathematical model described in equations (1)~(3) above.

[{'x': array([0., 0., 1., 0., 0.])}]
[-5.]]
[{'onehot': 0.0}]
  • decoded.solutions gives the state of the resulting computation. In this case, the decision variable sequence x is used, and what kind of {0, 1} number sequence it became is converted to numpy.ndarray and output.
  • In decoded.objects, the value of the objective function obtained from the result of the obtained calculation is obtained. In this case, the value -12.0 is obtained.
  • In decoded.constraint_violations, you can see the constraint violations from the obtained calculation results. Here only one constraint, onehot, is shown. If multiple constraints are added, the degree to which each constraint is violated is displayed.

It is also possible to extract only those solutions for which all constraints are satisfied (feasible solutions). To do so, use the following

print(decoded.feasibles().solution) # <- check the status of feasible solutions
print(decoded.feasibles().objective) # <- check the value of the objective function of the feasible solution