Digital-Twin Consistency Checking Based on Observed Timed Events With Unobservable Transitions in Smart Manufacturing

Smart factories manage digital twins (DTs) to evaluate the performance of various what-if production scenarios. This article presents a DT consistency-checking approach to maintain DT in high fidelity by checking whether each sensed timed event from the physical manufacturing plant is under its corresponding DT-based estimations in runtime. The approach targets DTs developed using time colored Petri net (TCPN). To build the candidates of the next observable event with observable time margins, we considered the stochastic property of the plant, frequent external actuation caused by a new order, machine maintenance, etc., as well as intermediate unobservable state transitions reaching the sensible events. Based on the considerations, we propose an iterative method to build the virtual estimates for streaming physical events using efficiently evolved state-class graphs (SCGs). We also propose a TCPN partitioning method to accelerate the SCG-evolution and make DT maintenance easier by supporting the isolation of inconsistent subnets being diagnosed. We applied the approach to a USB flash-drive factory to prove the concept and evaluated the performance under various situations to show speedups of the SCG evolution, that is the crucial overhead of the estimation.

performance, such as throughput and cycle time, after virtually commissioning various what-if scenarios. For reliable performance estimation, it is crucial to maintain high-fidelity DT that reflects the dynamics of part (or material) flows in the high-mix/low-volume manufacturing lines. For that, DTs should be developed considering various production behaviors of the target plant's components, such as machines, conveyors, human operators, and automated guided vehicles (AGVs). Due to a modeling error, an unreported physical change, or an unmodeled exception, a DT-based prediction might become inaccurate. This article focuses on a DT consistency check that verifies DT using physical events sensed from the plant along with their observed times (called timed events). Sufficient fidelity in DT can be guaranteed by iterative synchronization of DT with its physical counterpart throughout the manufacturing life cycle.
To reflect the production dynamics, DT is modeled using discrete-event models (DEMs). DEMs can describe the dynamics of the plant's operations using state transitions and exchanges of the events generated by the transitions. Examples of the events include part loading/unloading, AGV arrival/departure, machine process start/end, setup or operation mode changes, fault detection, and fault recovery start/end. Each physical timed event has corresponding virtual estimates, which comprises multiple event types and their observable time margin, as Fig. 1 shows. The consistency between the plant and its DT is examined iteratively by checking whether each event is consistent with its virtual estimates. If all state transitions and related events are observable and the plant's production is deterministic, building the estimates might be straightforward. However, computing the virtual estimates presents the following three main challenges.
This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ 1) Some state transitions that do not produce physical events (which can be monitored online) are unobservable due to the lack of monitoring or communication tools for legacy machines or human operators. 2) There are stochastic state transitions and state-duration variances due to faults in the machines, human operators, or AGVs. For example, a machine's arm can locate a part in the correct position after several attempts, or an obstacle can suspend a moving AGV for a while. 3) External interventions from the control facilities (CFs) to the plant can occur to handle new orders or manage machines for energy efficiency. We aim to compute the virtual estimates of the plant's next physical event based on reachable observable transitions, considering all possible sequences of intermediate unobservable transitions, stochastic state transitions and variable operation periods, and external interventions. For the reachability analysis, researchers have identified reachable transitions (or states) based on the states' time ranges and all possible transitions (including unobservable transitions) using time Petri net (TPN) [1], [2], [3], [4], [5]. In the studies, various state-class graphs (SCGs) and their evolution methods for the analysis have been developed. In this article, we propose a DT-modeling formalism of a time-colored Petri net (TCPN) to apply the TPN-based analysis methods to high-mix manufacturing plants. It is designed by extending TPN to use colors that distinguish tokens representing various types of parts. By extending a modern TPN-based SCG evolution method, we propose a runtime checking method to update SCG and compute the virtual estimates of the plant's next event based on the recent physical event (for further consistency checking), as Fig. 1 shows. We denote the CF's external interventions as actuation events (AEs) and the plant's sensed events as reactive events (REs). DT's observable transitions can be detected online using associated REs and AEs. Assuming AEs are unpredictable, only REs are used to check DT consistency. AEs are only used to revise SCG to update the virtual estimates.
We summarize the main contributions as follows. 1) We propose an approach to check the consistency of TCPN-based DTs during runtime using streaming timed physical events (monitored online). 2) We extend TPN-based SCG-evolution methods to estimate the next RE, considering AE's unpredictability. 3) We propose a TCPN partitioning method based on the AE's property (unpredictability), which improves the overall SCG evolution computational performance and supports the isolation of problematic subnets (which can examine other consistent subnets during the problem diagnosis). 4) We propose methods to update the SCG iteratively with newly observed event to maintain DT's consistency. The rest of the article is organized as follows. Section II introduces the related works about TCPN and TPN-based reachability analysis. Section III defines the formalism of TCPN. Section IV specifies the problems related to SCG in the runtime checking. Section V details the overall methods for constructing the SCG and estimating the next REs. Section VI shows the TCPN partitioning method and SCG revision according to the observed AE/RE. Section VII presents a group of experiments as a case study. Finally, Section VIII concludes this article.

