Petri Net Modeling for Ising Model Formulation in Quantum Annealing

Quantum annealing is an emerging new platform for combinatorial optimization, requiring an Ising model formulation for optimization problems. The formulation can be an essential obstacle to the permeation of this innovation into broad areas of everyday life. Our research is aimed at the proposal of a Petri net modeling approach for an Ising model formulation. Although the proposed method requires users to model their optimization problems with Petri nets, this process can be carried out in a relatively straightforward manner if we know the target problem and the simple Petri net modeling rules. With our method, the constraints and objective functions in the target optimization problems are represented as fundamental characteristics of Petri net models, extracted systematically from Petri net models, and then converted into binary quadratic nets, equivalent to Ising models. The proposed method can drastically reduce the difficulty of the Ising model formulation.

This study extends our MILP version to quantum annealing models, where quadratic functions of binary variables need to be formulated. As far as we know, this is the first study to apply Petri net modeling to the quantum annealing. In this paper, we introduce a new notation called binary quadratic nets, which denote Ising or QUBO models, respectively. Our method incrementally targets binary quadratic nets from problem-domain Petri net models that represent target optimization problems. This paper is organized as follows: Section 2 summarizes the basic knowledge on quantum annealing and Petri nets. Section 3 introduces a new class of Petri nets, binary quadratic nets, for representing Ising or QUBO models. Section 4 proposes a method for formulating Ising or QUBO models from problem-domain Petri nets and presents some examples. Finally, in Section 5, we present some concluding remarks and areas of future tasks.

Preliminaries
This section provides the basic definitions and notations of Ising models and Petri nets.

Ising and QUBO Models
Quantum annealing is a metaheuristic for solving combinatorial optimization problems, where we try to find a value +1 or −1 for each of the Ising variables s = (s 1 , s 2 , ..., s N ) such that the following Hamiltonian of the Ising models gives the lowest energy state: where h i is the magnetic field coefficient at site s i , and J i,j is the interaction coefficient between s i and s j . Ising variables correspond to discrete variables in metaheuristics. The lowest energy state, called the ground state of H, provides the optimal solution. Figure 2 depicts a two-dimensional array Ising model. A circle represents a spin and is connected to four neighbor spins. The arrows ↑ and ↓ represent the up spin (+1) and down spin (-1), respectively.
Up own Figure 1: 2D Array Ising Model A Hamiltonian may include an objective function, where the lower energy on the terms of the function contributes to a smaller objective value in minimization problems. For maximization problems, we reverse the sign of the Ising model parameters for the objective function and then minimize it. Constraints in optimization problems should also be composed of sets of terms in the Hamiltonian, where the lowest energy on the terms leads to the satisfaction of the corresponding constraints. Therefore, the annealing process needs to obtain feasible solutions that satisfy all constraints with a good quality of the objective value by exploring the lowest energy state of the total Hamiltonian.
Note that we can easily replace the type of decision variables from {+1, −1} to {0, 1}, and vice versa.

