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 $\mathbf{d} = [d_0, d_1, \dots, d_{N-1}]$.

\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 $N$ binary variables $\mathbf{x} = (x_0, x_1, \dots, x_{N-1})$ expressed in equation (3). If we only minimize equation (1), then $\mathbf{x} = \mathbf{0}$. However, the restriction in equation (2) requires that only one of $x_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 $\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 $\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 $\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 model$\min \sum_n d_n x_n,~\mathrm{s.t.} \sum_n x_n =1, x_n \in \{0, 1\}$
Instance$\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 $A$ can be defined as follows.

import jijmodeling as jm# Define variable AA = jm.Placeholder('A')

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

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

We can declare only that $B_{ij}$ is a 2-dimensional list, without specifying the number of elements in $B_{ij}$. By defining it in this way, we can have the flexibility to change the number of elements in $B_{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 BB = jm.Placeholder('B', dim=2)# Extract the number of rows N in BN = B.shape[0].# Extract the number of columns M of BM = 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 DC = jm.Placeholder('C')D = jm.Placeholder('D')# Define 2-dimensional variable column E with rows C and columns DE = 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\}$ can be defined as follows.

import jijmodeling as jm# Define a binary variable xx = 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 $y$.

# Define variables F, G, HF = jm.Placeholder('F')G = jm.Placeholder('G')H = jm.Placeholder('H')# Define a 3-dimensional binary variable sequence y for F x G x Hy = 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 $L \leq a \leq U$ as follows.

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

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

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

For such a definition, $L_{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 $u$ and lower $\ell$, such as $\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 $d_n$ for an instance or the decision variable $x_n$ for optimization, we need an index that specifies how many elements they are. Here is how this is defined.
The index $n$ used for a one-dimensional binary variable $x_n$ or an instance variable $d_n$ is defined as follows using jm.Element.

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

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

caution

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

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

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

and then write $M = \{ 2, 3, 5, 7, 11\}$ as an instance, you have defined $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 EE = 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]]$ 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 $E$ in a graph, $(uv) \in E$, we can define the vertex $(uv)$ that forms that edge.

Basic Summation operation​

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

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

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

# Define variables to be usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N))n = jm.Element('i', (0, N))# Calculate the sumsum_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.

$\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, \dots, N-1$.
Such a simple sum can be written in abbreviated form as follows

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

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

$\sum_{v \in V} v x_v$

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

# Define variables to be usedV = 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 Vsum_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.

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

This can be implemented in JijModeling as follows.

# Define variables to be usedQ = 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, jsum_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 $\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.

$\sum_{i

This can be implemented using JijModeling as follows

# Define variables to be usedI = 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<Usum_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.

$\mathrm{jm.Sum((\underbrace{i}_{index}, \overbrace{i

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.

$\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 usedd = 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≠Nsum_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.

$\mathrm{jm.Sum((\underbrace{i}_{subscript}, \overbrace{(i

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

$\sum_{\substack{i>L \\ i!=N}} \sum_{j

A formula such as the following can be implemented as follows

# Define variables to be usedR = 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<isum_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 $i$ is not yet defined in the $j part. This can be written in a formula like $\sum_{\substack{i>L \\ i!=N}} \left( \sum_{j, but $\sum_{jL \\ 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 = \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 $\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 usedN = 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>isum_j = jm.Sum((j, j>i), x[i, j])# Compute sum for iH = 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.

Equality Constraints​

Let us insert the following constraint condition into the problem.

$\sum_i x_i = 1$

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

# Define variables to be usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N))# Create a problemproblem = jm.Problem('constraint_example_01')# Add constraint: Σ x = 1problem += 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.

$\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 usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N))y = jm.Binary('y', shape=(N))U = jm.Placeholder('U')# Create problemproblem = jm.Problem('constraint_example_02')# Add constraint: Σ x < Σ y + Uproblem += 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 $N$ constraints for subscripts as follows.

$\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 usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N, N))i = jm.Element('i', (0, N))# Create problemproblem = 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.

$\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 usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N, N, N))i = jm.Element('i', (0, N))k = jm.Element('k', (0, N))# Create problemproblem = 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.

$\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 usedA = 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 problemproblem = 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 $k$ is not yet defined in the $j>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 $x_{0, 0}=0$ in a binary variable sequence $x_{i, j}$, use the following

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

When such an obvious constraint is added, $x_{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,

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

would be implemented as follows.

# Define variables to be usedN = jm.Placeholder('N')x = jm.Binary('x', shape=(N, N))j = jm.Element('j', (0, N))# Create problemproblem = 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 $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 instanceinstance_data = {'d': [1, 10, -5, 0.45, -0.1]}# Set up samplersampler = jz.JijSASampler(token='*** your API key ***', url='https://api.jijzept.com')# Calculate using the samplerresponse = 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 PyQUBOpyq_obj = problem.to_pyqubo(instance_data, {})# Convert for OpenJij using PyQUBO functionalitybqm = pyq_obj.compile().to_bqm(feed_dict={'onehot': 5.0})# Prepare SA samplersampler = oj.SASampler()# Solve using SAresponse = 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 stateprint(decoded.objectives) # <- check the value of the objective functionprint(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 solutionsprint(decoded.feasibles().objective) # <- check the value of the objective function of the feasible solution