II. RELATED WORKS
Smart manufacturing systems can be generally structured in three layers: 1) enterprise resource planning (ERP); 2) manufacturing-execution system (MES); 3) industrial control system (ICS) [6]. At the top level, ERP focuses on integrating organizational functions to provide forecasting and planning, inventory management, and accounting functionalities. In the middle level, MESs have been used to collect and manage the information from manufacturing plants and helps decision makers understand the current status and improve productivity [7]. The decision making requires the evaluations of various what-if manufacturing scenarios using the plant's high-fidelity DT. Our work is meant to increase DT's fidelity by checking the consistency between DT and monitored physical events. At the lowest level, ICSs are used to control and monitor specific equipment.
In previous works, researchers have examined the consistency between DT and physical information to detect abnormal behaviors in the ICS's security domain [8], [9], [10]. They considered DT as a reference model and noticed the adversary's attacks by detecting the inconsistency between physical and virtual states. Compared to the works from the security domain with our method, we aim to check DT's consistency considering unobservable transitions, stochastic state transitions, stochastic operation periods, and users' interventions (which are AEs).
Petri net (PN) has been widely used to describe plants' eventdriven behaviors for various purposes, including throughput measurement [11], [12], supervisory control [13], [14], and deadlock-free scheduling [15], [16], because of its advantages in modeling concurrent and synchronization operations. PN has been extended to include various features and functions for specific modeling purposes, such as colored PN (CPN) [17], [18], queueing PN [19], and time PN (TPN) [20], [21] or timed PN [22]. CPN adds an attribute to tokens as color to distinguish between multiple production sequences for different products in a complex plant. Some CPNs can also support the primitives to model specific data manipulation. TPNs or timed PNs help analyze the performance of timed systems.
The main distinction between TPNs and timed PNs is whether each enabled transition fires within a given interval or a given delay value, respectively. We represent the plant's stochastic state durations as intervals based on measured means and standard variations (see Section III). Then, we propose a new modeling formalism of TCPN to: 1) combine TPN's interval-based transition firing; 2) employ CPN's colors (for high-mix manufacturing) and AEs (for users' interventions); 3) reference previous TPN-based reachability-analysis techniques [to build the virtual estimates by identifying all reachable states (or transitions) that generate the REs]. The TCPN's formalism was not previously defined to the best of our knowledge.
Berthomieu and Diaz [1] were the first to propose an abstracted TPN state, called state classes (SCs), to represent the finite or infinite sets of the TPN's timed states. They presented a graph called a SCG-based on the SCs. The SCG's nodes are the set of SCs, and each edge shows the consequence of a labeled transition execution (or firing) from one source SC to another. Some studies have focused on better abstractions (reducing the number of SC nodes) by proposing a relaxed SCG, contracted SCG, or fault-diagnosis graph (which adds the sequence of transition information to the edges and removes the unnecessary SCs) [2], [3]. In [4], Basile et al. presented a modified SCG (MSCG) that extends edge labels using additional transition information related to timing variables and constraints. Using the MSCG, they proposed a fault-diagnosis method that requires solving a linear programming problem to detect the occurrence of a fault-related transition for a given sequence of observed events. In [5], He et al. showed that the previous MSCG evolution might miss some reachable states. Therefore, they consolidated the MSCG evolution method by defining deficient SCs that can generate new paths (series of SCs) using new subsequent SCs.
SCGs show all reachable SCs caused by possible observable and unobservable transition firings with their range-based time intervals. Researchers have typically utilized SCGs for fault diagnoses that check whether any fault transition might occur based on a given sequence of observed events. The proposed approach shares the existing TPN-based methods' principles in terms of SCG evolution but requires a revision based on TCPN. Moreover, when we use existing SCG-related methods for runtime DT checking, we need additional procedures, which are: 1) building the virtual estimates based on the SCG; 2) updating SCGs iteratively for streaming physical events. An SCG should be revised for the subsequent event estimation in runtime checking whenever a new physical event arrives. For the SCG update, we compute consistent SCs (which become roots for the next SCG) based on the previous SCG and each physical timed event. Among various types of existing SCGs and their evolution methods, our work uses MSCG (in [4] and [5]) because MSCG's timing information annotated with its edges is essential to derive consistent SCs.
For runtime checking, speedup of SCG-evolution (which is the main estimation overhead) is essential because early inconsistency guarantees better problem diagnosis by dispatching human operators to a problematic physical location. Thus, it is recommended to build the virtual estimates using the updated SCG before a new physical event arrives. For the evolution speedup, we resolve computational inefficiency problems of the latest MSCG evolution method in [5] (discussed in Section IV-A) and propose a TCPN partitioning method.
The following section defines a TCPN's formalism for the formal description of our approach and the plant's TCPN-based DT modeling.