Petri Net Fundamentals
A Petri net PN = (P, T, F, M 0 ) is a directed bipartite graph N = (P, T, F ) with initial marking M 0 . The bipartite graph has two types of vertex sets: a set of places P = {P 1 , P 2 , ..., P |P | } and a set of transitions T = {T 1 , T 2 , ..., T |T | }. The arc function F : (P × T ) ∪ (T × P ) → N defines the connection from places to transitions, and vice versa. The returned values of F indicate the number of removed tokens from the starting place when (P × T ) → N or the number of generated tokens to the ending place when (T × P ) → N upon the firing of the transition.
Instead of an arc function, a matrix representation is often used for the connection. Here, Pre and Post show the incident matrix of size |P | × |T | to indicate the connection from places to transitions and from transitions to places, respectively. It should be noted that Pre(P i , T j ) = F ((P i , T j )) and Post(P i , T j ) = F ((T j , P i )).
Tokens are located in places such that the token distribution on P represents the status of the modeled system. The vector M k of size |P | represents the number of tokens in each place P i ∈ P at step k. In addition, M 0 , called the initial marking, is the marking at step 0, which corresponds to the initial status of the system. Therefore, the Petri net P N = (N, M 0 ) represents the structure and initial states of the modeled system.
denote the set of input and output places of T j , respectively. Similarly, • P i and P • i are the sets of input and output transitions of P i , respectively. An enabled transition T j can fire. Firing implies the occurrence of an event in the system. The firing of transition T j removes Pre(P i , T j ) tokens from each P i ∈ • T j and adds Post(P i , T j ) tokens to each P i ∈ T • j . Firing changes the token distribution, which indicates the change in status of the system by the corresponding event. The following mathematical form shows the change in status caused by the firing of transitions at step k: where X k is a vector of size |T | showing the number of times each T j fires at step k, and X k is called the firing count vector at step k.
For a quantitative analysis of the dynamical behavior of a system, time was introduced to Petri nets [van der Aalst, 1996]. The timing methods can be categorized into three types: firing duration (FD), holding duration (HD), and enabling duration (ED) [van der Aalst, 1996]. The FD assigns time to transitions to represent the firing duration. The HD is referred to as place time Petri nets, where tokens are unavailable for firing for a particular period after being located in the place. In the last method, i.e., the ED, a transition cannot be fired for a given period after being enabled. Although other types can be allowed with our method, this paper focuses on timed Petri nets with the FD.
A timed Petri net is a six-tuple T P N = (P, T, F, T S, F D, M 0 ), where T S is the set of time values, and F D : T → T S is a function that returns the firing duration time of transition T j ∈ T . In this study, we assume that T S = N and is a set of natural numbers. A timestamp is also attached to tokens to record the token generation time. In the timed Petri net, the transition T j is enabled at time k when each input place P i of T j has more than or equal to F (P i , T j ) tokens, and its time stamp is no more than k. By firing T j at time k, the marking is changed according to the same rule as the Petri net described above, except that we attach the time stamp k + F D[T j ] to each output token.
A colored Petri net is an extension of ordinary Petri nets, where we can introduce values (colors) to tokens, guard functions to transitions for their firing conditions, and arc functions to output arcs of transitions to calculate the colors of the output tokens [Jensen et al., 2007]. Therefore, tokens become much more informative and allow us flexible and efficient modeling.
This extended Petri net is denoted by CPN = (Σ, P, T, F, V, C, G, E, M 0 ), where P , T , and F represent a set of places, transitions, and arcs, respectively. Σ shows a set of colors, and C : P → Σ is a color function for the places. Only tokens with colors specified by C(P i ) are located in place P i . In addition, V is a set of arc variables, where Type(v) ⊆ Σ for v ∈ V . Moreover, Type(v) is the type of variable v. Here, E denotes an arc function where E((T j , P i )) returns tokens on C(P i ) for each output place P i based on the binding values to each of the arc variables. The term G represents a guard function that returns a Boolean value to determine whether the corresponding transition is enabled. A marking M k at step k is a mapping of multiple sets of C(P i ) to each place. Transition T j is enabled at marking M k with binding b when for each input place P i of T j , where E((P i , T j ))(b) denotes the result of the arc function for arc (P i , T j ) when we bind b(v) for each variable v in E((P i , T j )). At marking M k , M k+1 is generated by firing T j with binding b: Colored Petri nets are extremely powerful in the sense that complicated systems can be easily modeled.

Binary Quadratic Nets
This section introduces a new Petri net class called binary quadratic nets to represent the Ising and QUBO models. The Ising model is a mathematical model in statistical mechanics. The model is a graph in which the vertices correspond to spins, and the edges interact between spins. Petri nets are also suitable for expressing Ising models because their fundamental components, places, and transitions can naturally represent the states of the spins and their interactions. A token in a place can show the corresponding spin state when we introduce colors {−1, +1} to tokens. For QUBO models, the color type of tokens becomes {0, 1}.
In our Petri net modeling-based Ising and QUBO model formulation, Petri net models representing target optimization problems were converted into the corresponding binary quadratic nets. In other words, our proposed method is a net transformation from problem-domain Petri nets into binary quadratic nets. Binary quadratic nets are entirely equivalent to Ising or QUBO models.

