Model abstraction for discrete-event systems by binary linear programming with applications to manufacturing systems

Model abstraction for finite state automata is helpful for decreasing computational complexity and improving comprehensibility for the verification and control synthesis of discrete-event systems (DES). Supremal quasi-congruence equivalence is an effective method for reducing the state space of DES and its effective algorithms based on graph theory have been developed. In this paper, a new method is proposed to convert the supremal quasi-congruence computation into a binary linear programming problem which can be solved by many powerful integer linear programming and satisfiability (SAT) solvers. Partitioning states to cosets is considered as allocating states to an unknown number of cosets and the requirement of finding the coarsest quasi-congruence is equivalent to using the least number of cosets. The novelty of this paper is to solve the optimal partitioning problem as an optimal state-to-coset allocation problem. The task of finding the coarsest quasi-congruence is equivalent to the objective of finding the least number of cosets. Then the problem can be solved by optimization methods, which are respectively implemented by mixed integer linear programming (MILP) in MATLAB and binary linear programming (BLP) in CPLEX. To reduce the computation time, the translation process is first optimized by introducing fewer decision variables and simplifying constraints in the programming problem. Second, the translation process formulates a few techniques of converting logic constraints on finite automata into binary linear constraints. These techniques will be helpful for other researchers exploiting integer linear programming and SAT solvers for solving partitioning or grouping problems. Third, the computational efficiency and correctness of the proposed method are verified by two different solvers. The proposed model abstraction approach is applied to simplify the large-scale supervisor model of a manufacturing system with five automated guided vehicles. The proposed method is not only a new solution for the coarsest quasi-congruence computation, but also provides us a more intuitive understanding of the quasi-congruence relation in the supervisory control theory. A future research direction is to apply more computationally efficient solvers to compute the optimal state-to-coset allocation problem.