III. TIME COLORED PETRI NET
TCPN is a combination of TPN and CPN for: 1) TPN's transition-associated time intervals; 2) CPN's token distinction. In our work, the tokens represent high-mix parts or other part-flow control signals (e.g., machine availability, delivery request/grants, etc.). Transitions represent part-specific or joint operations (with their range-based time duration) for machine operations, flow-control decisions, AGV movements, etc. We simplify a transition's logically enabling condition of TCPN as the required number of colors in the transition's backward places, similar to TPN's enabling condition (the required number of tokens). For that, we represent a token, tk, as a pair of colors (for the transition's token distinction) and data (which can be empty), i.e., tk ∈ C × (D ∪ ∅), where C is a set of colors and D is a set of possible manufacturing-related data. Then, it requires an additional color decision when the tokens are fired for the subsequent transition (if needed).
Based on the concept, we define a TCPN's formalism as follows: Definition 1: The net of TCPN is represented as N = P, T, En, Pre, Post, Q, L , where 1) P is a set of places.
2) T = T re ∪ T ae ∪ T si is a set of transitions, where T re /T ae is a set of observable transitions using their associated REs/AEs and T si is a set of silent transitions. 3) En : P × T → (C × N 0 ) |C| is an enabling-condition function to specify the required numbers of colors for a given transition t ∈ T from a place p ∈ P to be logically enabled. 4) Pre : {M } × T → {X * } is a function to show one or multiple candidates of inflowing tokens when a transition t ∈ T is fired at the given marking M ; a marking M : P → {K} is a vector that assigns each place a multiset of tokens; K is a multiset of tokens, and X : P → {K} is a vector of tokens inflowing to t from each place. 5) Post : {X} × T → {Y } is a function to compute a vector of outflowing tokens Y : P → {K} based on a given vector X when transition t ∈ T fires. 6) Q is a timing function that defines the set of static closed time intervals only for predictable transitions (T re ∪ T si ) as Q : T ae → ∅ and Q : assigns an RE's or AE's type (using a letter, λ) to each observable transition, t ∈ T re ∪ T ae , or the empty string to each silent transition t ∈ T si . LetM be a conversion from marking M after projecting the tokens onto the color domain, as Fig. 2 shows. We de-noteM (p) as the numbers of token colors in place p. A transition t(∈ T ) is logically enabled if sufficient colors are present inM , i.e.,M ≥ En(·, t), where En(·, t) is the vector of (En(p 1 , t), . . . , En(p m , t)), p i is the ith place, and m = |P |. The inequality ≥ denotes the relationship: ∀p i ∈ P,M (p i ) ≥ En(p i , t). We denote the set of transitions logically enabled at M as A(M ), i.e., A(M ) = {t ∈ T |M ≥ En(·, t)}. This color-based enabling description can make the same-colored tokens with different data move through a shared transition, which helps represent a common process of high-mix parts.
The function Pre is developed considering the corresponding actual system. Let us suppose that a machine randomly selects two parts (which require a t 1 -related operation) from its buffer, as Fig. 2 shows. Then, Pre is developed for Pre(M, If the machine selects two parts by the arrival order and token tk 2 arrived earlier than tk 1 , then tokens should contain the arrival-order information in their data and The function Post can forward each token tk in a given X to the next place after the token's color or data conversion. Post can also drop or generate tokens on purpose. When a token vector, X ∈ Pre(M, t), is moved by the firing of a transition t, M is changed to the next marking, where and ⊕ are the element-wise set operations of \ and ∪ on two vectors, respectively.
The function Q does not return any interval for AE-related transitions because of its unpredictable occurrence. Therefore, an AE-related transition, t ∈ T ae , is fired only by its real AE observation. For the reactive transitions (which are T re ∪ T si ), Q indicates the intervals of those transitions using two rational numbers, namely Q(t) = (l, u), where l ≥ 0 and u ≥ l. In our work, a logically enabled transition, t ∈ T \ T ae , must be fired within the interval as a constraint unless it becomes disabled. The constraint is called strong time semantic (STS) in [4]. The time interval of Q is decided based on its corresponding operation spans, which are reported online or measured manually. The interval of a transition t ∈ T \ T ae is represented as , whereŝ i and σ i are the means and standard variation of measured spans of a t-related operation, respectively, and k is a user's empirical constant, such as 3 (for the three-sigma rule of thumb).
We assume that each observable transition in T re ∪ T ae does not share its label with any of the other transitions. This assumption makes us find a specific transition corresponding to an event type λ, i.e., t = L -1 (λ). Suppose that a machine processes two steps: a common step (denoted by the transition t 1 ) and a variable token-dependent step (among two tasks symbolized by t 2 and t 3 ) and we can only detect the completion of the second step using an infrared sensor, as

IV. PROBLEM STATEMENT IN SCG-BASED DT CHECKING
Our approach is intended to check DT's consistency during runtime via three methods: 1) evolving the SCG (to build the next RE's virtual estimates) based on TCPN; 2) splitting the TCPN (for fast SCG evolution and subnet isolation); 3) reupdating the SCG based on each physical event. In this section, we describe each methods' SCG-related problems.