Formal Definition
Definition 1 (Binary Quadratic Net) A binary quadratic net is a colored Petri net denoted by BQN = (Σ,P ,T ,F ,Ĉ, w):Σ There is one token in each place, where a place corresponds to a spin variable in the Ising models and a binary variable in the QUBO models, respectively. For simplicity, we denote the value of a variable associated with placep i by M (p i ), that is, M (p i ) ∈ {−1, +1} in the Ising models or M (p i ) ∈ {0, 1} in the QUBO models. Note that in the QUBO model, M (p i ) = 0 indicates that the placep i includes a token with color 0.
Current quantum annealing platforms have various physical topologies, chimera graphs for D-waves, and threedimensional lattices for CMOS annealers [OKU et al., 2019]. Thus, we need to embed the logical Ising model to a specific graph topology to run the annealing process. However, we do not treat the embedding process in this study because we can use existing embedding algorithms for each platform. Moreover, Fujitsu Digital Annealer allows fully connected topologies. For reference, Fig. 2 depicts the binary quadratic net model for the 2D array Ising model shown in Fig.1. To combine binary quadratic nets and Ising models, we introduce a measure to represent the state of binary quadratic nets with a marking M. We define this as an energy function of binary quadratic nets corresponding to the Hamiltonian in the Ising models.
Definition 2 (Energy Function of Binary Quadratic Net) For a binary quadratic net BQN = (Σ,P ,T ,F ,Ĉ,ŵ) with a marking M, the energy function is defined as follows: The first summed terms show the energy derived from token existence at each place. The second summed term represents the energy from the interaction between tokens on both sides of each transition. The weight parameters w(p i ) and w(t i,j ) are defined forp i andt i,j , respectively. The energy function in Definition 2 is equivalent to the Ising model shown in (1), and the energy function is uniquely defined from the corresponding binary quadratic net. Although it is straightforward, we simply summarize this fact as a proposition for the readability of the remaining parts.
Proposition 1 For an Ising model, we have an equivalent binary quadratic net BQN in the sense that the energy function H BQN (M) is equivalent to the Hamiltonian of the Ising model.
In our binary quadratic nets, transitions represent interactive relations between tokens on both sides. There may be numerous interaction types. We choose interaction types carefully depending on the target optimization problems. In the Appendix, we summarize the primitive interaction for QUBO and the Ising model, respectively. Table A indicate inconsistency and a tautology, respectively. Table A.2 lists the interactions for the Ising model converted from those in Table A.1 by applying the following relation.
where M (p i ) ising and M (p i ) qubo denote the marking in the Ising and QUBO models, respectively.

Binary Quadratic Net Examples
As we describe in Proposition 1, binary quadratic net models are equivalent to the Ising models. We show the binary quadratic net models for well-known graph partitioning and minimum vertex cover problems. These problems are straightforwardly modeled with binary quadratic nets and formulated as marking problems in such nets.
Example 1 Minimum Vertex Cover: For a given undirected graph G = (V, E), with vertex set V and edge set E, a vertex cover set satisfies the condition such that every edge in E is incident on a vertex in the cover set. The problem is to find the minimum vertex cover set. The problem is known as NP-hard[Garey and Johnson, 1990].
We choose the QUBO model for the problem, where we generate the place set and transition set of the binary quadratic net under the one-to-one correspondence with the vertex set and the edge set, respectively. We express a vertex cover by a marking, where we place a token with a color of 1 to show that the corresponding vertex is in the vertex cover, and where a color of 0 is not in the cover. To satisfy the feasibility of the vertex cover, we need to avoid a token with a color of 0 in both input places of each transition because the situation does not cover the corresponding edge. We formulate the feasibility of the condition corresponding to I qubo 8 in Table A.1 as follows: To minimize the objective function, we attempt to reduce the number of color 1 tokens in M. Therefore, the following penalty function is suitable because only color 1 tokens increase the penalty.
Based on the superposition principle, we have the total formulation based on our binary quadratic net models by combining the binary quadratic nets corresponding to (14) and (15): where H constraint (M) counts the number of transitions, both of which have a 0 color token. In addition, H cost (M) denotes the number of places colored with 1 to minimize the vertex cover. Moreover, A and B show the parameters for a trade-off between the constraint and the objective function. The formulation is equivalent to the QUBO model in [Lucas, 2014].
Consider the graph shown in Fig. 3 as a problem instance of the vertex cover problem. We can formulate the problem instance as a marking problem in the corresponding binary quadratic net in Fig. 4. Marking M = (1, 0, 1, 0, 0, 1, 0, 1) is a feasible solution to the instance.
Example 2 Graph Partitioning: Graph partitioning is an optimization problem, which is explained as follows: Let us consider an undirected graph G = (V, E) with vertex set V and edge set E, where |V | is the number of vertices. The problem is partitioning V into two subsets whose sizes are equal to |V |/2, minimizing the number of edges between the subsets. The problem is known as NP-hard[Garey and Johnson, 1990].
We reduce the problem to a marking problem in binary quadratic nets, where each vertex in the given graph corresponds to a place in the binary quadratic net, and the marking such that M (p i ) ∈ {−1, +1}, ∀p i ∈P shows the partition. The objective function minimizes the number of edges connecting the two partitioned groups. The transitions have a one-to-one correspondence with the edges. To minimize the objective function, we want to reduce the different color tokens in pairs of places that share the output transition. We design the objective function with I Ising 6 in Table A.2 as follows: 1,2 1,2 t tt t 2,6 2,6 t tt t 1,6 1,6 t tt t 6,7 6,7 t tt t Figure 4: Binary Quadratic Net for Vertex Cover Problem Instance shown in Fig.3 Here, I Ising 6 outputs 1 only when both places have different color tokens, that is, XOR.
Second, we consider the constraint of the graph-partitioning problem and the equality in the vertex size of both partitioned groups. Concerning {−1, +1} logic, we can design the following energy function for the constraint because minimizing the function leads to a satisfaction of the constraint: The following penalty function shows the total energy function for graph partitioning and is equivalent to the Ising model presented in [Lucas, 2014].