Introduction
Discrete-event (dynamic) systems are a typical class of computer-integrated manmade system with discrete state spaces and event-driven characteristics. They can be modeled by a state-transition structure such as deterministic or nondeterministic finite automata (DFAs or NFAs) and Petri nets. Both formalisms are used for designing supervisory controllers of a DES. [1][2][3][4][5] On the one hand, the supervisory control theory developed in Wonham and Cai, 1 Ramadge and Wonham 3 is regarded as a seminal work for the treatment of supervisor synthesis for a DES such that the closed-loop behavior is observable, 6,7 controllable, nonblocking, 8 and maximally permissive. On the other hand, a Petri net (PN) is also a common mathematical tool for modeling and analyzing a DES, which often avoids explicit state exploration but derives certain network properties through the net structure. [9][10][11][12][13][14] Both methods are popular for synthesizing satisfactory supervisors.
Due to the theoretical maturity, DES have received much attention from academy and industry communities. [15][16][17][18][19][20][21][22] However, DES control methods have not yet been extensively applied to industrial scale problems, owing to the state explosion problem: the number of reachable states increases exponentially with that of concurrent components and variables. 23 The intractable number of states may quickly deplete the computer memory and terminate the computation if the given system is structurally complex. Therefore, effectively exploring the huge state space with limited memory and time is a significant and persistent research subject.
To tackle the state explosion problem, many reduction methods have been developed, including equivalence reduction 24,25 and process equivalence, 26 model abstraction, 27 symbolic computation and state tree structures, 28,29 and partial order reduction, 30 which all benefit from the proper design and analysis of the system architecture.
Model abstraction based on equivalence relations is an efficient way to reduce the state size of finite state automata. An equivalence relation on the state set of a finite state automaton partitions the state set as a collection of equivalence classes (or cosets). The set of cosets defines a simpler finite state automaton with fewer states, which is called the quotient automaton 31,32 of the original one with respect to the equivalence relation. In many applications, the observable language of the quotient automaton must be equivalent to the original automaton. The identity implies that if two states of the original automaton are equivalent, the same observable event sequences must be defined at both states. This property is guaranteed by the quasi-congruence equivalence, 1 which is similar to the concept of bisimulation. 33 It is necessary to thoroughly investigate the results on the quasi-congruence relation of DES in references, 1,2,34,35 and use them to effectively reduce the computational complexity of the formal verification and control synthesis of DES. Most contemporary approaches 31,32,[36][37][38][39] are based on graph theory. However, we adopt a new approach to the quasi-congruence relations by binary linear programming in the paper, and apply it to reduce the state space.
Wong and Wonham 31 solve the problem on how to design an observer by modifying the given causal reporter map that is not an observer. The method is to compute the coarsest equivalence relation, which is finer than that of the given map. A polynomial-time algorithm is developed to compute an observer of an automaton, and it can also be used to verify the observer property of natural projection. However, this algorithm cannot be applied to natural projection directly.
Feng and Wonham 32 adopt the algorithm in Wong and Wonham 31 to the natural projection, where the minimum event extension is particularly interesting, whose computation, however, is shown to be NP-hard. In this case, an acceptable event extension can be found by developing a polynomial-time algorithm.
These methods are based on graph theory and state-transition data structure. The novelty of this paper is to solve the optimal partitioning problem as an optimal state-to-coset allocation problem. Partitioning states to cosets is considered as allocating states to an unknown number of cosets and the requirement of finding the coarsest quasi-congruence is equivalent to using the least number of cosets. Then the problem can be solved by optimization methods, which leverage many powerful and easy-to-use software tools, such as CPLEX (the CPLEX optimizer can be available on the website https://www.ibm.com/analytics/cplex-optimizer) 40 and SAT solvers (the SAT solvers can be available on the website https://www.satlive.org/solvers). 41 These computation tools exploit modern computing technologies, such as graphics processing unit (GPU) computing, 42 parallel computing, 43 and cloud computing. 44 Large scale problems of finding the coarsest quasi-congruence may be solved by modern optimization technologies.
To this end, this paper formalizes the problem of finding the coarsest quasicongruence as a binary linear programming problem. The design variables form a binary allocation matrix between the states and cosets. The constraints are relationships between states and cosets, and the objective function is to minimize the number of cosets. The main contributions of this paper are summarized below.
1. The number of design variables determines the computational complexity.
In Section 4, a property of the allocation problem is discovered and proved to remove around half of the binary design variables. 2. A novel approach on how to compute the supremal quasi-congruence is proposed, which is based on a combination of propositional logic and binary representation in Section 5.
3. The method proposed in Section 5 is converted to a binary linear programming problem, so as to save computing time. In this transformation, we formulate a few techniques of converting logic constraints on finite automata into binary linear constraints which are formalized and proved to be correct in Section 6. 4. The correctness and efficiency of our method are verified by different solvers and examples in Section 8.
Compared with existing studies, the motivation of this work is to propose an optimization method for quasi-congruence computation of DFAs, and to increase the generality of optimal quasi-congruence computation algorithms. The computational results of this method are the same as that of using the computation software TCT (the software TCT can be available on the website https://www.control.utoronto.ca/DES/Research). 45 for supervisory control design.
The work is structured as follows. The preliminaries pertaining to automata and natural projection are recalled in the following section. Some functions and concepts that are necessary for the quasi-congruence calculation, and the coarsest quasi-congruence are reported in Section 3. Section 7 touches on the main algorithm implementation of the proposed approach. Finally, Section 9 concludes this research.

Preliminaries
This section recalls some basic concepts of deterministic finite automata used in the paper. Then we review the definition of quasi-congruence relation and some operations on it. More details on automata can be found in references. 1,2,31,32,[46][47][48] Fundamental of automata A deterministic finite automaton (DFA) 1 is represented as a five-tuple, Q = fq 1 , q 2 , . . . , q n g is the finite state set and n = Q j j is the number of states in the automaton. S is a nonempty finite event set with S = S o [ S u , where S o is the observable event set and S u is the unobservable event set of G. d : Q3S ! Q is the state transition function which is considered as a partial function. For an event s 2 S, if the function d(q, s) is defined, we write it as d(q, s)!. q 1 2 Q is the initial state. Q m is the set of marker states with Q m Q.
Note that the set of all finite strings over S is denoted as S * . The empty string e denotes the string with length 0.

Natural projection
Given two event sets: S is a finite alphabet and S o is the set of observable events with S o S, the corresponding natural projection function is P : S * ! S * o , such that P(e) = e and Next, the natural projection P can be extended to the image function, 32 P : P wr (S * ) ! P wr (S * o ) such that Note that P wr (A) is the power set of set A. Mathematically, the inverse function of natural projection P does not exist. Here, the notation of the function P À1 is expressed as the inverse image function, which can be defined as P À1 : P wr (S * o ) ! P wr (S * ) such that Equivalence relation and quasi-congruence relation This section defines three important functions and some concepts associated with the quasi-congruence computation.

Equivalence relation
Suppose that E(Q) is the lattice of equivalence relations on the state set Q. Given an equivalence relation p 2 E(Q), it must meet the following properties. 1,49 (p is symmetric.) 3. (8q, q 0 , q 00 2 Q) qpq 0^q0 pq 00 ) qpq 00 (p is transitive.) In this paper, qpq 0 is also written as q [ q 0 (mod p). For any state q 2 Q, the coset (or equivalence class) of q with respect to the equivalence relation p is denoted as ½q p .
By reflexivity q 2 ½q p , every coset is nonempty. In this paper, Q=p represents the set of all cosets.

Functions
In a DFA G, we can decide the possibly reachable states from any state q i 2 Q based only on observation of the events in S o .
Definition 1. Let s o 2 S * o be an observable string and q i 2 Q a state. The function D is defined as and Note that if d(q i , s) is undefined for every string s with P(s) = s o , then D(q i , s o ) = [. In addition, if s o = s 2 S o , we can simplify the notation of D(q i , s) to D s (q i ), which denotes the set of reachable states from state q i 2 Q via an observable event s. Similarly, if s o = e, it can be written as D e (q i ) which denotes the set of all reachable states from state q i 2 Q via a sequence of unobservable events in S À S o . Moreover, the set of all reachable marker states can be computed by For concise presentation, we define a new event m 6 2 S, and define D m as D m and In addition, we denote n = f1, 2, Á Á Á , ng.
Definition 2. Let p be an equivalence relation, q i 2 Q and s 2 S o . Denote the canonical projection function as P p : Q ! Q=p : q i 7 !½q i p . We can further define its extension projection function as Definition 3. Function P p 8D s is the composition function of D s and P p , At a state q i 2 Q, this composed function P p 8D s can be further expressed as The composite p8D s (s 2 S 0 o ) defines an equivalence relation on Q as follows. Suppose i 6 ¼ j and i, j 2 n.
Quasi-congruence relation Assume that (Q, D) is a nondeterministic dynamic system where Q is the set of system states and D : Q ! P wr (Q) is its state transition function. According to references, 1,31,32 we can define the quasi-congruence relation as follows.
Definition 4. An equivalence relation p on Q is a quasi-congruence for (Q, D), if it satisfies the following three equivalent conditions.
The coarsest quasi-congruence Definition 5. Given a DFA G and the observable event set S o S, the coarsest quasi-congruence 31,50,51 with respect to Q and S o 0 is

Quotient automaton
For a state set Q and the coarsest quasi-congruence relation r 2 E(Q), we have the canonical projection P r : Q ! Q=r as defined in Definition 2. Renaming the cosets in Q=r, we get a new state set Q 0 and a bijection r : Q=r ffi Q 0 . The composition of P r and r is denoted as a function g = r8P r , with equivalence kernel ker g = r.
Definition 6. Given a DFA G = (Q, S, d, q 1 , Q m ), an observable event set S o S, and an equivalence relation r 2 E(Q), the quotient automaton 1,51 is where Q 0 ffi Q=r, q 0 The quotient automaton may be nondeterministic.

Matrix representations
Define the reachable matrix of each event s (s 2 S o 0 ) as a Boolean matrix R s , which is an n3n square matrix and n = jQj. For any i, j 2 n, Let p 2 E(Q) be a quasi-congruence relation on the state set Q. Assume that p partitions Q into m n cosets, then Q=p = fC 1 , C 2 , Á Á Á , C m g and C j represents a coset.
The quasi-congruence relation p is representable by an allocation matrix X n3m , which is also Boolean, that is, To decrease the number of unknown variables, we show that through permutation of the cosets, the allocation matrix is reducible to a lower triangular matrix as formalized in Proposition 1.

Proposition 1.
There must exist a permutation of the cosets fC 1 , C 2 , Á Á Á , C m g such that the allocation matrix X n3m is lower triangular, namely Proof. According to the definition of the allocation matrix X n3m , if m = 1, then there is only one coset C 1 , so the proposition is true. If 1\m n, then we perform the following operations iteratively for each state q i 2 Q in an ascending order of i 2 n. If q i 2 C j and i ø j, then we need not any permutation of columns, that is, x ij = 1, and the rest of the elements in the i-th row are 0. Otherwise, if q i 2 C j and i\j, then x ii = 0 and x ij = 1. In this case, since the previous operations have ensured that (8s 2 i À 1) (8t 2 m) s\t ! x st = 0, the first i À 1 rows of the allocation matrix is already lower triangular. Now that the i-th row violates the lower triangular property, we permute cosets C i and C j and denote the new allocation matrix as X 0 . After the permutation, x 0 ii = 1 and x 0 ij = 0. The i-th row satisfies the lower-triangular property. We need to further prove that the first i À 1 rows of the new allocation matrix still satisfy the lower triangular property.
For all s 2 i À 1 and t 2 m, if t 6 2 fi, jg, then These arguments prove that the first i À 1 rows of X 0 indeed satisfy the lower triangular property.
Set X = X 0 and iterate the operations from the next state q i+1 . The procedure terminates when the allocation of state q m also satisfies the lower-triangular property.
Finally, we can construct the lower triangle allocation matrix X n3m as shown in Table 1.

Quasi-congruence calculation
In this part, we propose a new approach of calculating the supremal quasicongruence using the Boolean allocation matrix and logical propositions.
Let p be an equivalence relation on the state set Q, n = Q j j. Suppose that the number of cosets is m, m n. Since the minimal number of cosets is unknown, we define the Boolean variable y j to represent if coset C j is not empty, where j 2 m.
Then the decision variables are x ij 2 f0, 1g, 1 i n, 1 j m, j i; y k 2 f0, 1g, 1 k m: In total, there are mÁ(2ÁnÀm + 1) 2 + m decision variables. The objective is to find the coarsest quasi-congruence, which has the least number of cosets among all quasi-congruences of the given DFA G based on an observable event set S o . Therefore the objective function for the optimization problem is: To represent quasi-congruences, these decision variables must satisfy the following constraints from C1 to C3.

C1:
For each coset C k , y k = 1 iff at least one state is allocated to it, that is, C2: Each state must belong to one and only one coset, that is, Table 1. Lower triangular allocation matrix X n3m .
C3: For all the state pairs q i , q j 2 Q and all s 2 S 0 o , If two states q i , q j belong to the same coset C k (1 k m), then they must meet the above constraints C3. Equation (20) is further transformed into propositions consisting of binary variables in two steps.
Step 1: The reachable state sets of states q i and q j can be obtained by standard algorithms, assuming that they are represented as Based on the definition of reachable matrix in (13), the two state subsets are equivalent to two row vectors of length n, namely R s (i, : ) and R s ( j, : ), respectively. By Step 2: we need to compute P p 8D s (q i ) and P p 8D s (q j ). The formulas can be simply implemented by Boolean matrix multiplication, denoted by V . The rule is the same as normal matrix multiplication, but replaces multiplication by logic^and addition by logic _. Given a state q i 2 Q, s 2 S 0 o , P p 8D s (q i ) can be represented as Thus, it follows that the equality P p 8D s (q i ) = P p 8D s (q j ) is equivalent to the new vector equality as follows.
Remark 1. R s (i, : ) is a 13n row vector and X is an n3m matrix, then the Boolean multiplication is a 13m row vector. The notation of X (i p , : ) means that all the elements of the row correspond to state q i p in the allocation matrix X n3m .
Therefore, we can further get the logical formulas of this constraint C3: (a) If R s (i, : ) = R s ( j, : ), then the vector equality (21) is always true. So the two states q i and q j must belong to the same coset, that is, ), then we further consider two cases.
One case is that one of the two row vectors R s (i, : ) and R s ( j, : ) is zero but not both. Then q i and q j cannot belong to the same coset. Consequently, X (i, : ) 6 ¼ X ( j, : ). For all k min (i, j, m), When neither of the two reachable state sets is empty, that is, both R s (i, : ) and R s ( j, : ) are not zero, the following condition derived from (20) must be satisfied.
Transforming logical propositions to linear constraints Many constraints in Section 5 are logical constraints. Solving this problem is a Constraint Satisfaction Problem (CSP) or SAT problem. 52,53 A well-known solution to the CSP is to convert the logical constraints into binary linear constraints and to solve the problem by binary linear programming. 54 There are three main logical propositions that need to be converted equivalently. In this section, we will show how to convert them into equivalent linear inequalities.
Proposition 2. If y 2 f0, 1g and x i 2 f0, 1g (i 2 n) are binary variables, then the following two expressions are equivalent.
Proof. ( ) ) For the logic proposition, we need to consider two scenarios: If y = 0, then all x i = 0, i 2 n. Then these values satisfy the linear inequalities, because y n = 0, n Á y = 0, and P n i = 1 x i = 0.
If y = 1, then there exists at least one x i = 1. Therefore, y n = 1 n and n Á y = n, and 1 P n i = 1 x i n. x i n Á y.
( ( ) If the linear inequalities on the right-hand side hold, then we need to prove that the left logical proposition is also true. x i = 0. Under this condition, all x i = 0, i = f1, 2, . . . , ng, thus _ n i = 1 x i = 0. Clearly, it can be seen from the above derivation that the logic proposition is also satisfied in the two cases.
Accordingly, we can conclude that this equivalent transformation is always true. Ã Proposition 3. If x i 2 f0, 1g (i 2 m), and y j 2 f0, 1g (j 2 n) are binary variables, then the following equivalence holds, where N ø max (m, n) is a positive integer. x i ø N . The two linear inequalities at the right-hand side are both true.
If the logical expressions are 0, then all x i and all y j are 0. Then the linear inequalities are also obviously true.
(() If the right-hand side of Proposition 3 holds, then we need to prove that the left logical formula is also true. There are two cases as below.
1. If all y j = 0, then , all x i = 0. Obviously, in this case, the logical formula _ m i = 1 x i = _ n j = 1 y j is satisfied.
(2) If there exists at least one y j = 1, then Since x i 2 f0, 1g (i 2 m) and N ø max (m, n) is a positive integer. Thus, P m i = 1 x i ø 1 N .0, i.e., there must exist at least one x i = 1. Clearly, these values also satisfy the logical formula _ m i = 1 x i = _ n j = 1 y j . Therefore, it can be seen from the above derivation that the Proposition 3 is always true. Ã Next, we introduce the following lemma for establishing Proposition 4.
Proof. ( ) ) According to Proposition 3, we can first equivalently convert the logical implication into the statement as below.
Second, based on Lemma 1, we can translate the logical implication (27) into the following form.
The next step is to show that if the two propositions in (28) hold, then the two inequalities of (26) also hold. The proof has two cases.
Because N is a positive integer, À 1 N \0. Further, we can transform the above inequality set (29) into the following two inequality pairs.
Summing the inequality pairs in (30) and (31), we have the following two inequalities.
Since m, n.0 and N ø max (m, n), we have À1 À m N and ÀN Àn. From this we derive the following pairs.
Recalling z 1 + z 2 À 2 À1, then both inequalities in (26) are also satisfied, regardless of the values of x i and y j . ( ( ) We prove this direction in two cases.
1. If z 1 = 0 or z 2 = 0, the proposition z 1 = 1^z 2 = 1 is always false. The logical proposition at the left-hand side is always true regardless of the values of x i and y j in this case. 2. If z 1 = 1^z 2 = 1, then z 1 + z 2 À 2 = 0. In this case, we can convert the set of inequalities in Proposition 4 into the set of inequalities below.
The inequalities (35) are equivalent to the right-hand side inequalities by Proposition 3. Therefore, we can know that the equivalent conversion of Proposition 4 is always true. Ã

Main algorithm
According to Proposition 2, it is easy to convert the constraints of C1 into the following pairs of linear inequalities. For any k 2 m, Moreover, the constraints of C2 are already linear equality constraints. Finally, all linear inequality constraints of C3 are obtained by Algorithm 1. In the following Algorithm 1, n is the number of states, m is the maximal number of cosets, (R s ) n3n represents the reachable matrices for each event s 2 S 0 o , and X n3m (m n) is the allocation matrix we supposed.
This algorithm computes the linear inequalities for all pairs of states q i and q j according to (20). Line 4 determines whether the two states cannot be equivalent. In that case, the algorithm adds the inequalities specified by lines 5-7, skips the remaining part, and iterates to the next value of j. Otherwise, the algorithm continues into the loop of lines 11-24. Lines 12-13 check whether the reachable state sets from both states q i and q j are equivalent. If they are, no constraint is added; otherwise, the algorithm continues to execute lines 15-22 and output the linear inequalities of lines 17-22. In the algorithm, line 6 corresponds to formula (22) in Section 5. Line 12 represents case (a) of the constraints C3. Lines 14-23 implement the second case in case (b) of C3. Finally, lines 17-22 correspond to the linear inequalities from the transformation of the logical formula (23) according to Proposition 4.

Model abstraction
This section presents the complete procedure of computing the model abstraction.
Input: A finite automaton G = (Q, S, d, q 1 , Q m ) and an observable event set S o S.
Step 1: Compute the reachable state set function D s for all s 2 S 0 o . According to the definition of reachable matrix in (13), the reachable matrix of each event s can be obtained as (R s ) n3n . for t = 1 to m do 19 Step 2: Compute the coarsest quasi-congruence by solving the optimization problem (17)- (20) using the proposed method in the paper, based on all the reachable matrices (R s ) n3n and s 2 S 0 o . Step 3: Compute the quotient automaton for G with the coarsest quasi-congruence, according to Definition 6.
Output: The quotient automaton is the abstraction of G.

An illustrative example
This section elaborates the complete calculation process of the method with the simple example in Figure 1. The automaton M2 has the state set Q = fq 1 , q 2 , . . . , q 6 g, the event set S = fa, b, gg, the initial state q 1 and the marker state set fq 6 g. Given an observable event set S o = fa, gg, we calculate the coarsest quasi-congruence through the following steps. First, we need to calculate its reachable matrices for all observable events, including (R a ) 636 , (R g ) 636 , and (R m ) 636 as in Tables 2 to 4. Second, let m = 6 and the corresponding lower triangular allocation matrix X 636 is in Table 5. Third, on the basis of Propositions 2-4 and the constraints C1-C3 in Section 5, we can list all linear inequalities constraints for M2 as follows. Here, n = 6 and m = 6. C1: y k = _ n i = k x ik , (k 2 m) , 1 n Á y k À x kk À x (k + 1)k À Á Á Á À x nk 0 x kk + x (k + 1)k + Á Á Á + x nk À n Á y k 0 x ij = 1, (i 2 n) , : For all pairs of i 2 n À 1 and i + 1 j n, for all s 2 S 0 o , we need to consider the following three cases in turn. Firstly, if one of the two row vectors R s (i, : ) and R s ( j, : ) is zero but not both, then states q i and q j cannot belong to the same coset. Then we add the following inequalities.
x 21 + x 41 1 Secondly, if the condition R s (i, : ) = R s ( j, : ) is true for all reachable matrices (s 2 S 0 o ), then the two states q i and q j must belong to a coset, and no constraint (R a ) 636 q 1 q 2 q 3 q 4 q 5 q 6 Table 3. Reachable matrix (R g ) 636 of M2.
needs to be added. From Tables 2 to 4, we can see that ½R a (1, : ) = R a (3, : )^½R g (1, : ) = R g (3, : )^½R m (1, : ) = R m (3, : ), so we can determine that q 1 and q 3 must belong to a coset, and there are no constraints to be added in this case. The final case is that R s (i, : ) 6 ¼ R s ( j, : ) and both are not zero vectors. Since M2 is simple, the third case does not exist in this example. If there is such a case in more complex automata, we also need to add the following pairs of linear inequalities according to (23).
Note that the symbol Z s (i, t) in the formula (38) represents the new vector multiplication, which can be calculated according to formula (21) and expressed as follows.
Here, s 2 S 0 o , i 2 n and t 2 m. On this basis, we can compute the corresponding product matrix (Z s ) n3m . For example,  Table 5. Supposed matrix X 636 of M2.
x 61 x 62 x 63 x 64 Thus, if k 2 m, then Z a (2, : ) = R a (2, : )^X ( : , k) = Therefore, this example is formalized as the following binary linear programming problem. The objective function for the optimization problem is: subject to the following constraints.
x kk + x (k + 1)k + Á Á Á + x nk À n Á y k 0, k 2 m ð43Þ x 2k + x jk 1, j 2 f1, 3, 4, 5, 6g, 1 k min (2, j) ð45Þ x 4k + x jk 1, j 2 f1, 2, 3, 5, 6g, 1 k min (4, j) ð46Þ x ij 2 f0, 1g, i j 2 n ð47Þ In total, this simple example has 27 variables, 6 linear equalities and 35 linear inequalities. The number of linear inequalities obtained by C3 is 23. The problem is solved by MATLAB with a PC with Intel(R) i7-4600U CPU @2.10 GHz, 2.70 GHz, and 16.0 GB installed memory (RAM). The computation time is 0.0135 s. Table 6 shows the final result of the coarsest quasi-congruence relation for M2, and its state partition is p = ffq 1 , q 3 g, fq 2 g, fq 4 g, fq 5 , q 6 gg, as shown in Figure 2. All the states in a dotted box represent a coset. Figure 3 further illustrates the simplified quotient automaton of M2. The abstracted model is a nondeterministic finite automaton, because there are unobservable transitions from states 0 and 2 to state 3. This suggests that the selected observable event set is not proper for simplifying the original DFA as a smaller DFA. In this case, we may need to consider re-selecting an observable event set. 32 Finally, the coarsest quasi-congruence is also computed by the standard method in TCT and the result confirms the correctness of the new BLP approach.

Comparisons of computation methods
This section compares the computation results of three different methods: the classical graph-partition method 1,45 implemented in TCT, and the binary linear programming method proposed in this paper which are respectively implemented by mixed integer linear programming (MILP) 55 in MATLAB and binary linear programming in CPLEX (BLP). Table 6. Allocation matrix X 636 of M2.   Table 7 is a comparison of computing the supremal quasi-congruence relations on the basis of different observable event sets of the automaton M2 by using the stated three methods. Column 1 lists different sets of observable events. Column 2 lists the result using the graph-based partition method implemented in TCT 45 for each S o . Column 3 lists the assumed m value, where each m is equal to 6 for comparison of calculation results. Column 4 shows the size of the linear inequalities C3 obtained by Algorithm 1. Since the size of constraints C3 (Number of linear inequalities 3 Number of variables) is the main factor that affects the computation time of the method proposed in the paper, it is necessary to analyze the relationship between its size and running time. Columns 5-6 show the running time using two different methods, the MILP method in MATLAB and the BLP method implemented by calling CPLEX Studio 12.8 in MATLAB. Note that when the running time of BLP method is less than 0.01 s, the system automatically reduces it to 0.00 s. The computation results are summarized as follows.
The computational results using MILP and BLP methods are the same as that of using the computation software TCT. On different sets of observable events, when m is invariant, the scale of linear inequalities C3 is not necessarily the same. In general, however, the smaller the constraint size is, the shorter the run time will be. Moreover, the smaller the number of observable events is, the smaller the number of cosets will be. Table 8   For the same automaton and the same observable event set, a tight estimate of m is important. A smaller m generally results in less computation time. On the other hand, if m is less than the real minimum, the binary linear programming problem is infeasible. As the state space increases, the calculation time and memory consumption of MILP method also increase, mainly due to the increase of the number of linear inequality constraints C3. When the automaton is not large, the MILP method is efficient. The BLP method is more efficient than MILP method under the same conditions.

An application
The proposed model abstraction method is applied to simplify the supremal nonblocking supervisor of a manufacturing system with five automated guided vehicles (AGVs), 1,51 as shown in Figure 4. The simplified supervisor model exposes the control logic hidden in the original large automaton model. The manufacturing system has two input stations IPS1 and IPS2 for parts of types 1, 2; three workstations WS1, WS2, and WS3, each containing a conveyor belt; one completed parts station CPS; and five AGVs. The five AGVs alternately load and unload parts on their own fixed circular routes.
The system model consists of five automaton models of the five AGVs, eight automaton models of control specifications. Figure 5 shows the automaton models of the five AGVs. The states corresponding to the physical zones are illustrated as gray rectangles. The solid and dotted lines represent the loading and unloading of each AGV respectively. Each of the four shared zones is modeled as an automaton (Zone j , j = 1, Á Á Á , 4) that restricts only one AGV in the zone. Each of the three workstations is modeled as an automaton (WS k , k = 1, Á Á Á , 3). WS1 receives type 1 parts via AGV3 and type 2 parts via AGV4, and assembles both into the complete products, which are transferred to CPS via AGV5. WS2 receives type 1 parts via AGV1 and processes them. The finished type 1 parts are transferred to WS1 via AGV3. WS3 receives type 2 parts via AGV2 and processes them. The finished type 2 parts are transferred to WS1 via AGV4. The dotted zone of IPS is modeled as an automaton, which restricts only one AGV inside the zone.
More details of the AGVs system can be found in Wonham and Cai, 1 Feng and Wonham. 51 The plant model is the synchronous product of the five AGV models. Note that the two values in parentheses respectively represent the number of states and transitions in the automaton model obtained by each calculation.
On this basis, the supremal nonblocking supevisor of the AGVs system is ZWISUP :¼ Supcon(AGV, ZWISUP) (4406, 11338): Namely, the full AGVs system consists of 4406 states and 11338 transitions. The automaton model of the supervisor is too complex to be comprehended. We apply the proposed model abstraction method to simplify the supervisor model so that the control logic becomes clearer. Deadlock arises in many manufacturing systems when the number of being-processed parts in the system is more than a limit. We investigate this limit through the model abstraction method. Type 1 parts enter the system via event f11g and type 2 parts via event f21g. The completed products exit the system via event f51g. So we choose the observable event set as f11, 21, 51g, then compute the supremal quasi-congruence of the supervisor model ZWISUP using the method proposed in the paper. Finally, we obtain a simplified supervisor model for AGVs system with 22 cosets and 45 transitions. Figure 7 shows the quotient automaton of ZWISUP with respect to an observable event set f11, 21, 51g. The simplified model reveals a constraint between events Figure 6. Eight specification automaton models for the AGVs system. rence of event 12 (11). The reason is that after three successive occurrences of event 11 (12), there are three type 1 (2) parts in the manufacturing system. One is in WS1 waiting for the assembly with the other type of part, one is in AGV3 (AGV4) waiting to be transferred to WS1, and the last one is in WS2 (WS3). If event 11 (12) would occur once more before the other event occurs, AGV1 (AGV2) would have to park within the shared zone IPS, because WS2 (WS3) is full. Then AGV2 (AGV1) cannot enter the input station IPS2 (IPS1). The final product can never be assembled at WS1 and deadlock thus arises.

Conclusions
In this paper, an optimization process is designed to automatically convert the computation of the coarsest quasi-congruence into a binary linear problem. Partitioning states to cosets is formalized as a binary allocation matrix X n3m , where n is the number of states, m is an estimated number of cosets. Evidently m n. The task of finding the coarsest quasi-congruence is equivalent to the objective of finding the least number of cosets. To reduce the computation time, the translation process is also optimized by introducing fewer decision variables and simplifying constraints in the formulas. The computational efficiency and correctness of this method are verified by two different solvers. Therefore, the proposed method is not only a new solution for the coarsest quasi-congruence computation, but also provides us a more intuitive understanding of the quasi-congruence relation in the supervisory control theory. A future direction is to apply more computationally efficient solvers, such as satisfiability modulo theories (SMT) and SAT, to compute the optimal state-tocoset allocation problem. Another direction is to apply the same method to compute model abstraction of finite automata with variables, namely the extended finite automata which are easily applicable to the real industrial scale problems. Systems, Man, and Cybernetics Society and served as a member of IFAC Technical Committee on Discrete-Event and Hybrid Systems, from 2011 to 2014. He was a recipient of the Alexander von Humboldt Research Grant, Alexander von Humboldt Foundation, Germany. He is listed in Marquis Who's Who in the World (27th ed., 2010). He is the Founding Chair of Xi'an Chapter of IEEE Systems, Man, and Cybernetics Society. He serves as a frequent reviewer for more than 90 international journals, including Automatica, a number of IEEE Transactions, and many international conferences.