A. Background and Problems of SCG Evolution
As discussed in Section II, we uses the MSCG in [4] and [5] as a default SCG. First, we introduce the background of MSCG and its evolution. Then, we describe the latest evolution algorithm's computational inefficiencies. Last, we clarify additional requirements for runtime TCPN checking.
The SCG is a directed graph whose nodes are called SC. Each SC is associated with the net's state and represented as a pair of a reachable marking (M ) and a set of inequalities (Θ) that define the timing constraints of logically enabled transitions at M , i.e., C k = (M k , Θ k ), where C k is the kth SC. Each inequality in Θ means the remaining time before its associated transition's (t's) firing if t does not become logically disabled due to a token preemption by another transition firing. The inequality can depend on a certain number of variables, denoted Δ variables, which consider how long a transition has been enabled. Let θ i ∈ Θ be transition t i 's timing constraints. Then, we generally represent θ i as is the time spent of the th previous SC C (-) ; and d i is the number of intermediate SCs after t i was newly enabled. In the case of SC C 5 (which is one of C 2 's next SCs) in Fig. 4, Δ (-) and d 3 of θ 3 (∈ Θ 5 ) are Δ 2 (that is, C 2 's time spent) and 1, respectively. This indicates that transition t 3 was newly logically enabled at C 2 . No timing variable in θ 5 ∈ Θ 5 , i.e., d 5 = 0 means that transition t 5 is newly logically enabled at C 5 .
Each edge e of SCG has an extensible label where t i is the transition whose firing leads to the target SC's marking; L(t i ) is t i 's label; and Δ k is a constraint of the source SC's time spent. The constraint bounds l * k and u * k are functions of Δ variables, as Fig. 4 shows. As in [5], the label can include a function, f eq : Δ n → Δ k , if an existing SC can be reused as the edge's next after one or multiple existing time-spent variables {Δ k } (depicted in the edge's source) are renamed with their pairs {Δ n } of the reusable SC. In Fig. 4, C 5 's forwarding edge contains f eq to reuse C 3 after renaming Δ 5 as Δ 3 . The reuse can reduce the evolution overhead because the reused SC will not be further explored if it is not a deficient SC. Definition 1 in [5] defines the deficient SC, which could generate a new path from the deficient SC when a new deficient SC reuse (for a new path reaching the deficient SC) lessens a firing constraint of a logically enabled but previously preempted transition (called deficient transition) enough to be fireable. Fig. 4 depicts the deficient SC example, which is C 5 , because it has the deficient transition θ 3 preempted by the prior transition t 5 based on max(θ 5 ) < min(θ 3 ). However, if C 5 is reused and its new backward edge contains a time spent, such as Δ 2 ∈ [1, 2], then t 3 can be fired first as an imminent transition, which will generate a new subsequent SC.
Based on the concept of SCs' deficiency and reuse, the previous SCG evolution algorithms in [5] have the following two computational inefficiencies: 1) Whenever an SC (from a queue) is called, imminent transitions (which can be fired before the others) are computed based on all the previous paths reaching the SC (as Algorithm 2 in [5]). 2) All deficient SCs are called whenever the queue for the waiting SCs (to be evolved) is empty to find possible new paths (as Algorithm 1 in [5]). Because deficient SCs can be called multiple times, computing the imminent transitions based on all previous paths can lead to serious computational overhead. The iterative checking of all deficient SCs is inefficient because most deficient SCs would not lead to new paths as the SCG evolution cycle elapsed. To prevent each SC's all-path checking, the proposed method schedules each SC with its new paths and timing information, defined as follows: Definition 2: Given an SC C k in the queue, C k has its own new paths Π + k and related previous SCs' time- Π} is a set of new paths reaching C k from the roots; Π is the set of full paths (that comprise all SC sequences from each root to terminal SCs that do not have any next reachable SC) based on SCG; π a is the ath path in Π; C (1) is the root on π a ; and C (n Based on Φ + k , we define the path-dependent newly observable time range of C k as follows.
Definition 3: Given a path π [1:n] a = C (1) . . . C k ∈ Φ + , C k is newly observable in a range of I acc (π (2) The range of C k 's newly observable time based on a given path π [1:n] a is the same as the time range of C k 's last transition t (n-1) , where (t (n-1) , . . .) = L(e (n-1) ); e (n-1) is the edge between C (n-1) and C k .
To prevent the iterative invocation of deficient SCs, we schedule particular deficient SCs when they are engaged in a new path as a necessary condition. In addition, the proposed evolution method considers: 1) TCPN (in Definition 1); 2) different stopping criteria sufficient to find all candidates of the next observable transitions (related to REs), without reaching all possible states.

B. Net Partitioning and Event Estimation Problems
First, we introduce why TCPN partitioning is required for fast SCG evolution and the requirements of TCPN partitioning. Then, we present preliminary notations for virtual estimates.
The number of SCs in an SCG increases exponentially with TPN system complexity (related to the net structure and tokens in the initial marking), as in [4]. Let C 0 be original TPN's (N 0 's) system complexity, and let size(C 0 ) be the expected number of SCs induced by C 0 . If the net and marking of N 0 can be split into two nets, N 1 and N 2 , with distributed markings as |C 0 | = |C 1 | + |C 2 |, then the overall SC nodes are greatly reduced, for example, size(C 1 ) + size(C 2 ) size(C 0 ) based on the exponential node growth by the system complexity. The net splitting for the complexity distribution should support the independent evolutions of multiple subnets' SCGs that lead to the same virtual estimates (which are computed using the orignal net's SCG). The detailed partitioning method that satisfies this requirement will be described in Section VI-A).
After partitioning, there can be multiple subnets {N m }, where N m = (P m , T m , . . .). Each subnet N m has its own SCG G m , which leads its own virtual estimates. Each subnet's virtual estimates consist of a set of observable RE's types Λ m and their observable time intervals T m . To compute the estimates, each subnet N m manages its own full paths Π m of SCG G m (see Definition 2), path-related time constraints Φ m = {I(π a ) | π a ∈ Π m }, and last event's observed time τm . The estimation methods will be described in Section V-B.