Binary Quadratic Net Construction from Problem Domain Petri Nets
In the previous section, we modeled combinatorial optimization problems directly, a minimum vertex cover problem, and a graph partitioning problem with binary quadratic nets and formulated them as marking problems. However, they are exceptional cases. In general, we can meet the conceptual gap between combinatorial optimization problems and target binary quadratic nets. Note that the problem is more serious in the direct formulation of the Ising or QUBO models.
Our approach attempts to minimize the gap by converting problem-domain Petri net models, represented by timed Petri nets and colored Petri nets, into the corresponding binary quadratic nets. This section proposes a method for constructing target binary quadratic nets from problem-domain Petri net models.

Incremental Construction based on Superposition Principle
A binary quadratic net can be composed incrementally by combining binary quadratic subnets, each corresponding to a constraint or an objective function in the original optimization problem. This superposition principle of the net structure and weight values simplifies the binary quadratic net construction. The composition is straightforward, where we add the weight values on the places and transitions if the same places and transitions in the subnets are to be combined; otherwise, add new places and transitions with their weight values. In the previous section, we observed this process in the graph partitioning example (Example 2). Even though we did not consider the weight values, we combined the two subnets. The following definition formally represents the binary quadratic net composition.
Definition 3 (Binary Quadratic Net Composition) For two given binary quadratic nets BQN h = (Σ,P h ,T h ,F h ,Ĉ h , w h ) and BQN k = (Σ,P k ,T k ,F k ,Ĉ k ,ŵ k ) such that the model types of both nets, Ising or QUBO, are the same, the new binary quadratic net BQN = (Σ,P ,T ,F ,Ĉ,ŵ) is composed based on the following superposition principle: The following proposition is a straightforward but essential property for validating our method.
Proposition 2 (Superposition Property) Let BQN be composed from BQN h and BQN k . The following properties hold.
where A is a scalar, and H A·BQN (M) is the energy function of BQN, such that we replace the weight functionŵ with A ·ŵ.
Proof 1 In Definition 3, the energy function is defined by the net structure and weight functionŵ. Based on the composition rules in (21), (22), and (23), all structural properties of BQN h and BQN k are transformed into BQN. Based on the rule for the weight function (25), we can confirm thatŵ(p i ) andŵ(t i,j ) can be divided intoŵ h (p i ) + w k (p i ) andŵ h (t i,j ) +ŵ k (t i,j ) for the common placep i and transitiont i,j , respectively. Therefore, the additivity in (26) holds.
The homogeneity property with a degree of 1 is given by Definition 2 if we replace the weight functionŵ with A ·ŵ.
Owing to the additivity and homogeneity properties in Proposition 2, we have the following corollary: Corollary 1 Binary quadratic nets can be composed by the incremental application of Definition 3, in which we can scale the weight functionŵ in subnets with a constant factor.
In our approach, we construct a target binary quadratic net based on the incremental compositions of binary quadratic subnets. Each subnet is converted from a property of the Petri net model representing the optimization problems. We call the Petri net problem-domain Petri net. The properties of the problem-domain Petri net models are expressed with marking or firing sequences. To focus on markings or firing sequences, we employ marking-based or firing-based constructions. In the following subsections, we assume that binary quadratic nets are QUBO models unless otherwise stated, but can be converted into Ising models by using the conversion rule (58).

