- Email: [email protected]
A Novel Method for Reachability Determination in Petri Nets
Proceedings of the 13th IFAC Symposium on Information Control Problems in Manufacturing Moscow, Russia, June 3-5, 2009
A Novel Method for Reachability Determination in Petri Nets ⋆ Duan Li ∗ Xiaoling Sun ∗∗ Jianjun Gao ∗∗∗ Shenshen Gu ∗∗∗∗ Xiaojin Zheng † ∗
Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) ∗∗ Department of Management Science, Fudan University, Shanghai, China (e-mail: [email protected]) ∗∗∗ Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) ∗∗∗∗ Department of Mechanical and Automation Engineering, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) † Department of Mathematics, Shanghai University, Shanghai, China (e-mail: [email protected]) Abstract: Reachability is one of the most important behavioral properties of Petri nets and the past four decades have witnessed great efforts in developing various implementable methodologies in determining reachability of Petri nets. We propose in this paper a novel method for solving the fundamental equation in the reachability analysis, which has been known to be NP-complete. More specifically, by adopting a revised version of the cell enumeration method for an arrangement of hyperplanes in discrete geometry, we develop an efficient solution scheme to identify firing count vector solution(s) to the fundamental equation on a bounded integer set, with a complexity bound of O((nω)n−m ), where n is the number of transitions, m is the number of places and ω is the upper bound of the number of firings for every transition. Keywords: Petri nets, reachability analysis, cell enumeration, linear Diophantine equations on bounded integer set, integer programming, discrete optimization. 1. INTRODUCTION Petri net, introduced by Petri in his seminal work “Communication with Automata” (Petri (1962)), has been a promising mathematical formalism for modeling, analyzing and designing discrete event systems, especially by its remarkable capability in modeling process synchronization, asynchronous events, concurrent and distributed operations and conflicts or resource sharing. The past four decades have witnessed innumerable successful applications of Petri nets in various areas, such as i) modeling and analysis of communication protocols and networks, production systems, sequence controllers and software development, and ii) performance evaluation of the modeled systems. The Petri net is a particular kind of bipartite directed graph consisting of three types of elements: places, transitions, and directed arcs connecting places and transitions, while each place may hold a non-negative number of tokens. See Fig. 1 1 for an illustrative example of Petri nets with 4 places and 6 transitions. The distribution of tokens ⋆ This research is partially supported by University Grants Council, Hong Kong, under CUHK414207, and Shanghai Leading Academic Discipline Project (B210).
978-3-902661-43-2/09/$20.00 © 2009 IFAC
Fig. 1. Example of a Petri net on places, called Petri net marking, defines the state of the modeled system. A marking for a Petri net with m places is represented by an (m × 1) vector M, where Mj , j = 1, 2, . . . , m, are nonnegative integers representing the number of tokens in the corresponding places. A general Petri net is characterized by a quintuple (P, T, I, O, M0 ) (see Zurawski and Zhou (1994)), where P = {p1 , p2 , . . . , pm } is a finite set of m places, T = {t1 , t2 , . . . , tn } is a finite set of n transitions, DI : (P × m×n T ) → Z+ is an input function that defines weights
276
10.3182/20090603-3-RU-2001.0100
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
associated with directed arcs from places to transitions, m×n m ) is the set of (m × n) dimensional mawhere Z+ (Z+ trices (m dimensional vector) with all entries being in Z+ , m×n DO : (P × T ) → Z+ is an output function that defines weights associated with directed arcs from transitions to m places, and M0 : P → Z+ is the initial marking. For the Petri net given in Figure 1, we have 0 1 0 0 1 0 0 0 1 0 0 1 DI = , 0 0 0 4 0 0 2 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 , DO = 0 0 1 0 1 0 0 0 0 1 0 5 T and M0 = (1, 1, 0, 0) . The change of the distribution of the tokens represents the dynamics of the modeled system, while the distribution of tokens on places may change according to the following enabling rule and firing rule (see Zurawski and Zhou (1994)): Enabling Rule: A transition t is said to be enabled if each input place p of t contains at least the number of tokens equal to the weight of the arc connecting p to t. Firing Rule: (a) An enabled transition t may or may not fire depending on the additional interpretation, and (b) A firing of an enabled transition t removes form each input place p the number of tokens equal to the weight of the directed arc connecting p to t, and deposits in each output place p the number of tokens equal to the weight of the arc connecting t to p. Reachability is one of the most important behavioral and structural properties of Petri nets. Given both an initial state M0 and a target terminal state M , a natural question to ask is whether or not we have a sequence of firing rules such that the system can reach the specific target state within finite steps. There are two primary approaches in investigating reachability: Reachability graph analysis and reachability algebraic analysis. Reachability graph analysis is based on the creation and investigation of a reachability graph or a reduced reachability graph. As this approach relies on exhaustively generating all the reachable markings from a given initial marking, it suffers the state explosion phenomenon, while transferring the reachability graph to a reduced counterpart is NP-hard (see Kostin (2003)). The second approach is based on methods of linear algebra. It is well known (see Murata (1977) and Yen (2006)) that a necessary condition for reachability of marking M from an initial marking M0 of a Petri net is that there exists a nonnegative integer vector solution of the following system of linear equations, Ax = b,
(1) T
where b = M − M0 , x = (x1 , x2 , . . . , xm ) is a firing count vector and A is an (m × n) dimensional incidence matrix given by A = [a(pi , tj )] = [DO (pi , tj ) − DI (pi , tj )]. If we set in Figure 1 the final marking to be M = (2, 0, 1, 3)T , the fundamental equation for this example is then given by
2 0 0 −2
−1 1 0 0
0 −1 1 0
0 0 −4 1
−1 0 1 0
0 1 −1 −1 x= . 0 1 5 3
When the Petri net is acyclic, the above condition becomes both necessary and sufficient (Matsumoto et al. (1997)). Equation (1), termed the fundamental equation in Petri nets, reveals that the investigation of reachability can be achieved, to certain degree, by finding out nonnegative integer solution(s) to a system of linear equations. The first step in the reachability algebraic analysis is to find firing count vectors, which can be performed either by solving the fundamental equation directly, or by solving the corresponding integer programming problem. Unfortunately, finding nonnegative integer solutions to the fundamental equation in Petri nets is NP-complete (Matsumoto et al. (2002)). The second step in the reachability algebraic analysis is to translate a firing count vector into a firing sequence, if there is a reachable path, while such an algorithm (see, e.g., Kostin (2003)) is of a combinatorial nature. In essence, solving the fundamental equation is an integer programming problem (see, for example, Schrijver (1986) Nemhauser and Wolsey (1988) and Li and Sun (2006)). Generally speaking, the capability of the existing methods for finding the nonnegative integer solution to the fundamental equation (1) is still very limited. In particular, the Petri nets community has been largely dependent on some existing softwares to obtain nonnegative integer solutions to the fundamental equation. More efficient solution schemes that invoke the state-of-the-art from some other related subjects, especially from discrete optimization and integer programming, deserve further investigation. There have been a few attempts in the Petri net communities in this direction. For example, Matsumoto et al. (2002) applied Groebner Bases method from integer programming to solve reachability problems in Petri nets. The solution approach using Groebner Bases depends heavily on symbolic computing power, and its capability in solving large-scale problem is thus not promising. Although the original analytical reachability analysis is to find nonnegative integer solution(s) to the fundamental equation, a more practical problem is to find integer solution(s) to the fundamental equation on a bounded integer set, X = {x ∈ Z n | 0 ≤ xi ≤ ui , i = 1, . . . , n}, where Z n denotes the set of all integer points in Rn . This confinement is reasonable as there always exist resource constraints in real-world applications of Petri nets. For example, in chemical processes, the operation of reactors (transitions in a corresponding Petri net) is expensive, thus the number of such operations should not exceed a given upper bound. It is well known that linear Diophantine equations, Ax = b, x ∈ Z n , are polynomially solvable (see, e.g., Schrijver (1986) and Nemhauser and Wolsey (1988)), while linear Diophantine equations on a bounded integer set, Ax = b, x ∈ X, where A is an m × n integral matrix with rank(A) = m, and b ∈ Rm being integral, are NP-complete, as the special case of linear Diophantine equations with m = 1, x ∈ {−1, 1}n and b = 0 is NPcomplete (Garey and Johnson (1979)).
277
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
We propose in this paper a novel method for solving the fundamental equation in the reachability analysis over a bounded integer set X of firing count vectors. Our method is based on a recognition that whether or not there exists a solution to the fundamental solution is equivalent to whether or not the distance from X to the affine solution set of Ax = b in Rn attains zero. Furthermore, from recent results in Sun et al. (2008), finding the distance from an integer set to an affine set can be efficiently achieved by the cell enumeration method for an arrangement of hyperplanes in discrete geometry. This cell enumeration approach provides a promising platform for designing an efficient method with a complexity of O((nω)n−m ), where ω = maxi=1,...,n ui , to find integer solutions to the fundamental equation on a bounded integer set . 2. SOLVING THE FUNDAMENTAL EQUATION BY CELL ENUMERATION Let U = (U1 , . . . , Un−m ) be an orthogonal basis for the null space of A and x0 be a special solution to Ax = b in Rn , then C = {x ∈ Rn | Ax = b} can be expressed by n
C = {x ∈ R | x = x0 +
n−m
(n − m)-dimensional polyhedral set formed by the halfspaces induced by hij ’s and can be characterized by a n 1 n u j=1 j -dimensional sign vector, sign(E) = (w , . . . , w ), where wj = (w1,j , . . . , wuj ,j ) is specified by
+, if gij (π) > 0 wij = (5) , i = 1, . . . , uj , −, if gij (π) < 0 with π being an interior point of E. For fixed j, since hij (i = 1, . . . , uj ) are parallel hyperplanes, wj = (w1,j , . . . , wuj ,j ) must take the following form: i
(+, . . . , +, −, . . . , −),
i = 0, 1, . . . , uj .
The one-to-one map between all distinct ϕ(y) for y ∈ C and the sign vectors of all the cells of the hyperplane arrangement can be established as follows: w = (w1 , . . . , wn ) ←→ ϕ = (ϕ1 , . . . , ϕn ), where for j = 1, . . . , n, wj = (−, . . . , −) ←→ ϕj = 0, wj = (+, . . . , +) ←→ ϕj = uj ,
zi Ui ,
i
w = (+, . . . , +, −, . . . , −) ←→ ϕj = i, j
i=1
zi ∈ R, i = 1, . . . , n − m}.
(2)
Consider the Euclidian distance from C to X = {x ∈ Z n | 0 ≤ xi ≤ ui , i = 1, . . . , n}: δ = dist(C, X) = min{x − y | x ∈ X, y ∈ C}.
(3)
It is obvious that Ax = b has a solution in X if and only if δ = 0. Furthermore, when δ = 0, any x∗ ∈ X that achieves the zero distance is an integer solution to the linear equations Ax = b in X. Note that combining (2) and (3) yields
where i ∈ {1, . . . , uj − 1}.
It has been known that the number of cells of the hyperplane arrangement generated by (4) is bounded by O((nω)n−m ), where ω = maxi=1,...,n ui (see Zaslavsky (1975) and Edelsbrunner (1987)). Moreover, using the cell enumeration methods in Avis and Fukuda (1996) and Sleumer (1999), we can find all the cells of the hyperplane arrangement generated by (4) in O((nω)n−m ) time. Listing all the distinct ϕ(y) for y ∈ C as ϕ1 , . . . , ϕk , the distance δ = dist(C, X) can then be calculated as follows: δ = min δj ≡ min dist(C, ϕj )
δ = dist(C, X)
j=1,...,k
= min min{y − x | 0 ≤ xi ≤ ui , xi ∈ Z, i = 1, . . . , n}
= min (U U T − In )(ϕj − x0 ).
y∈C
j=1,...,k
= min y − ϕ(y), y∈C
where, for j = 1, . . . , n, ϕ(y)j is determined by the following equation: 1 0, yj ∈ (−∞, ], 2 1 1 ϕ(y)j = i, yj ∈ (i − , i + ), i ∈ {1, . . . , uj − 1}, 2 2 uj , yj ∈ [uj − 1 , ∞). 2
Let’s consider the hyperplane arrangement generated by n the following j=1 uj hyperplanes in Rn−m :
+
l=1
1 zl Ulj − (i − ) = 0} 2
(6)
Therefore, finding a solution of the linear system Ax = b over the integer set X or checking its infeasibility can be done in O((nω)n−m ) time. Especially, when r = n − m is fixed (0 ≤ r ≤ n − 1), the linear system Ax = b over X is polynomially solvable. Let us consider now the example in Fig. 1 to illustrate the solution algorithm for the fundamental equation using the cell enumeration method. For this instance, we assume that the bounded integer set is given by X = {x ∈ Z 6 | 0 ≤ xi ≤ 2, i = 1, . . . , 6}. Since the incidence matrix is of rank 4 in this example, solving the fundamental equation in R6 yields
hij = {z ∈ Rn−m | gij (z) := (x0 )j n−m
j=1,...,k
C = {x ∈ R6 | Ax = b} (4)
for i = 1, . . . , uj , j = 1, . . . , n. Note that for each j, hij (i = 1, . . . , uj ) are parallel hyperplanes. A cell E of the hyperplane arrangement corresponding to hij ’s is an
= {x ∈ R6 | x = x0 + z1 U1 + z2 U2 , zi ∈ R, i = 1, 2}, where x0 = (0.2912, −0.2138, 0.0100, −0.2984, −0.2038, 0.776)T and
278
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
δ5 = dist(X, Cell 5) = (U U T − I6 )(ϕ5 − x0 ) = 0.
15
g14 g26
For Cell 61 labelled in Figure 2, all points are above the 12 hyperplanes. Thus, the sign vector for Cell 61 is ((++)(++)(++)(++)(++)(++)) and the closest integer point in {0, 1, 2}6 to Cell 61 is ϕ61 = (2, 2, 2, 2, 2, 2). Furthermore, the distance between X and Cell 61 is given by
g24 Cell 61
10
g21 g15
g25
z2
g11 5
Cell 9
g16
δ61 = dist(X, Cell 61) = (U U T − I6 )(ϕ61 − x0 ) = 1.578.
g22 g12 0
Cell 5 g23 g13 −5 −10
−8
−6
−4
−2
0
z1
2
4
6
8
10
Fig. 2. Cell arrangement for the example problem. 0.3508 0.3944 −0.1605 0.6882 −0.2713 0.5636 U = (U1 , U2 ) = . 0.1477 0.1661 0.8622 0.1006 0.1108 0.1245 Note that uj = 2, j = 1, . . ., 6. Setting each xj , j = 1, . . ., 6, equal to 12 and 32 , respectively, generates the following 12 hyperplanes on C defined in (4),
g11 (z) = 0.3508z1 + 0.3944z2 = 0.2088, j=1: g21 (z) = 0.3508z1 + 0.3944z2 = 1.2088,
g12 (z) = −0.1605z1 + 0.6882z2 = 0.7138, j=2: g22 (z) = −0.1605z1 + 0.6882z2 = 1.7138,
g13 (z) = −0.2713z1 + 0.5636z2 = 0.4900, j=3: g23 (z) = −0.2713z1 + 0.5636z2 = 1.4900,
g14 (z) = 0.1477z1 + 0.1661z2 = 0.7984, j=4: g24 (z) = 0.1477z1 + 0.1661z2 = 1.7984,
g15 (z) = 0.8622z1 + 0.1006z2 = 0.7038, j=5: g25 (z) = 0.8622z1 + 0.1006z2 = 1.7038,
g16 (z) = 0.1108z1 + 0.1245z2 = −0.2762, j=6: g26 (z) = 0.1108z1 + 0.1245z2 = 0.7238. As illustrated in Fig. 2, these 12 hyperplanes partition set C into various cells, each of which corresponds to a unique sign vector defined in (5), sign(E) = (w1 , . . . , w6 ), where wj = (w1j , w2j ). Note that the arrow attached to each hyperplane points to the positive half space generated by this hyperplane. Table 1 lists all the 61 cells identified with their sign vectors and the corresponding distances to X = {0, 1, 2}6 . Let us consider Cells 5 and 61 labeled in Fig. 2. As any interior point π of Cell 5 satisfies g11 (π) > 0, g21 (π) < 0, g12 (π) < 0, g22 (π) < 0, g13 (π) < 0, g23 (π) < 0, g14 (π) < 0, g24 (π) < 0, g15 (π) > 0, g25 (π) < 0, g16 (π) > 0, and g26 (π) < 0, the sign vector for any interior point π of Cell 5 is thus specified by ((+−)(−−)(−−)(−−)(+−)(+−)). Furthermore, for all points in Cell 5, the closest integer point in X = {0, 1, 2}6 is clearly given by ϕ5 = (1, 0, 0, 0, 1, 1) and the distance between X and Cell 5 is given by
We observe from Table 1 that there exist two solutions, (1, 0, 0, 0, 1, 1)T and (1, 1, 1, 0, 0, 1)T , to the fundamental equation of this example, as they both achieve the zero distance. We need to emphasize that our solution scheme enables us to identify all solutions to the fundamental equation on a finite integer set. Note also that the total number of cells in this example is 61, which is much smaller than the upper bound of the cell numbers O((6 × 2)(6−4) ). This conclusion applies for general situations as observed from randomly generated numerical experiments presented in the next section. The flow chart given in Fig. 3 demonstrates the structure of an algorithm in converting a (feasible) firing count vector into a firing sequence, which is similar to the one in Kostin (2003). For the current instance, firing count vector (1, 0, 0, 0, 1, 1)T leads to a firing sequence of {t5 → t6 → t1 } and firing count vector (1, 1, 1, 0, 0, 1)T leads to a firing sequence {t2 → t3 → t6 → t1 }. As both firing count vectors are realizable, a decision maker would prefer firing count vector (1, 0, 0, 0, 1, 1)T to firing count vector (1, 1, 1, 0, 0, 1)T as the former one requires less number of firings, and thus smaller expenses (when assuming the same cost for firings at different transitions).
3. PRELIMINARY NUMERICAL RESULTS We present some preliminary numerical results for the solution method using cell enumeration developed in Section 2. The solution algorithm was coded by C and ran on a Pentium PC (2G Hz and 1G RAM). Our implementation of the cell enumeration method is based on the reverse search procedure proposed in Avis and Fukuda (1996) and Sleumer (1999). In our implementation, the LP solver in CPLEX 9.0 is called to solve the linear programming subproblems in the reverse search procedure. In all the test problems, we set ω = ui , i = 1, . . . , n. For each triple of (m, n, ω), 10 test problems are randomly generated. The numerical results are reported in Table 2, where |X| denotes the total number of integer points in X, the number of cells is the average for 10 test problems (after rounding off) and the CPU time is the average for 10 test problems. We point out that the purpose of implementing our proposed solution method in numerical experiments at this stage is to check the feasibility of the proposed solution method. Sophisticated implementation of the cell enumeration procedure that integrates the reverse searching technique with ZRAM library (see Ferrez et al. (2005)) will further significantly improve the performance of the method.
279
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
Table 1. A full list of cells, their sign vectors and the closest integer points. Cell j (w1 , w2 , w3 , w4 , w5 , w6 ) ϕj 1 ((+−)(−−)(−−)(−−)(−−)(+−)) (1, 0, 0, 0, 0, 1)T 2 ((−−)(−−)(−−)(−−)(−−)(+−)) (0, 0, 0, 0, 0, 1)T 3 ((+−)(+−)(−−)(−−)(−−)(+−)) (1, 1, 0, 0, 0, 1)T 4 ((+−)(−−)(+−)(−−)(−−)(+−)) (1, 0, 1, 0, 0, 1)T 5∗ ((+−)(−−)(−−)(−−)(+−)(+−)) (1, 0, 0, 0, 1, 1)T 6 ((−−)(−−)(+−)(−−)(−−)(+−)) (0, 0, 1, 0, 0, 1)T 7 ((−−)(−−)(−−)(−−)(+−)(+−)) (0, 0, 0, 0, 1, 1)T 8 ((−−)(−−)(−−)(−−)(−−)(−−)) (0, 0, 0, 0, 0, 0)T ((+−)(+−)(+−)(−−)(−−)(+−)) (1, 1, 1, 0, 0, 1)T 9∗ 10 ((+−)(+−)(−−)(−−)(+−)(+−)) (1, 1, 0, 0, 1, 1)T 11 ((+−)(−−)(−−)(−−)(++)(+−)) (1, 0, 0, 0, 2, 1)T 12 ((−−)(+−)(+−)(−−)(−−)(+−)) (0, 1, 1, 0, 0, 1)T 13 ((−−)(−−)(+−)(−−)(−−)(−−)) (0, 0, 1, 0, 0, 0)T 14 ((−−)(−−)(−−)(−−)(++)(+−)) (0, 0, 0, 0, 2, 1)T 15 ((−−)(−−)(−−)(−−)(+−)(−−)) (0, 0, 0, 0, 1, 0)T 16 ((+−)(++)(+−)(−−)(−−)(+−)) (1, 2, 1, 0, 0, 1)T 17 ((+−)(+−)(++)(−−)(−−)(+−)) (1, 1, 2, 0, 0, 1)T 18 ((+−)(+−)(+−)(−−)(+−)(+−)) (1, 1, 1, 0, 1, 1)T 19 ((++)(+−)(−−)(−−)(+−)(+−)) (2, 1, 0, 0, 1, 1)T 20 ((+−)(+−)(−−)(−−)(++)(+−)) (1, 1, 0, 0, 2, 1)T 21 ((++)(−−)(−−)(−−)(++)(+−)) (2, 0, 0, 0, 2, 1)T 22 ((−−)(+−)(++)(−−)(−−)(+−)) (0, 1, 2, 0, 0, 1)T 23 ((−−)(+−)(+−)(−−)(−−)(−−)) (0, 1, 1, 0, 0, 0)T 24 ((−−)(−−)(++)(−−)(−−)(−−)) (0, 0, 2, 0, 0, 0)T 25 ((−−)(−−)(−−)(−−)(++)(−−)) (0, 0, 0, 0, 2, 0)T 26 ((++)(++)(+−)(−−)(−−)(+−)) (2, 2, 1, 0, 0, 1)T 27 ((+−)(++)(++)(−−)(−−)(+−)) (1, 2, 2, 0, 0, 1)T 28 ((+−)(++)(+−)(−−)(+−)(+−)) (1, 2, 1, 0, 1, 1)T 29 ((++)(+−)(+−)(−−)(+−)(+−)) (2, 1, 1, 0, 1, 1)T 30 ((++)(+−)(−−)(−−)(++)(+−)) (2, 1, 0, 0, 2, 1)T 31 ((++)(−−)(−−)(+−)(++)(+−)) (2, 0, 0, 1, 2, 1)T 32 ((−−)(++)(++)(−−)(−−)(+−)) (0, 2, 2, 0, 0, 1)T 33 ((−−)(+−)(++)(−−)(−−)(−−)) (0, 1, 2, 0, 0, 0)T 34 ((++)(++)(++)(−−)(−−)(+−)) (2, 2, 2, 0, 0, 1)T 35 ((++)(++)(+−)(−−)(+−)(+−)) (2, 2, 1, 0, 1, 1)T 36 ((++)(+−)(+−)(−−)(++)(+−)) (2, 1, 1, 0, 2, 1)T 37 ((++)(+−)(−−)(+−)(++)(+−)) (2, 1, 0, 1, 2, 1)T 38 ((++)(−−)(−−)(+−)(++)(++)) (2, 0, 0, 1, 2, 2)T 39 ((−−)(++)(++)(−−)(−−)(−−)) (0, 2, 2, 0, 0, 0)T 40 ((++)(++)(++)(+−)(−−)(+−)) (2, 2, 2, 1, 0, 1)T 41 ((++)(++)(++)(−−)(+−)(+−)) (2, 2, 2, 0, 1, 1)T 42 ((++)(++)(+−)(+−)(+−)(+−)) (2, 2, 1, 1, 1, 1)T 43 ((++)(++)(+−)(−−)(++)(+−)) (2, 2, 1, 0, 2, 1)T 44 ((++)(+−)(+−)(+−)(++)(+−)) (2, 1, 1, 1, 2, 1)T 45 ((++)(+−)(−−)(+−)(++)(++)) (2, 1, 0, 1, 2, 2)T 46 ((++)(−−)(−−)(++)(++)(++)) (2, 0, 0, 2, 2, 2)T 47 ((++)(++)(++)(+−)(+−)(+−)) (2, 2, 2, 1, 1, 1)T 48 ((++)(++)(++)(+−)(−−)(++)) (2, 2, 2, 1, 0, 2)T 49 ((++)(++)(+−)(+−)(++)(+−)) (2, 2, 1, 1, 2, 1)T 50 ((++)(+−)(+−)(+−)(++)(++)) (2, 1, 1, 1, 2, 2)T 51 ((++)(++)(−−)(+−)(++)(++)) (2, 2, 0, 1, 2, 2)T 52 ((++)(+−)(−−)(++)(++)(++)) (2, 1, 0, 2, 2, 2)T 53 ((++)(++)(++)(+−)(++)(+−)) (2, 2, 2, 1, 2, 1)T 54 ((++)(++)(++)(+−)(+−)(++)) (2, 2, 2, 1, 1, 2)T 55 ((++)(++)(++)(++)(−−)(++)) (2, 2, 2, 2, 0, 2)T 56 ((++)(++)(+−)(+−)(++)(++)) (2, 2, 1, 1, 2, 2)T 57 ((++)(++)(−−)(++)(++)(++)) (2, 2, 0, 2, 2, 2)T 58 ((++)(++)(++)(+−)(++)(++)) (2, 2, 2, 1, 2, 2)T 59 ((++)(++)(++)(++)(+−)(++)) (2, 2, 2, 2, 1, 2)T 60 ((++)(++)(+−)(++)(++)(++)) (2, 2, 1, 2, 2, 2)T 61 ((++)(++)(++)(++)(++)(++)) (2, 2, 2, 2, 2, 2)T ∗ denotes the cells corresponding to the integer solutions of Ax
δj 0.496 0.532 0.780 0.708 0.000 0.890 0.849 0.929 0.000 0.708 0.496 0.849 1.205 1.285 1.232 0.708 0.780 0.496 0.890 0.941 0.532 1.259 1.232 1.805 1.632 0.890 0.496 0.941 0.532 0.702 0.843 1.285 1.569 0.532 0.702 0.582 0.860 1.103 1.632 0.843 0.582 0.860 0.828 0.692 1.053 1.863 0.692 1.103 0.803 0.876 1.415 1.784 0.859 0.876 1.863 0.894 1.974 0.901 1.654 1.608 1.578 = b.
<
1
<
< <
1
<
1
1
Fig. 3. Flow chart of finding firing sequence from firing count vector Table 2. Numerical results for randomly generated problems. m 19 19 18 18 24 24 23 23 6 8 19 5 4 19
4. CONCLUSION An integer solution to the fundamental equation must be also a real solution at the same time. Therefore, we confine our search efforts on set C, the set of all real solutions to the fundamental equation, while set C can be partitioned into a finite number of cells by a family of hyperplanes, hij ’s. Furthermore, for each cell, e.g., Cell i,
1
n 20 20 20 20 25 25 25 25 8 10 20 6 6 20
ω 1 2 1 2 1 2 1 2 5 5 5 10 10 10
|X| 220 320 220 320 225 325 225 325 68 610 620 116 116 1120
number of cells 40 47 766 1633 43 65 1157 2487 1021 1631 111 59 1573 160
CPU Time (seconds) 4.3 10.3 113.7 211.1 5.1 11.2 147.4 536.1 153.3 479.9 29.5 10.1 533.3 320.8
it is straightforward to identify its closest integer point, ϕi , in X. Thus, the existence of an integer solution to the fundamental solution on a bounded integer set X reduces to a problem whether or not there is a ϕi that is exactly falling on C (or its distance to C is zero). The efficiency of the solution algorithm developed in this
280
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
paper therefore depends entirely on the efficiency of the cell enumeration scheme which is of a complexity bound O((nω)n−m ). Preliminary numerical results have revealed that the proposed solution method is promising when n−m is small. One of our future research topics is to explore possible ways to improve the efficiency of the algorithm when n − m is large. Our solution method developed in this paper for solving the fundamental equation of Petri nets is applicable in finding integer solution(s) to general linear equations on a finite integer set. In contrast to the rich literature on general linear Diophantine equations in Z n , there exist only a few methods developed for linear Diophantine equations on a bounded integer set. Based on Lov´ asz lattice basis reduction, Aardal et al. (2000) proposed an algorithm to solve bounded linear Diophantine equations. Utilizing the results of Kertzner (1981) to express explicitly integer solutions to a linear equation of two variables, Ramachandran (2006) developed a two-phase algorithm for bounded linear Diophantine equations based on the Euclidean algorithm. All the existing literatures apply only to situations with integral matrix A, while our method proposed in this paper is applicable to situations with real-valued matrix A. REFERENCES
Petri, C.A. (1962). Kommunikation mit Automaten. Ph.D. thesis, Schriften des Rheinisch-Westfalischen Instituts fur Instrumentelle Mathematik, Bonn, Germany. Ramachandran, P. (2006). Use of extended euclidean algorithm in solving a system of linear diophantine equations with bounded variables. Lecture Notes in Computer Science, 4076, 182–192. Schrijver, A. (1986). Theory of Linear and Integer Programming. John Wiley & Sons, New York. Sleumer, N. (1999). Output-sensitive cell enumeration in hyperplane arrangements. Nordic Journal of Computing, 6, 137–161. Sun, X.L., Liu, C.L., Li, D., and Gao, J.J. (2008). On duality gap in binary quadratic optimization. In Working paper series. Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong. Yen, H.C. (2006). Introduction to petri net theory. In Recent Advances in Formal Languages and Applications, Studies in Computational Intelligence, 25, 343–373. Zaslavsky, T. (1975). Facing up to arrangements: Facecount formulas for partitions of space by hyperplanes. Mem. Amer. Math. Soc., 1, 1–101. Zurawski, R. and Zhou, M.C. (1994). Petri nets and industrial applications: A tutorial. IEEE Transactions on Industrial Electronics, 41, 567–583.
Aardal, K., Hurkens, C.A.J., and Lenstra, A.K. (2000). Solving a system of linear diophantine equations with lower and upper bounds on the variables. Mathematics of Operations Research, 25, 427–442. Avis, D. and Fukuda, K. (1996). Reverse search for enumeration. Discrete Applied Mathematics, 65, 21–46. Edelsbrunner, H. (1987). Algorithms in Combinatorial Geometry. Springer. Ferrez, J.A., Fukuda, K., and Liebling, T.M. (2005). Solving the fixed rank convex quadratic maximization in binary variables by a parallel zonotope construction algorithm. European Journal of Operational Research, 166, 35–50. Garey, M.R. and Johnson, D.S. (1979). Computers and Intractability: A Guide to the Theory of NP-completeness. Freeman, San Francisco. Kertzner, S. (1981). The linear diophantine equation. American Mathematical Monthly, 88, 200–203. Kostin, A.E. (2003). Reachability analysis in t-invariantless petri nets. IEEE Transactions on Automatic Control, 48, 1019–1024. Li, D. and Sun, X.L. (2006). Nonlinear Integer Programming. Springer, New York. Matsumoto, T., Miyano, Y., and Jiang, Y. (1997). Some useful sufficient criteria for the basic reachability problem in general petri nets. In Proceedings of the 36th IEEE Conference on Decision and Control, volume 4, 4104–4109. San Diego, CA, USA. Matsumoto, T., Takata, M., and Moro, S. (2002). Reachability analyses in petri nets by groebner bases. In Proceeding of SICE 2002, 841–846. Murata, T. (1977). State equation, controllability, and maximal matchings of petri nets. IEEE Transactions on Automatic Control, 22, 412–416. Nemhauser, G.L. and Wolsey, L. (1988). Integer and Combinatorial Optimization. John Wiley & Sons, New York.
281
A Novel Method for Reachability Determination in Petri Nets ⋆ Duan Li ∗ Xiaoling Sun ∗∗ Jianjun Gao ∗∗∗ Shenshen Gu ∗∗∗∗ Xiaojin Zheng † ∗
Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) ∗∗ Department of Management Science, Fudan University, Shanghai, China (e-mail: [email protected]) ∗∗∗ Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) ∗∗∗∗ Department of Mechanical and Automation Engineering, Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: [email protected]) † Department of Mathematics, Shanghai University, Shanghai, China (e-mail: [email protected]) Abstract: Reachability is one of the most important behavioral properties of Petri nets and the past four decades have witnessed great efforts in developing various implementable methodologies in determining reachability of Petri nets. We propose in this paper a novel method for solving the fundamental equation in the reachability analysis, which has been known to be NP-complete. More specifically, by adopting a revised version of the cell enumeration method for an arrangement of hyperplanes in discrete geometry, we develop an efficient solution scheme to identify firing count vector solution(s) to the fundamental equation on a bounded integer set, with a complexity bound of O((nω)n−m ), where n is the number of transitions, m is the number of places and ω is the upper bound of the number of firings for every transition. Keywords: Petri nets, reachability analysis, cell enumeration, linear Diophantine equations on bounded integer set, integer programming, discrete optimization. 1. INTRODUCTION Petri net, introduced by Petri in his seminal work “Communication with Automata” (Petri (1962)), has been a promising mathematical formalism for modeling, analyzing and designing discrete event systems, especially by its remarkable capability in modeling process synchronization, asynchronous events, concurrent and distributed operations and conflicts or resource sharing. The past four decades have witnessed innumerable successful applications of Petri nets in various areas, such as i) modeling and analysis of communication protocols and networks, production systems, sequence controllers and software development, and ii) performance evaluation of the modeled systems. The Petri net is a particular kind of bipartite directed graph consisting of three types of elements: places, transitions, and directed arcs connecting places and transitions, while each place may hold a non-negative number of tokens. See Fig. 1 1 for an illustrative example of Petri nets with 4 places and 6 transitions. The distribution of tokens ⋆ This research is partially supported by University Grants Council, Hong Kong, under CUHK414207, and Shanghai Leading Academic Discipline Project (B210).
978-3-902661-43-2/09/$20.00 © 2009 IFAC
Fig. 1. Example of a Petri net on places, called Petri net marking, defines the state of the modeled system. A marking for a Petri net with m places is represented by an (m × 1) vector M, where Mj , j = 1, 2, . . . , m, are nonnegative integers representing the number of tokens in the corresponding places. A general Petri net is characterized by a quintuple (P, T, I, O, M0 ) (see Zurawski and Zhou (1994)), where P = {p1 , p2 , . . . , pm } is a finite set of m places, T = {t1 , t2 , . . . , tn } is a finite set of n transitions, DI : (P × m×n T ) → Z+ is an input function that defines weights
276
10.3182/20090603-3-RU-2001.0100
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
associated with directed arcs from places to transitions, m×n m ) is the set of (m × n) dimensional mawhere Z+ (Z+ trices (m dimensional vector) with all entries being in Z+ , m×n DO : (P × T ) → Z+ is an output function that defines weights associated with directed arcs from transitions to m places, and M0 : P → Z+ is the initial marking. For the Petri net given in Figure 1, we have 0 1 0 0 1 0 0 0 1 0 0 1 DI = , 0 0 0 4 0 0 2 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 , DO = 0 0 1 0 1 0 0 0 0 1 0 5 T and M0 = (1, 1, 0, 0) . The change of the distribution of the tokens represents the dynamics of the modeled system, while the distribution of tokens on places may change according to the following enabling rule and firing rule (see Zurawski and Zhou (1994)): Enabling Rule: A transition t is said to be enabled if each input place p of t contains at least the number of tokens equal to the weight of the arc connecting p to t. Firing Rule: (a) An enabled transition t may or may not fire depending on the additional interpretation, and (b) A firing of an enabled transition t removes form each input place p the number of tokens equal to the weight of the directed arc connecting p to t, and deposits in each output place p the number of tokens equal to the weight of the arc connecting t to p. Reachability is one of the most important behavioral and structural properties of Petri nets. Given both an initial state M0 and a target terminal state M , a natural question to ask is whether or not we have a sequence of firing rules such that the system can reach the specific target state within finite steps. There are two primary approaches in investigating reachability: Reachability graph analysis and reachability algebraic analysis. Reachability graph analysis is based on the creation and investigation of a reachability graph or a reduced reachability graph. As this approach relies on exhaustively generating all the reachable markings from a given initial marking, it suffers the state explosion phenomenon, while transferring the reachability graph to a reduced counterpart is NP-hard (see Kostin (2003)). The second approach is based on methods of linear algebra. It is well known (see Murata (1977) and Yen (2006)) that a necessary condition for reachability of marking M from an initial marking M0 of a Petri net is that there exists a nonnegative integer vector solution of the following system of linear equations, Ax = b,
(1) T
where b = M − M0 , x = (x1 , x2 , . . . , xm ) is a firing count vector and A is an (m × n) dimensional incidence matrix given by A = [a(pi , tj )] = [DO (pi , tj ) − DI (pi , tj )]. If we set in Figure 1 the final marking to be M = (2, 0, 1, 3)T , the fundamental equation for this example is then given by
2 0 0 −2
−1 1 0 0
0 −1 1 0
0 0 −4 1
−1 0 1 0
0 1 −1 −1 x= . 0 1 5 3
When the Petri net is acyclic, the above condition becomes both necessary and sufficient (Matsumoto et al. (1997)). Equation (1), termed the fundamental equation in Petri nets, reveals that the investigation of reachability can be achieved, to certain degree, by finding out nonnegative integer solution(s) to a system of linear equations. The first step in the reachability algebraic analysis is to find firing count vectors, which can be performed either by solving the fundamental equation directly, or by solving the corresponding integer programming problem. Unfortunately, finding nonnegative integer solutions to the fundamental equation in Petri nets is NP-complete (Matsumoto et al. (2002)). The second step in the reachability algebraic analysis is to translate a firing count vector into a firing sequence, if there is a reachable path, while such an algorithm (see, e.g., Kostin (2003)) is of a combinatorial nature. In essence, solving the fundamental equation is an integer programming problem (see, for example, Schrijver (1986) Nemhauser and Wolsey (1988) and Li and Sun (2006)). Generally speaking, the capability of the existing methods for finding the nonnegative integer solution to the fundamental equation (1) is still very limited. In particular, the Petri nets community has been largely dependent on some existing softwares to obtain nonnegative integer solutions to the fundamental equation. More efficient solution schemes that invoke the state-of-the-art from some other related subjects, especially from discrete optimization and integer programming, deserve further investigation. There have been a few attempts in the Petri net communities in this direction. For example, Matsumoto et al. (2002) applied Groebner Bases method from integer programming to solve reachability problems in Petri nets. The solution approach using Groebner Bases depends heavily on symbolic computing power, and its capability in solving large-scale problem is thus not promising. Although the original analytical reachability analysis is to find nonnegative integer solution(s) to the fundamental equation, a more practical problem is to find integer solution(s) to the fundamental equation on a bounded integer set, X = {x ∈ Z n | 0 ≤ xi ≤ ui , i = 1, . . . , n}, where Z n denotes the set of all integer points in Rn . This confinement is reasonable as there always exist resource constraints in real-world applications of Petri nets. For example, in chemical processes, the operation of reactors (transitions in a corresponding Petri net) is expensive, thus the number of such operations should not exceed a given upper bound. It is well known that linear Diophantine equations, Ax = b, x ∈ Z n , are polynomially solvable (see, e.g., Schrijver (1986) and Nemhauser and Wolsey (1988)), while linear Diophantine equations on a bounded integer set, Ax = b, x ∈ X, where A is an m × n integral matrix with rank(A) = m, and b ∈ Rm being integral, are NP-complete, as the special case of linear Diophantine equations with m = 1, x ∈ {−1, 1}n and b = 0 is NPcomplete (Garey and Johnson (1979)).
277
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
We propose in this paper a novel method for solving the fundamental equation in the reachability analysis over a bounded integer set X of firing count vectors. Our method is based on a recognition that whether or not there exists a solution to the fundamental solution is equivalent to whether or not the distance from X to the affine solution set of Ax = b in Rn attains zero. Furthermore, from recent results in Sun et al. (2008), finding the distance from an integer set to an affine set can be efficiently achieved by the cell enumeration method for an arrangement of hyperplanes in discrete geometry. This cell enumeration approach provides a promising platform for designing an efficient method with a complexity of O((nω)n−m ), where ω = maxi=1,...,n ui , to find integer solutions to the fundamental equation on a bounded integer set . 2. SOLVING THE FUNDAMENTAL EQUATION BY CELL ENUMERATION Let U = (U1 , . . . , Un−m ) be an orthogonal basis for the null space of A and x0 be a special solution to Ax = b in Rn , then C = {x ∈ Rn | Ax = b} can be expressed by n
C = {x ∈ R | x = x0 +
n−m
(n − m)-dimensional polyhedral set formed by the halfspaces induced by hij ’s and can be characterized by a n 1 n u j=1 j -dimensional sign vector, sign(E) = (w , . . . , w ), where wj = (w1,j , . . . , wuj ,j ) is specified by
+, if gij (π) > 0 wij = (5) , i = 1, . . . , uj , −, if gij (π) < 0 with π being an interior point of E. For fixed j, since hij (i = 1, . . . , uj ) are parallel hyperplanes, wj = (w1,j , . . . , wuj ,j ) must take the following form: i
(+, . . . , +, −, . . . , −),
i = 0, 1, . . . , uj .
The one-to-one map between all distinct ϕ(y) for y ∈ C and the sign vectors of all the cells of the hyperplane arrangement can be established as follows: w = (w1 , . . . , wn ) ←→ ϕ = (ϕ1 , . . . , ϕn ), where for j = 1, . . . , n, wj = (−, . . . , −) ←→ ϕj = 0, wj = (+, . . . , +) ←→ ϕj = uj ,
zi Ui ,
i
w = (+, . . . , +, −, . . . , −) ←→ ϕj = i, j
i=1
zi ∈ R, i = 1, . . . , n − m}.
(2)
Consider the Euclidian distance from C to X = {x ∈ Z n | 0 ≤ xi ≤ ui , i = 1, . . . , n}: δ = dist(C, X) = min{x − y | x ∈ X, y ∈ C}.
(3)
It is obvious that Ax = b has a solution in X if and only if δ = 0. Furthermore, when δ = 0, any x∗ ∈ X that achieves the zero distance is an integer solution to the linear equations Ax = b in X. Note that combining (2) and (3) yields
where i ∈ {1, . . . , uj − 1}.
It has been known that the number of cells of the hyperplane arrangement generated by (4) is bounded by O((nω)n−m ), where ω = maxi=1,...,n ui (see Zaslavsky (1975) and Edelsbrunner (1987)). Moreover, using the cell enumeration methods in Avis and Fukuda (1996) and Sleumer (1999), we can find all the cells of the hyperplane arrangement generated by (4) in O((nω)n−m ) time. Listing all the distinct ϕ(y) for y ∈ C as ϕ1 , . . . , ϕk , the distance δ = dist(C, X) can then be calculated as follows: δ = min δj ≡ min dist(C, ϕj )
δ = dist(C, X)
j=1,...,k
= min min{y − x | 0 ≤ xi ≤ ui , xi ∈ Z, i = 1, . . . , n}
= min (U U T − In )(ϕj − x0 ).
y∈C
j=1,...,k
= min y − ϕ(y), y∈C
where, for j = 1, . . . , n, ϕ(y)j is determined by the following equation: 1 0, yj ∈ (−∞, ], 2 1 1 ϕ(y)j = i, yj ∈ (i − , i + ), i ∈ {1, . . . , uj − 1}, 2 2 uj , yj ∈ [uj − 1 , ∞). 2
Let’s consider the hyperplane arrangement generated by n the following j=1 uj hyperplanes in Rn−m :
+
l=1
1 zl Ulj − (i − ) = 0} 2
(6)
Therefore, finding a solution of the linear system Ax = b over the integer set X or checking its infeasibility can be done in O((nω)n−m ) time. Especially, when r = n − m is fixed (0 ≤ r ≤ n − 1), the linear system Ax = b over X is polynomially solvable. Let us consider now the example in Fig. 1 to illustrate the solution algorithm for the fundamental equation using the cell enumeration method. For this instance, we assume that the bounded integer set is given by X = {x ∈ Z 6 | 0 ≤ xi ≤ 2, i = 1, . . . , 6}. Since the incidence matrix is of rank 4 in this example, solving the fundamental equation in R6 yields
hij = {z ∈ Rn−m | gij (z) := (x0 )j n−m
j=1,...,k
C = {x ∈ R6 | Ax = b} (4)
for i = 1, . . . , uj , j = 1, . . . , n. Note that for each j, hij (i = 1, . . . , uj ) are parallel hyperplanes. A cell E of the hyperplane arrangement corresponding to hij ’s is an
= {x ∈ R6 | x = x0 + z1 U1 + z2 U2 , zi ∈ R, i = 1, 2}, where x0 = (0.2912, −0.2138, 0.0100, −0.2984, −0.2038, 0.776)T and
278
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
δ5 = dist(X, Cell 5) = (U U T − I6 )(ϕ5 − x0 ) = 0.
15
g14 g26
For Cell 61 labelled in Figure 2, all points are above the 12 hyperplanes. Thus, the sign vector for Cell 61 is ((++)(++)(++)(++)(++)(++)) and the closest integer point in {0, 1, 2}6 to Cell 61 is ϕ61 = (2, 2, 2, 2, 2, 2). Furthermore, the distance between X and Cell 61 is given by
g24 Cell 61
10
g21 g15
g25
z2
g11 5
Cell 9
g16
δ61 = dist(X, Cell 61) = (U U T − I6 )(ϕ61 − x0 ) = 1.578.
g22 g12 0
Cell 5 g23 g13 −5 −10
−8
−6
−4
−2
0
z1
2
4
6
8
10
Fig. 2. Cell arrangement for the example problem. 0.3508 0.3944 −0.1605 0.6882 −0.2713 0.5636 U = (U1 , U2 ) = . 0.1477 0.1661 0.8622 0.1006 0.1108 0.1245 Note that uj = 2, j = 1, . . ., 6. Setting each xj , j = 1, . . ., 6, equal to 12 and 32 , respectively, generates the following 12 hyperplanes on C defined in (4),
g11 (z) = 0.3508z1 + 0.3944z2 = 0.2088, j=1: g21 (z) = 0.3508z1 + 0.3944z2 = 1.2088,
g12 (z) = −0.1605z1 + 0.6882z2 = 0.7138, j=2: g22 (z) = −0.1605z1 + 0.6882z2 = 1.7138,
g13 (z) = −0.2713z1 + 0.5636z2 = 0.4900, j=3: g23 (z) = −0.2713z1 + 0.5636z2 = 1.4900,
g14 (z) = 0.1477z1 + 0.1661z2 = 0.7984, j=4: g24 (z) = 0.1477z1 + 0.1661z2 = 1.7984,
g15 (z) = 0.8622z1 + 0.1006z2 = 0.7038, j=5: g25 (z) = 0.8622z1 + 0.1006z2 = 1.7038,
g16 (z) = 0.1108z1 + 0.1245z2 = −0.2762, j=6: g26 (z) = 0.1108z1 + 0.1245z2 = 0.7238. As illustrated in Fig. 2, these 12 hyperplanes partition set C into various cells, each of which corresponds to a unique sign vector defined in (5), sign(E) = (w1 , . . . , w6 ), where wj = (w1j , w2j ). Note that the arrow attached to each hyperplane points to the positive half space generated by this hyperplane. Table 1 lists all the 61 cells identified with their sign vectors and the corresponding distances to X = {0, 1, 2}6 . Let us consider Cells 5 and 61 labeled in Fig. 2. As any interior point π of Cell 5 satisfies g11 (π) > 0, g21 (π) < 0, g12 (π) < 0, g22 (π) < 0, g13 (π) < 0, g23 (π) < 0, g14 (π) < 0, g24 (π) < 0, g15 (π) > 0, g25 (π) < 0, g16 (π) > 0, and g26 (π) < 0, the sign vector for any interior point π of Cell 5 is thus specified by ((+−)(−−)(−−)(−−)(+−)(+−)). Furthermore, for all points in Cell 5, the closest integer point in X = {0, 1, 2}6 is clearly given by ϕ5 = (1, 0, 0, 0, 1, 1) and the distance between X and Cell 5 is given by
We observe from Table 1 that there exist two solutions, (1, 0, 0, 0, 1, 1)T and (1, 1, 1, 0, 0, 1)T , to the fundamental equation of this example, as they both achieve the zero distance. We need to emphasize that our solution scheme enables us to identify all solutions to the fundamental equation on a finite integer set. Note also that the total number of cells in this example is 61, which is much smaller than the upper bound of the cell numbers O((6 × 2)(6−4) ). This conclusion applies for general situations as observed from randomly generated numerical experiments presented in the next section. The flow chart given in Fig. 3 demonstrates the structure of an algorithm in converting a (feasible) firing count vector into a firing sequence, which is similar to the one in Kostin (2003). For the current instance, firing count vector (1, 0, 0, 0, 1, 1)T leads to a firing sequence of {t5 → t6 → t1 } and firing count vector (1, 1, 1, 0, 0, 1)T leads to a firing sequence {t2 → t3 → t6 → t1 }. As both firing count vectors are realizable, a decision maker would prefer firing count vector (1, 0, 0, 0, 1, 1)T to firing count vector (1, 1, 1, 0, 0, 1)T as the former one requires less number of firings, and thus smaller expenses (when assuming the same cost for firings at different transitions).
3. PRELIMINARY NUMERICAL RESULTS We present some preliminary numerical results for the solution method using cell enumeration developed in Section 2. The solution algorithm was coded by C and ran on a Pentium PC (2G Hz and 1G RAM). Our implementation of the cell enumeration method is based on the reverse search procedure proposed in Avis and Fukuda (1996) and Sleumer (1999). In our implementation, the LP solver in CPLEX 9.0 is called to solve the linear programming subproblems in the reverse search procedure. In all the test problems, we set ω = ui , i = 1, . . . , n. For each triple of (m, n, ω), 10 test problems are randomly generated. The numerical results are reported in Table 2, where |X| denotes the total number of integer points in X, the number of cells is the average for 10 test problems (after rounding off) and the CPU time is the average for 10 test problems. We point out that the purpose of implementing our proposed solution method in numerical experiments at this stage is to check the feasibility of the proposed solution method. Sophisticated implementation of the cell enumeration procedure that integrates the reverse searching technique with ZRAM library (see Ferrez et al. (2005)) will further significantly improve the performance of the method.
279
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
Table 1. A full list of cells, their sign vectors and the closest integer points. Cell j (w1 , w2 , w3 , w4 , w5 , w6 ) ϕj 1 ((+−)(−−)(−−)(−−)(−−)(+−)) (1, 0, 0, 0, 0, 1)T 2 ((−−)(−−)(−−)(−−)(−−)(+−)) (0, 0, 0, 0, 0, 1)T 3 ((+−)(+−)(−−)(−−)(−−)(+−)) (1, 1, 0, 0, 0, 1)T 4 ((+−)(−−)(+−)(−−)(−−)(+−)) (1, 0, 1, 0, 0, 1)T 5∗ ((+−)(−−)(−−)(−−)(+−)(+−)) (1, 0, 0, 0, 1, 1)T 6 ((−−)(−−)(+−)(−−)(−−)(+−)) (0, 0, 1, 0, 0, 1)T 7 ((−−)(−−)(−−)(−−)(+−)(+−)) (0, 0, 0, 0, 1, 1)T 8 ((−−)(−−)(−−)(−−)(−−)(−−)) (0, 0, 0, 0, 0, 0)T ((+−)(+−)(+−)(−−)(−−)(+−)) (1, 1, 1, 0, 0, 1)T 9∗ 10 ((+−)(+−)(−−)(−−)(+−)(+−)) (1, 1, 0, 0, 1, 1)T 11 ((+−)(−−)(−−)(−−)(++)(+−)) (1, 0, 0, 0, 2, 1)T 12 ((−−)(+−)(+−)(−−)(−−)(+−)) (0, 1, 1, 0, 0, 1)T 13 ((−−)(−−)(+−)(−−)(−−)(−−)) (0, 0, 1, 0, 0, 0)T 14 ((−−)(−−)(−−)(−−)(++)(+−)) (0, 0, 0, 0, 2, 1)T 15 ((−−)(−−)(−−)(−−)(+−)(−−)) (0, 0, 0, 0, 1, 0)T 16 ((+−)(++)(+−)(−−)(−−)(+−)) (1, 2, 1, 0, 0, 1)T 17 ((+−)(+−)(++)(−−)(−−)(+−)) (1, 1, 2, 0, 0, 1)T 18 ((+−)(+−)(+−)(−−)(+−)(+−)) (1, 1, 1, 0, 1, 1)T 19 ((++)(+−)(−−)(−−)(+−)(+−)) (2, 1, 0, 0, 1, 1)T 20 ((+−)(+−)(−−)(−−)(++)(+−)) (1, 1, 0, 0, 2, 1)T 21 ((++)(−−)(−−)(−−)(++)(+−)) (2, 0, 0, 0, 2, 1)T 22 ((−−)(+−)(++)(−−)(−−)(+−)) (0, 1, 2, 0, 0, 1)T 23 ((−−)(+−)(+−)(−−)(−−)(−−)) (0, 1, 1, 0, 0, 0)T 24 ((−−)(−−)(++)(−−)(−−)(−−)) (0, 0, 2, 0, 0, 0)T 25 ((−−)(−−)(−−)(−−)(++)(−−)) (0, 0, 0, 0, 2, 0)T 26 ((++)(++)(+−)(−−)(−−)(+−)) (2, 2, 1, 0, 0, 1)T 27 ((+−)(++)(++)(−−)(−−)(+−)) (1, 2, 2, 0, 0, 1)T 28 ((+−)(++)(+−)(−−)(+−)(+−)) (1, 2, 1, 0, 1, 1)T 29 ((++)(+−)(+−)(−−)(+−)(+−)) (2, 1, 1, 0, 1, 1)T 30 ((++)(+−)(−−)(−−)(++)(+−)) (2, 1, 0, 0, 2, 1)T 31 ((++)(−−)(−−)(+−)(++)(+−)) (2, 0, 0, 1, 2, 1)T 32 ((−−)(++)(++)(−−)(−−)(+−)) (0, 2, 2, 0, 0, 1)T 33 ((−−)(+−)(++)(−−)(−−)(−−)) (0, 1, 2, 0, 0, 0)T 34 ((++)(++)(++)(−−)(−−)(+−)) (2, 2, 2, 0, 0, 1)T 35 ((++)(++)(+−)(−−)(+−)(+−)) (2, 2, 1, 0, 1, 1)T 36 ((++)(+−)(+−)(−−)(++)(+−)) (2, 1, 1, 0, 2, 1)T 37 ((++)(+−)(−−)(+−)(++)(+−)) (2, 1, 0, 1, 2, 1)T 38 ((++)(−−)(−−)(+−)(++)(++)) (2, 0, 0, 1, 2, 2)T 39 ((−−)(++)(++)(−−)(−−)(−−)) (0, 2, 2, 0, 0, 0)T 40 ((++)(++)(++)(+−)(−−)(+−)) (2, 2, 2, 1, 0, 1)T 41 ((++)(++)(++)(−−)(+−)(+−)) (2, 2, 2, 0, 1, 1)T 42 ((++)(++)(+−)(+−)(+−)(+−)) (2, 2, 1, 1, 1, 1)T 43 ((++)(++)(+−)(−−)(++)(+−)) (2, 2, 1, 0, 2, 1)T 44 ((++)(+−)(+−)(+−)(++)(+−)) (2, 1, 1, 1, 2, 1)T 45 ((++)(+−)(−−)(+−)(++)(++)) (2, 1, 0, 1, 2, 2)T 46 ((++)(−−)(−−)(++)(++)(++)) (2, 0, 0, 2, 2, 2)T 47 ((++)(++)(++)(+−)(+−)(+−)) (2, 2, 2, 1, 1, 1)T 48 ((++)(++)(++)(+−)(−−)(++)) (2, 2, 2, 1, 0, 2)T 49 ((++)(++)(+−)(+−)(++)(+−)) (2, 2, 1, 1, 2, 1)T 50 ((++)(+−)(+−)(+−)(++)(++)) (2, 1, 1, 1, 2, 2)T 51 ((++)(++)(−−)(+−)(++)(++)) (2, 2, 0, 1, 2, 2)T 52 ((++)(+−)(−−)(++)(++)(++)) (2, 1, 0, 2, 2, 2)T 53 ((++)(++)(++)(+−)(++)(+−)) (2, 2, 2, 1, 2, 1)T 54 ((++)(++)(++)(+−)(+−)(++)) (2, 2, 2, 1, 1, 2)T 55 ((++)(++)(++)(++)(−−)(++)) (2, 2, 2, 2, 0, 2)T 56 ((++)(++)(+−)(+−)(++)(++)) (2, 2, 1, 1, 2, 2)T 57 ((++)(++)(−−)(++)(++)(++)) (2, 2, 0, 2, 2, 2)T 58 ((++)(++)(++)(+−)(++)(++)) (2, 2, 2, 1, 2, 2)T 59 ((++)(++)(++)(++)(+−)(++)) (2, 2, 2, 2, 1, 2)T 60 ((++)(++)(+−)(++)(++)(++)) (2, 2, 1, 2, 2, 2)T 61 ((++)(++)(++)(++)(++)(++)) (2, 2, 2, 2, 2, 2)T ∗ denotes the cells corresponding to the integer solutions of Ax
δj 0.496 0.532 0.780 0.708 0.000 0.890 0.849 0.929 0.000 0.708 0.496 0.849 1.205 1.285 1.232 0.708 0.780 0.496 0.890 0.941 0.532 1.259 1.232 1.805 1.632 0.890 0.496 0.941 0.532 0.702 0.843 1.285 1.569 0.532 0.702 0.582 0.860 1.103 1.632 0.843 0.582 0.860 0.828 0.692 1.053 1.863 0.692 1.103 0.803 0.876 1.415 1.784 0.859 0.876 1.863 0.894 1.974 0.901 1.654 1.608 1.578 = b.
<
1
<
< <
1
<
1
1
Fig. 3. Flow chart of finding firing sequence from firing count vector Table 2. Numerical results for randomly generated problems. m 19 19 18 18 24 24 23 23 6 8 19 5 4 19
4. CONCLUSION An integer solution to the fundamental equation must be also a real solution at the same time. Therefore, we confine our search efforts on set C, the set of all real solutions to the fundamental equation, while set C can be partitioned into a finite number of cells by a family of hyperplanes, hij ’s. Furthermore, for each cell, e.g., Cell i,
1
n 20 20 20 20 25 25 25 25 8 10 20 6 6 20
ω 1 2 1 2 1 2 1 2 5 5 5 10 10 10
|X| 220 320 220 320 225 325 225 325 68 610 620 116 116 1120
number of cells 40 47 766 1633 43 65 1157 2487 1021 1631 111 59 1573 160
CPU Time (seconds) 4.3 10.3 113.7 211.1 5.1 11.2 147.4 536.1 153.3 479.9 29.5 10.1 533.3 320.8
it is straightforward to identify its closest integer point, ϕi , in X. Thus, the existence of an integer solution to the fundamental solution on a bounded integer set X reduces to a problem whether or not there is a ϕi that is exactly falling on C (or its distance to C is zero). The efficiency of the solution algorithm developed in this
280
13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009
paper therefore depends entirely on the efficiency of the cell enumeration scheme which is of a complexity bound O((nω)n−m ). Preliminary numerical results have revealed that the proposed solution method is promising when n−m is small. One of our future research topics is to explore possible ways to improve the efficiency of the algorithm when n − m is large. Our solution method developed in this paper for solving the fundamental equation of Petri nets is applicable in finding integer solution(s) to general linear equations on a finite integer set. In contrast to the rich literature on general linear Diophantine equations in Z n , there exist only a few methods developed for linear Diophantine equations on a bounded integer set. Based on Lov´ asz lattice basis reduction, Aardal et al. (2000) proposed an algorithm to solve bounded linear Diophantine equations. Utilizing the results of Kertzner (1981) to express explicitly integer solutions to a linear equation of two variables, Ramachandran (2006) developed a two-phase algorithm for bounded linear Diophantine equations based on the Euclidean algorithm. All the existing literatures apply only to situations with integral matrix A, while our method proposed in this paper is applicable to situations with real-valued matrix A. REFERENCES
Petri, C.A. (1962). Kommunikation mit Automaten. Ph.D. thesis, Schriften des Rheinisch-Westfalischen Instituts fur Instrumentelle Mathematik, Bonn, Germany. Ramachandran, P. (2006). Use of extended euclidean algorithm in solving a system of linear diophantine equations with bounded variables. Lecture Notes in Computer Science, 4076, 182–192. Schrijver, A. (1986). Theory of Linear and Integer Programming. John Wiley & Sons, New York. Sleumer, N. (1999). Output-sensitive cell enumeration in hyperplane arrangements. Nordic Journal of Computing, 6, 137–161. Sun, X.L., Liu, C.L., Li, D., and Gao, J.J. (2008). On duality gap in binary quadratic optimization. In Working paper series. Department of Systems Engineering and Engineering Management, Chinese University of Hong Kong. Yen, H.C. (2006). Introduction to petri net theory. In Recent Advances in Formal Languages and Applications, Studies in Computational Intelligence, 25, 343–373. Zaslavsky, T. (1975). Facing up to arrangements: Facecount formulas for partitions of space by hyperplanes. Mem. Amer. Math. Soc., 1, 1–101. Zurawski, R. and Zhou, M.C. (1994). Petri nets and industrial applications: A tutorial. IEEE Transactions on Industrial Electronics, 41, 567–583.
Aardal, K., Hurkens, C.A.J., and Lenstra, A.K. (2000). Solving a system of linear diophantine equations with lower and upper bounds on the variables. Mathematics of Operations Research, 25, 427–442. Avis, D. and Fukuda, K. (1996). Reverse search for enumeration. Discrete Applied Mathematics, 65, 21–46. Edelsbrunner, H. (1987). Algorithms in Combinatorial Geometry. Springer. Ferrez, J.A., Fukuda, K., and Liebling, T.M. (2005). Solving the fixed rank convex quadratic maximization in binary variables by a parallel zonotope construction algorithm. European Journal of Operational Research, 166, 35–50. Garey, M.R. and Johnson, D.S. (1979). Computers and Intractability: A Guide to the Theory of NP-completeness. Freeman, San Francisco. Kertzner, S. (1981). The linear diophantine equation. American Mathematical Monthly, 88, 200–203. Kostin, A.E. (2003). Reachability analysis in t-invariantless petri nets. IEEE Transactions on Automatic Control, 48, 1019–1024. Li, D. and Sun, X.L. (2006). Nonlinear Integer Programming. Springer, New York. Matsumoto, T., Miyano, Y., and Jiang, Y. (1997). Some useful sufficient criteria for the basic reachability problem in general petri nets. In Proceedings of the 36th IEEE Conference on Decision and Control, volume 4, 4104–4109. San Diego, CA, USA. Matsumoto, T., Takata, M., and Moro, S. (2002). Reachability analyses in petri nets by groebner bases. In Proceeding of SICE 2002, 841–846. Murata, T. (1977). State equation, controllability, and maximal matchings of petri nets. IEEE Transactions on Automatic Control, 22, 412–416. Nemhauser, G.L. and Wolsey, L. (1988). Integer and Combinatorial Optimization. John Wiley & Sons, New York.
281