C. SCG Update Problem
If a physical event is observed, its dedicated SCG handles the physical event. The proposed net-splitting method enables subnets to share observable transitions whose backward and forward places are located in different subnets. If the physical event's transition is shared, the number of related SCGs can be two. After checking DT's consistency using the event (if the event is an RE), a related SCG G m is updated by three processes: 1) finding the consistent SCs; 2) revising the SCs to meet a root's property (defined below); 3) starting the SCG evolution based on the roots. In [4] and [5], an initial marking (as a root) is given. In our approach, an initial marking should be given only if a previous SCG does not exist. However, more than one root is iteratively derived during runtime using the previous SCG and each physical event.
An RE-consistent SC has a backward edge whose transition is related to the RE and can occur at the RE-observed time. An AE-consistent SC is an SC that has a marking where an AE-related transition is logically enabled and can stay at the AE-observed time. We will formally define the consistent SCs in Sections VI-B and VI-C. As in (1), each SC's timing constraint θ i can depend on previous SCs' time constraints ({Δ (-) }, where > 0). Thus, we should examine each SC's timing consistency considering all paths Π m using their time-spent information Φ m .
After finding the consistent SCs, we should revise the SCs so that each SC's constraint θ i does not include previous SCs' constraints {Δ (-) } because there are no more previous nodes of the SC. For that, we define roots' property as follows.
Definition 4: Given a root C k = (M k , Θ k ), each constraint θ i ∈ Θ k is a pair of two constants, i.e., d i = 0.
When the physical event is related to a shared transition t i between two subnets, a subnet that has t i 's forward places (denoted by t • i ) could need the inflowing token information from the other subnet's SCG to compute t i -induced next marking. To specify the exchanging token information between subnets, we extend the label of SCG's each edge e to contain inflowing tokens X j , such as L(e) = t i , . . . , X j .

V. DT-BASED EVENT ESTIMATION
This section presents the methods for the proposed DT consistency checking. First, we discuss the evolution of SCG using TCPN in Section V-A. Next, the estimation of next REs are discussed in Section V-B.
The overall architecture is illustrated in Fig. 5. Before checking, the modelers should develop the TCPN N based on plantspecific knowledge. Then, they should collect operation timing data related to the transitions manually or automatically to build the TCPN's time function Q. Depending on the transition's observability, N can be split into multiple subnets {N m }, where N m = (P m , T m , . . .) is the mth subnet. The subnets only share some observable transitions, and we will discuss the detailed methods of TCPN partitioning in Section VI-A. Each subnet represents the production dynamics in a specific location of the plant. The subnet N m is used to compute its own SCG G m based on initial consistent SCs to estimate the subnet's next RE. The event monitor (EM) converts the streaming physical events and the external interventions into related REs and AEs, respectively.
When an λ-type RE is observed at time instant τ , the EM's specific procedures are as follows.
1) The EM selects an SCG G m , whose net N m includes the λ-related transition t i , that is ∃t i ∈ T m s.t. t i = L -1 (λ). If t i is shared by two subnets, then N m satisfies • t i ∈ P m , where • t i is the set of the backward places of t i . 2) The EM checks the consistency condition using the estimated RE types Λ m and their observable times T m based on last event's observed time τ − m . The condition If it is consistent, the EM starts the G m evolution using identified consistent SCs for further RE estimation. If t i is shared, then t i 's forward SCG G n , s.t. t • i ∈ P n , is also reupdated based on the candidates of moving tokens from G m to G n . 4) If it is inconsistent, the EM activates an inconsistency diagnosis. If an AE arrives, the EM does not check the AE's consistency because users (in the CF) can generate the AE at any time in their purpose. The AE only reupdates its associated SCG based on identified consistent SCs. If the AE's transition is shared, then two SCGs are reupdated considering possible moving tokens.
For DT consistency checking, an inconsistency can stem from a TCPN error or any physical change or error, so cross-validation is required to diagnose the actual reason. We categorize the main reasons as follows: 1) a modeling error in the TCPN structure; 2) a new physical exception; 3) an operation-logic change that requires a TCPN revision; 4) a statistical change of an operation period; 5) a packet loss. Then, the EM initiates a manual diagnosis by looking into subsystems related to the inconsistent subnet.