Marking-based Construction
Let us denote a problem-domain Petri net for a target optimization problem by N = (P, T, F ) with P = {P 1 , P 2 , ..., P n } and T = {T 1 , T 2 , ..., T m }. Marking M k (P j ) represents a set of tokens in place P j at time step k. In addition, M k is a vector of size |P |, showing a token distribution on the Petri net at time step k. The marking trajectory of a Petri net model, M 0 , M 1 , ..., M K denotes the status changes triggered by the firing of transitions.
In the marking-based construction of binary quadratic nets, we represent a marking trajectory of the problem domain Petri net by the place set of the target binary quadratic net. If each place P i in the problem-domain Petri net has at most one token at any time and the maximum step is K, we initially prepare the following place set of the target binary quadratic net.P = {p k i |P i ∈ P, k = 0, 1, 2, ..., K} If M k (P i ) is a natural number, that is, more than one token or a color token with a natural number value is possible in each place in P , we need to prepare more places because each place in the binary quadratic nets has one token with color in {0, 1} or {−1, +1}.P = {p k i,n |P i ∈ P, n = 0, 1, 2, ..., N, k = 0, 1, 2, ..., K}, (29) where N is the possible maximum value of M k (P i ). At the same time, we require a one-hot constraint because only one place betweenp k i,n for n = 0, 1, ..., N for each i and k should be marked with 1, and the others should be marked with 0: As the energy function representation, we need the following: Note that we can obtain the corresponding binary quadratic net from an energy function.

Boundedness
The boundedness is an essential characteristic to ensure avoiding an overflow in the system behavior. In the Petri net theory, we can express this characteristic as follows.
where U i is the upper bound for P i . We can convert the constraint into a binary quadratic net and represent it as the energy function under the one-hot constraint (31): Function (33) is sufficient for the equality constraint M k (P i ) = U i , but not for the upper bound. Therefore, we improve the function by introducing ancilla places,û i,m , i = 1, 2, ..., , |P |, m = 0, 1, ..., U i . This technique is well known in the Ising model formulation [Lucas, 2014], [Tanahashi et al., 2019].
The upper bound constraint appears in numerous optimization problems. The knapsack constraint is a well-known example of this constraint for the knapsack place. We can also express the boundedness of specific places P i by removing the other places from the function.

Invariant
The boundedness shown in (32) denotes inequality constraints. The invariance leads to equality constraints based on markings. Note that the invariance is different from the structural invariance of the net theory. The behavioral invariance requires that the total weighted sum of the tokens be equal among the firing sequences. The sum may be the cost required for the resource to operate the system. The following constraint shows that the weighted sum becomes W for all k.
We can convert the constraint into a binary quadratic net and represent it as the energy function under the one-hot constraint (31):

Firing-based Construction
A firing count vector X k represents the firing counts of each transition in a Petri net model at step k, where X k (T j ) denotes the firing counts of transition T j at time step k. Elements of X k are usually one of {0, 1}; that is, transitions can fire only once at each time step. We call this single firing restriction single-server semantics in Petri net theory. However, any natural number values (allowing more than one) are also possible with infinite server semantics. For the single-server semantics, we generate a place setP in the target binary quadratic net, where M (p k i ) corresponds to Similarly, we can extend the single-firing semantics to the N times firing semantics. Note that an infinite number of firing times is possible mathematically but impossible practically. Thus, we restrict the maximum number reasonably to N .P = {p k i,n |T i ∈ P, n = 0, 1, 2, ..., N, k = 0, 1, 2, ..., K}, As the energy function representation, we need the same one-hot constraint (31):

Resource Conflict
If transitions T i and T j have a common input place and conflict with a single token, X k (T i ) = X k (T j ) = 1 cannot be allowed. Therefore, we need the following energy function: where C is a set of conflict transitions. In addition, C can be extracted from the problem domain Petri net.
In timed Petri net models, we need to consider the firing duration. Let N = (P, T, F, T S, F D) with P = {P 1 , P 2 , ..., P n } and T = {T 1 , T 2 , ..., T m } be a timed Petri net, where F D : T → N is a function that returns the firing duration. We then extend the resource conflict in the stepwise firing into the timed firing.
where C timed is the timed conflict set such that conflict transitions cannot fire until the firing duration is complete. We can obtain C timed from the given timed Petri net.

Firing Count
In some applications, we must specify the number of firing occurrences for each transition.
where F C i is the specified number of firings of T i .
We can convert the constraint into a binary quadratic net and represent it as the energy function under the one-hot constraint (31): Note that M (p k i,n ) in the binary quadratic net corresponds to X k (T i ) in the problem domain Petri net. Assuming that each transition T i should fire exactly once during X 0 , X 1 , ..., X K , the function (42) can be represented as follows: In practical cases, this constraint is commonly used.

Precedence Relation
Let us assume again that each transition T i should fire exactly once during X 0 , X 1 , ..., X K . We consider the precedence relation between the firing of transitions. If T i precedes structurally to T j , that is, T • i ⊆ • T j and • T i ⊆ T • j , the following penalty function is required.
where P rec is the set of precedence relations. Here, P rec can be extracted from the problem domain Petri net. Similar to the firing conflict, we consider timed Petri nets by introducing the firing duration.

Application Example 1 (Marking-based Construction): Traveling Salesman Problems
We consider the traveling salesman problem as an example of a marking-based construction. Let N = (P, T, F, T S, F D) with P = {P 1 , P 2 , ..., P n } and T = {T i,j |∀(P i , P j ) ∈ P × P } be a timed Petri net, where F D : T → N is a function that returns the firing duration. Each place P i corresponds to a city, and the transition T i.j denotes the path from P i to P j . The firing of T i,j indicates the movement from P i to P j . The initial marking M 0 includes only one token at P 1 .
Because the salesman visits each city only once from the problem definition, In addition, because the salesman should be at one place for each step, The objective function is to minimize the total time to visit all cities and return to the starting place.
By converting the constraints and objective function in the problem-domain Petri net into the binary quadratic nets, we have The total binary quadratic net is as follows: where A, B, and C denote the scale factors of the corresponding subnets. The formulation is equivalent to the QUBO model in [Lucas, 2014].
The following are the energy functions of the binary quadratic subnets corresponding to the precedence relation constraint between tasks in each job, the resource conflict constraint for each shared resource, and the firing count constraint for each task.
The precedence set P rec in (54) is easily constructed by extracting the precedence relation in each sequential system (job). We can extract the timed conflict set C timed in (55) from the connection between the transitions and resource places and the firing duration F D. MaxTime is the only parameter we need to set before the optimization process. This value indicates the delivery time deadline. The total binary quadratic net is as follows: Note that this formulation is used to obtain a feasible solution; however, we can extend this to the optimization by incorporating it with a binary search to find the minimum MaxTime. This formulation is equivalent to the QUBO model in [Venturelli et al., 2016]. Figure 7 shows the schedule obtained through the annealing process using a Fujitsu Digital Annealer.

Conclusions
This paper proposes a Petri net modeling approach to the Ising model formulation for quantum annealing. Although our method requires users to model their optimization problems with Petri nets, this process can be carried out in a relatively straightforward manner if we know the target problem and the simple Petri net modeling rules. Therefore, we can drastically relax the difficulty of the Ising model formulation. We implemented our method with Python incorporated using well-known Petri net tools, CPNTools and SNAKES. We can automatically generate the Ising models for optimization problems, such as scheduling problems, vehicle routing problems, portfolio optimization problems, and others once we model the target optimization problems with Petri nets.
We defined binary quadratic nets to represent the Ising model formulation. However, the binary quadratic net can also be used to analyze the quantum annealing process by attaching an additional subnet that simulates the annealing. This tool may contribute to parameter tuning for annealing, which is another task used to expand the quantum annealing technology mentioned in the Introduction.  are the inconsistency and tautology, respectively. The others also show important logical functions. Table A.2 shows the interaction primitives for the Ising model converted from Table A.1 by using the following relation.