A. SCG Evolution
In this section, we describe the SCG evolution procedures in detail, as depicted in Algorithm 1. A subnet N m 's SCG G m is computed based on root SCs, which are consistent with the recent physical timed event, meet the root property (see Definition 4), and denoted by S root . Root C k has one initial path π [1:1] = C (1) (= C k ) in its new path set Π + k and its emptied previous time-spent series, i.e., I(π [1:1] ) = , in Φ + k (see Definition 2). C k 's newly observable time range is zero, i.e., I acc (π [1:1] ) = [0:0]. SCG G m 's full path set Π m is initialized as C k ∈S root Π + k . Π m 's time-spent intervals, Φ m , consist of emptied time-spent information for each path in Π m . S root is stored in a set S wait that consists of waiting SCs (to be evolved) in the insertion order. Each SC C k = (M k , Θ k ) in S wait is evolved, as depicted in Algorithm 1.
First, the EM attempts to identify the set of imminent transitions (which can be fired first), denoted by T imm , among candidates T k as Line 4. T k is a set of logically enabled transitions, except for T ae (whose firing is triggered by physical AEs) and the previously identified imminent transitions (which are denoted by T − imm and specified on forward edges). Then, time-bound u min , s.t min{u • i | θ i = (·, u • i ) ∈ Θ k }, is computed to examine the following necessary proposition.
Proposition 1: Transition t i ∈ T k s.t. l • i > u min cannot be an imminent transition.
Proof: Let t j be a transition s.t. • u j = u min . If all transitions in T k \ {t j } are not fired until u min , then t j is eventually fired first based on the STS assumption. Therefore, a first transition always occurs within u min .
As (1) shows, each upper bound u • i can include any previous SCs' time-spent variables {Δ (-) }, and the intervals vary depending on the path, which changes u min and T imm accordingly. Thus, the necessary condition to be imminent is checked for each transition t i ∈ T k based on Proposition 1, considering all new previous paths in Π + k , as Lines 6-7. From (1), l • i (π If a previous path, π [1:n] a ∈ Π + m , makes transition t i imminent as Line 7, π [1:n] a is extended by t i -induced subsequent SC(s) from C k . In addition, SC C k adds π ) }, to save the path candidates for their extension. As mentioned in Section IV-C, the set of full paths, Π m , and its set of time-spent intervals, Φ m , are important to deriving the next RE's virtual candidates. If π [1:n] a ∈ Π m (i.e., |π a | = n), then π a and I(π a ) are removed from Π m and Φ m as Line 10. After examining all new paths using their own time-spent information, C k clears (Π + k , Φ + k ). Let S def be the set of deficient SCs. C k updates S def as Lines 12 and 13.
Using each new imminent transition t i , we derive the candidates of tokens inflowing to t i , {X j }. Then, the EM finds the next SC, C = (M , Θ ). Thus, M is derived as Line 16, and Θ is computed by considering newly disabled transitions (as Line 17), the time spent of C k (as Line 19), and newly enabled transitions (as Line 21).
If C is reusable using a function f eq , then we reuse existing C (in G m ) by making an edge e i from C k to C . If not, a new node for C is connected to C k by edge e i . If C is reusable, then edge e i 's label is extended by adding f eq . Unlike in the TPN-based SCG evolution, there can be multiple candidates for leaving tokens at the transition t i firing (i.e., |P re(M k , t i )| ≥ 1), so an SC can have multiple forward edges for the t i firing, as in SC C 3 in Fig. 6.
Using each new edge (caused by each transition t i ) from SC C k to C , C k computes the set of new full paths (denoted by Π • ) using the saved paths in D i . C is one of a newly added or a reused SC. If C was newly added, Π • is set by the If SC C is reused, C can already have its subsequent paths. Then, we extend the previous paths saved in D i based on C 's subsequent paths. We denote C 's subsequent paths as Π sub = {π sub | π sub = C . . . C ( ) ⊂ π a ∈ Π m , = |π a |}. A full path π a ∈ Π m can lead to multiple subpaths (reaching π a 's terminal) when C exists multiple times in π a . Using Π sub , we represent the new path set by Π • = {π [1:s] When merging a previous path π [1:n] a and a path π b ∈ Π sub , an SC C ( ) ∈ π b might not be accessible from π We denote Π • 's corresponding time-spent intervals as Φ • = {I(π a ) | π a ∈ Π • . The derivation method of Φ • based on the calculation of I(π a ) will be discussed in the following section. Based on (Π • , Φ • ), G m 's full path information (Π m , Φ m ) is updated as Lines 28. In our work, the next SC C evolves if the fired transition t i is silent. If t i is observable, C evolves after t i is confirmed as being consistent with the next RE and after being revised to satisfy the property in Definition 4. If C was not reused and scheduled in S wait for further evolution, new path information (Π • , Φ • ) is delivered to C as Line 32. If C was reused, C k traverses subsequent deficient SCs on new full paths to schedule them as Lines 34-36. This active scheduling of deficient SCs can avoid inefficiency in iterative examinations of all deficient SCs whenever there are no new SCs waiting for further evolution, as in [5]. Finally, a scheduled SC starts a new evolution.

B. Estimation of Next Reactive Events
The proposed approach utilizes full path and time-spent information of SCG G m , (Π m , Φ m ), to estimate the next timed RE. Let e ( ) a be the th edge on path π a ∈ Π m , connecting C ( ) and C ( +1) , where 1 ≤ < |π a |. Then, the set of the next REs' types from SCG G m is given by (5) Let Π m (λ) be the set of paths in SCG G m , leading to specific RE type λ. Based on subnet N m 's recent event-arrival time, τ − m , the expected observable time range of the λ-type RE, T m (λ), is given by Thus, the expected observable time of any RE related to N m this, we reformulate the lower and upper bounds as follows: where i is the transition index in e a , and so forth) to be recursively calculated as constant bounds.

A. Partitioning of TCPN
As discussed in Section IV-B, the sizes of TPN-based SCGs in [4] and [5] increase exponentially with TPN system complexity. Compared to the TPN-based evolution, the size of the proposed TCPN-based SCG grows more rapidly than that of a TPN-based SCG because each imminent transition t i of each TCPN-based SC C k can lead to multiple subsequent SCs-as many as Pre(M k , t i )| (see SC C 3 in Fig. 6). Thus, the speed up evolution of TCPN-based SCGs is essential for the runtime analysis. Hence, we propose partitioning TCPN to speed up the SCG evolution, described as follows.
We denote the set of the backward and forward transitions of place p as • p and p • , respectively. Let us look at the example in Fig. 7(a), where • p ⊂ T ob = T re ∪ T ae . Suppose that two subnets of N 1 = (P 1 , T 1 , . . .) and N 2 = (P 2 , T 2 , . . .) have the following properties: P 1 ∩ P 2 = ∅ and T 1 ∩ T 2 = • p . Based on the SC scheduling condition of Algorithm 1, any path from a root in an SCG is terminated when an observable transition is encountered. Thus, N 1 's SCG G 1 cannot affect the evolution of N 2 's SCG G 2 because any path reaching • p is stopped. Let T (i) ae and T (i) re be the sets of subnet N i 's AE-and RE-related transitions, respectively. Then, G 2 can also evolve independently to estimate an RE related to (T 2 ∩ T re )\ • p after assuming • p to be an unpredictable transition in T ae (i.e., • p ∈ T (2) ae ), as Fig. 7(b) shows. G 2 is influenced by G 1 only when • p -related physical events are observed, which triggers the G 2 reupdate. Generally, the set of subnets N from an original TCPN N 0 = (P 0 , T 0 , · · ·) must meet the following definition: Definition 5: The set of subnet N = {N }, where N = (P , T , · · ·) is the th subnet, 1 ≤ ≤ n, and n = |N |, should meet the following properties: To compute the subnets, we initially make temporary and minimum-sized subnets using each place in P 0 of the original net N 0 : e.g., N = {N } s.t. |N | = |P 0 |; |P | = 1, and T = • p ∪ p • , where p is the th place in P 0 . Each subnet iteratively merges other subnets in N when two paired subnets share any silent transition. If there are no longer any subnets to merge, then the remaining subnets are the final subnets.

B. SCG Update Based on Reactive Events
When detecting a new λ-type event at time τ (0) , the EM selects the SCG G m whose subnet N m contains transition t i s.
re . Then, the EM decides the consistency by checking (λ ∈ Λ m ) ∧ (dτ ∈ T m (λ)) from (5) and (6), where dτ = τ (0) − τ − m , not one subsequent SC in the TPN-based SCG. If it is inconsistent, the EM triggers the inconsistency diagnosis for the subnet. If it is consistent, the EM finds root SCs for the further SCG evolution by following two main steps.
1) Select the paths whose terminals are consistent SCs {C k = (M k , Θ k )} with the timed RE. 2) Revise each time constraint θ i ∈ Θ k of each consistent SC C k by resolving previous SCs' time-spent variables {Δ ( ) } to satisfy the property of Definition 4. Let Ω be a set of the pairs of consistent SC and its path. Then, Ω is represented by From a physical event, if we can get more information than type λ, which helps distinguish the processed tokens (annotated in the last edge e (n-1) ) based on their data, we can reduce the number of consistent SCs additionally, which lessens the evolution overhead.
For each (C (n) , π a ) ∈ Ω, where n = |π a |, the previous timespent terms in each time constraint θ i can be simply changed into a constant interval from ). However, we should consider one additional constraint, dτ .
Suppose that consistent path π a , s.t |π a | = 5 and (·, π a ) ∈ Ω, meets the following properties: I acc (π , respectively, which makes the interval of Δ (4) smaller than I acc (π [4:5] ). In a similar manner, dτ = 28 reduces the interval of Δ (4) to [4 : 4]. Considering the dτ constraint, we can generally formulate the interval of d i =1 Δ (n-) based on dτ and π a , as follows: From (9) and (10) Then, renewed consistent SCs are saved to the root set, S root , as initial SCs for the next G m .

C. SCG Update Based on Actuation Events
When observing a λ-type event at time τ (0) , the EM selects one or two SCG {G m } s.t. t i = L -1 (λ) ∈ T (m) ae ; two SCGs can be found when t i is shared by two subnets. Similarly, in Section VI-B, the EM aims to find the consistent SCs and their paths Ω based on dτ (= τ (0) − τ − m ) to derive the roots S root for the next SCG. Compared to the RE-consistent SCs (which are the terminals of specific paths in Π m ), any SCs {C (j) = (M (j) , Θ (j) )} on paths s.t. C (j) is observable at dτ and t i is logically enabled at M (j) are consistent SCs. Thus, we represent the pairs of consistent SCs and their belonging paths, Ω, by Ω = {(C (j) , π . If an SC C (j) is a path π a 's terminal and its last transition is silent (i.e., t (j-1) ∈ T si ), then u (j) a is ∞. If t i is shared and a subnet G m = (P m , . . .) has t i 's forward places (i.e., t • i ∈ P m ), then t i 's logical enabling is not checked. After deriving Ω, the EM performs the following two main steps to obtain the roots.
1) Renew Δ-dependent timing constraints of consistent SCs in Ω for Definition 4. 2) Update the markings and add new timing constraints of renewed SCs based on the t i firing. Compared to the renewals of RE-driven timing constraints, we should consider the elapsed time of each SC from previous transition t (j-1) up to τ (0) . For each (C (j) , π ) ∈ Ω, we add new jth time-spent variable Δ (j) to each constraint θ i ∈ Θ (j) , which symbolizes the interval of C (j) elapsed time. This leads to φ . Then, constraint θ i becomes Δ free using (9) and (10).
After the constraint renewal, the EM begins the marking update. If t i is not a shared transition, the marking of each If the size of {M } is larger than 1, multiple SCs are generated from C k and their constraints are also updated based on their markings.
If t i is shared by two subnets G m and G n , s.t. • t i ∈ P m and t • i ∈ P n , two subnets exchange token candidates inflowing to t i from G m . Let S (m) con and S (n) con be the constraintrenewed consistent SCs from G m and G n , respectively. If t i ∈ T ae , then each C k ∈ S (m) con is revised for the next markings {M | X j ∈ Pre(M k , t i ), M = M k X j }. Then, each SC C r ∈ S (n) con is revised or duplicated based on the next mark- and L(t i )-type RE arrives, each SC C r ∈ S (n) con is revised in the same manner in which the token candidates between the t i -shared subnets are exchanged.

VII. CASE STUDY
We applied the proposed approach to maintaining a DT of a factory located on the Singapore Institute of Manufacturing Technology (SIMTech) that has produced various research prototypes, such as customized USB devices [23]. The USB device requires a customer-selected logo and emits a scent, such as lavender, from a scent-filled cartridge when the device is active. Overall, the manufacturing processes consist of 17 steps: Steps 1 to 13, which are common to all products, and the following optional steps (including Step 14 for customer-specific sentences), as Fig. 8 shows. When a customer's order arrives, the CF generates an AE for specific parts of trays to be processed by custom steps. An operator in the CF occasionally generates an AE for some parts of a tray to complete specific steps (e.g., Steps 3-6) to reduce the makespan of the final product. If a tray does not undergo a process, it stays in a work-in-progress (WIP) conveyor or a machine's output dock(s).
Twelve machines are used for the production, and some machines, such as an infrared oven and a screen printer, are involved in more than one step. The processed trays in a machine's output dock are transferred by an automatic intelligence vehicle (AIV) to a destination's input dock or WIP. The TCPN for the factory consists of 220 places and 381 transitions, representing the tray unloading/loading by AIV, single-or multiple-part-specific processes of each machine, the delivery request and its grant, AIV movements along paths, and so on. The tokens represent the actual trays (for the part exchanges between machines), parts (for machines' processes), and synthetic control units (for the delivery requests, machines' availability, etc.). Each tray can carry up to 12 different types of parts. The machines, excluding the oven, process the input tray's specific parts based on the tray token's data. For example, if four parts of a tray token are required to be processed based on the tray token's part-related data, a Post function changes the tray token to four part tokens, as Fig. 9(a) shows. In the case of an injection-molding machine, depending on the part type (whether it is a reservoir or a bottom cartridge), its processing time varies; thus, part tokens' colors change for their subsequent processes. The part tokens would be merged into a tray token to be moved. As illustrated in Fig. 9(b), the target TCPN has common subnets for a tray-delivery request and its grant using an AIV. In the net, feasible requests are identified considering the availabilities of destinations' input docks. Among the tokens symbolizing feasible requests, a top-priority token is selected based on the tokens' data.
Using a daily RE/AE-observation history, we evaluated the DT-checking performance by varying the degree of the RE observation (OB). The degrees are categorized as follows.
1) OB1: All REs for all machines' tray process starts and trays' arrivals at the input and output docks. 2) OB2: All REs related to OB1, process endings/failures, and delivery requests. 3) OB3: All REs related to OB2 and AIV's idles/failures. OB3 represents all possible REs from the factory. Examples of silent transitions include order-related setups and  subprocesses of machines and AIV's multiple attempts (with extant maximums) to arrive at specific positions. The transition durations are measured manually or can be calculated based on observed events. Due to the stochastic property of the transition duration, there can be more than one next RE candidate with different time occurrence intervals.
We implemented the environment for the proposed TCPNbased modeling and SCG-based runtime checking using C++. If we intentionally create a DT inconsistency by considerably lowering the average of a specific duration, such as changing an infrared-oven operation from 2987 to 2680 s, an inconsistency can be detected, as Fig. 10 shows.
Aside from OB1, OB2, and OB3 enabled us to split the TCPN into 16 subnets for each of the machines, WIP conveyors, and AIVs, based on the partitioning condition (see Definition 5). Among the subnets, a subnet for AIV and delivery decision, with a place size of 145, is dominant. The overall workload consists of SCG evolution and finding the consistent SCs based on each RE/AE. We measured the performances using a machine with an Intel Xeon E5-1650 3.6 GHz CPU and 32 GB of memory. The results are displayed as Fig. 11 shows.
As the degree of event observability decreases (from OB3 to OB1), the length and number of paths created by firing possible silent transitions subsequently increases, which increases the numbers of the next RE candidates and computational overhead. Compared to the previous SCG evolution method in [5], the proposed SCG evolution method that supports active scheduling (for path-changed deficient SCs) and new path management (which prevents redundant computations in examining deficient transitions), thus decreasing the number of iterations.
TCPN partitioning requires additional consistent SC derivation and SCG evolutions whenever an RE related to a shared transition occurs, which consequently increases the overhead of consistent SC derivation. However, the average SCG size becomes significantly smaller, which drastically reduces the overall computation times. In our OB1 experimentation, the new evolution method reduced the overhead by about 8%. In the OB2 and OB3 cases, the final speedups (which are the ratios of the conventional evolution times to the active-scheduling-based evolution times after partitioning) were 18.97(= 18.21/0.96) and 8.28(= 5.69/0.68) times, respectively. Future work will discuss the proposed work's scalability by measuring the SCG evolution times and their speedups caused by TCPN partitioning, varying the plants' scales.

VIII. CONCLUSION
This work presents a novel approach for checking DT's consistency (for its high fidelity) by comparing streaming physical events with corresponding DT-based virtual estimates during runtime. For that, we aimed to build the virtual estimates of the plant's each physical event based on reachable observable transitions (that can produce physical events), considering all possible sequences of intermediate unobservable transitions, stochastic state transitions, stochastic operation periods, and external interventions. The previous TPN-based reachabilityanalysis methods facilitate deriving the reachable transitions and states by identifying all possible sequences of transitions based on their time ranges. To model the plant's high-mix manufacturing using CPN's colors and allow unpredictable users' interventions (AEs), we first designed the TCPN formalism. Then, we extended the previous SCG evolution method to 1) resolve some computational inefficiencies; 2) consider the TCPN formalism (from TPN); 3) alter the evolution stopping criteria sufficient to build the next event's virtual estimates. For better SCG evolution speedup and problematic-subnet isolation, we proposed a TCPN-partitioning method. We also proposed an iterative SCG update based on the previous SCG and each observed physical event for the next-RE estimation. We applied the approach to a USB device factory. Under a given experimentation set, the TCPN partitioning and revised SCG evolution algorithm accelerated the performance of the state-of-the-art evolution up to 18.97 times.