Automatic Model Generation and Data Assimilation Framework for Cyber-Physical Production Systems

Automatic model generation is defined as a method of creating simulation models from external data sources using interfaces of the simulation system and algorithms  [13]. These methods can be classified into three main categories [2, 10]:

• Parametric approaches: Models are initialized and configured based on parameters,
  
• Structural approaches: Model generation is based on data describing the structure of a system,
  
• Hybrid-knowledge-based approaches: combining existing methods (expert systems, neural network) with two of the above presented.
  

Our proposed framework applies an alternative approach by extracting the parameters and structure from the data to generate the model.

Son and Wysk [26] presented an automatic simulation model generation using a shop floor resource model and control model. The control model is constructed logically from a connectivity graph of the resource models. A component-based model generation with a data-driven approach is proposed by [10]. Data analysis results are used to extract the structural and behavioral preconditions as the basis for constructing the model and initializing the model state. The model is adjusted in an iterative process where the data from the real system and the simulation output data are compared and analyzed. However, these approaches require data describing processes and the topology of the system. By utilizing process mining techniques, our approach can extract the system structure from the data.

Krenczyk et al. [13] proposed a semi-automatic concept for model generation based on a parameter-based approach. Parameters will be used to generate different variants of production. Mages et al. [19] proposed a proof of concept to automate and hence reduce the cost of the process of simulation model creation, thereby allowing for the creation of a larger number of selectable solution variants. They used component-based synthesis using combinatory logic to synthesize a product line of varying simulation models for a given configuration to be executed and evaluated to find suitable solutions. However, these parameter-based approaches require the definition of the simulation model, which will be used to generate different variations of the production. Other than automatically extracting parameters from the data, our proposed approach is also able to generate the complete process model.

Friederich et al. [6] proposed the use of process mining techniques to support the extraction of reliability models from event data generated in smart manufacturing systems. More specifically, a stochastic Petri-net is extracted to analyze the overall system reliability as well as to test new system configurations. An online data-driven approach is not considered in their work; only an offline approach is used to extract the process model for further predictive analysis.

After generating the process model, it is important to validate the generated model against the system. This can also be used to select the best-estimated model during the data assimilation process. We focus on operational validation, which is defined as “determining that the model's output behavior has a satisfactory range of accuracy for the model's intended purpose over the domain of the model's intended applicability” [23].

Lugaresi et al. [17] presented a new method to perform online validation based on harmonic analysis which is able to achieve reliable results even when applied to a relatively small dataset. The proposed method uses information about the power spectral density of KPI trends as signals to detect potential deviations between the model and the real system. Morgan and Barton [21] proposed a methodology to discriminate between two systems based on the Fourier transform of the output trajectory. By analyzing the output in the frequency domain, the method can consider the temporal evolution of data, rather than comparing aggregated values such as the average. Lugaresi et al. [18] proposed a new technique to deal with a limited data set and analyze the evolution in time of the behavior. The simulation model validity is assessed by measuring the similarity between real system data and simulation model data using a sequence comparison technique. We proposed to use coefficient determination of the historical output between the real system and the simulation model to assess the fitness of the model. However, a detailed analysis of our approach compared to the existing methods has not been investigated in this paper. This will be one of the future directions for further investigation.

For discrete event simulation, conventional data assimilation methods, such as Kalman filters [28] and Extended Kalman filters [8], are inapplicable due to nonlinearity in the model, and non-Gaussian property of noise. Sequential Monte Carlo (SMC) methods, also called particle filters, can approximate arbitrary probability densities and have little or no assumption about the properties of the system model [1]. Particle filters estimate the posterior distribution of a state trajectory by a set of random samples with associated weights and compute estimates based on these samples and weights. As the number of samples becomes very large, this becomes an equivalent representation to the usual functional description of the posterior distribution [5, 28].

Hu and Wu [9] proposed a data assimilation framework for discrete event simulation. A mechanism is designed to synchronize the system state of a discrete event simulation with the data assimilation time; otherwise, the system state at the data assimilation time will be inaccurate. To address the low diversity of samples due to the deterministic or small degree of stochasticity in the discrete event model, sample rejuvenation is used to add stochastic perturbations to the system states represented by samples after the resampling step.

Xie and Verbraeck [30] presented a particle filter-based data assimilation framework for discrete event simulations. Instead of synchronizing the system state with the data assimilation, they proposed the interpolation of the discrete model states to obtain updated state values. The sequential importance sampling algorithm is applied to update the samples that approximate the posterior distribution of the state trajectory with variable dimensions.

However, the SMC method assumes that model components do not change over time (i.e., closed systems). Our proposed approach sequentially generates new models with the arrival of new data, hence it is not possible to directly apply the existing methods in our framework. The proposed approach can also deal with the diversity of samples. Variations in the samples are introduced when each sample is generated based on a different snapshot of the events and the initial state of the generated model is also randomly generated.

A general production system is modeled with a queue before each resource, as shown in Figure 1. Resources represent the Automated Guided Vehicles (AGVs), assembly cells, or various machines and tools in the factory. Each order will be processed by a sequence of resources in the system according to the process recipe. When the order arrives at the resource, it enters the queue waiting to be processed, emitting the Arrive event. If the resource is idle, the order in the queue will be processed, emitting the Produce event. Depending on the capacity of the resource, it is possible to process one or more orders concurrently. When the order is processed by the resource, the Route event is emitted to indicate that the order is completed and will be routed to the following resource for processing. Reentrant manufacturing is also supported where the orders can repeatedly pass through the same resource at different stages of the process routes.

The manufacturing system's data is collected and distributed to the company's information systems. Several information systems in the company may collect information about the physical system, such as Supervisory Control and Data Acquisition (SCADA), Programmable Logic Controller (PLC), Manufacturing Execution System (MES), or Enterprise Resource Planning (ERP) [6]. The data captured by such systems is then aggregated in event logs or used to generate key performance indicators (KPIs). A manufacturing KPI is a well-defined measurement to monitor, analyze and optimize manufacturing processes regarding their quantity, quality as well as different cost aspects. They give manufacturers valuable business insights to meet their organizational goals.

The general assumptions [15] about a process are slightly modified to fit the manufacturing domain:

• A manufacturing process is triggered by production orders.
  
• An order is characterized by a sequence of activities that begin and end by events.
  
• Each event has a corresponding timestamp and resource that is utilized and released.
  

Considering the above-mentioned assumptions, an event log EL is defined as a set of event entries:

Each event log entry is defined as a tuple e_i = (t, o, r, c), where t is the timestamp, o is the set of order identifiers, r is the resource identifier, and c is the event class. Order identifiers are unique strings representing individual production orders. Event classes are unique strings to mark the basic operations in a production system, i.e., c ∈ {Arrive, Produce, Route}.

Typical event logs, used by mainstream process mining techniques, require an event to be related to a case or an object. However, batch-processing machines can process several orders as a batch at the same time. To aggregate data regarding batch processing in the event logs, object-centric event logs (OCEL) are used to relax the assumption that an event is related to exactly one case [7]. Hence, o = {o_j|j = 0, 1,...} so that an event can be related to multiple orders to represent the batches.

A Petri-net model is used to describe production process models extracted from the event log. It is a class of discrete event dynamic systems, which can be used to simulate discrete event systems, such as the production process. Formally, a Petri-net is defined as a tuple N = (P, T, A) where

1. P and T are disjoint finite set of places (drawn as circles) and transitions (drawn as rectangles), respectively.
   
2. $A \subseteq (P \times T) \bigcup (T \times P)$ is a set of directed arcs.
   

Each transition is defined as T_i = (c, r), which corresponds to events (t, o, r, c) in the event log. The transition types are modeled after the components in the system. Transitions are labeled as a tuple (r, c) for the resource r and the event class c. Figure 2 shows a Petri-net model of a single resource that contains the following components as the transitions: Arrive, Produce, and Route. The Arrive transition represents the arrival of a token at an unbounded resource queue in the production process. An inhibitor arc is used to model the resource capacity by preventing a transition from firing based on the number of tokens in the place. The timed transition at the Produce represents the production time of the resource. Route is represented by an immediate transition that routes the token to the next resource. By simulating the Petri-net model as a discrete event simulation, it is possible to measure various system performances, e.g., KPIs such as the cycle time or throughput of the manufacturing system.

A set of data and process mining techniques are used to automatically generate the process model based on the system components. First, process discovery is used to generate a process model that generalizes the observed behavior from the event data without any apriori information [27]. The discovered process model may describe only the control flow of events or other operational aspects, e.g., processing time. Next, process enhancement using various data mining techniques is performed to improve the current process model to achieve a more complete model behavior from the events, e.g., queuing or batching operations [24]. Figure 4 shows the main steps in the data and process mining.

1. The incomplete processes are trimmed from the OCEL.
   
2. Object-centric directly-follows multigraphs (OCDFG) are generated by discovering Directly-Follows Graphs (DFG) from the OCEL.
   
3. Process discovery generates the Petri-net model from the flattened EL.
   
4. Production time distributions are fitted based on the activity durations from the EL.
   
5. Resource capacities are extracted from the OCDFG.
   
6. Decision flows in the Petri-net are constructed using decision tree classifiers based on the attributes of the events in the event log.
   

The steps will be explained in further detail in the subsequent subsections.

4.1.1 Removal of Incomplete Processes. For an online data assimilation procedure, the sampling of the events may occur for any sequence of events. In Figure 5, there are four orders in which each order is processed in sequence. Hence, the complete process for each order will be the following events: e₁ → e₂ → e₃. Orders o₁ and o₂ are only able to capture the incomplete processes, while orders o₃ and o₄ managed to capture the complete processes. Process mining on the partial processes will generate an incorrect model as the process does not begin on events e₂ or e₃. To remove the incomplete process events, the frequencies of the occurrences of the first events grouped by the resource and event class are determined from the event logs. As the event log aggregate more complete processes, the frequency of complete processes will be higher than the incomplete processes, i.e., highly frequent first events indicate complete processes while infrequent first events are incomplete processes. Hence, by filtering infrequent first events from the set of first events, incomplete processes can be removed from the event log before extracting the process model based on the event log. For example, the frequency count of events {e₁, e₂, e₃} is {2, 1, 1} in Figure 5. If we filter away infrequent events e₂ and e₃, event e₁ is left, which represents the first event of a complete process.

The behavior of discrete-event simulation is highly non-linear, non-Gaussian [31]. SMC methods can be applied to discrete-event simulation since they are able to approximate arbitrary probability densities and have little or no assumption about the properties of the system [1]. However, the sequential sampling process assumes a closed system, and cannot handle structural changes in the system [30].

A new data assimilation method is proposed to handle changes in the system, as shown in Figure 6. A set of samples $\mathcal {S}$ is used to approximate the system, where each sample consists of a generated simulation model, the corresponding simulation states, and the historical measurements of the model. At each sampling step, the events from the real system are used to generate a new sample. The sample is simulated to obtain the measurements. Instead of comparing the likelihood between the observations from the real system and simulated measurements at each sampling step, the coefficient of determination (R²) between the historical observations and measurements is used to evaluate the fitness of the generated models. R² measures how well the observed outcomes are replicated by the model, based on the proportion of total variation of outcomes exhibited by the model. The proposed method then continuously generates new samples with the arrival of new events and observations, while discarding samples with low fitness.

For example in Figure 6, the event log EL₁ is used to generate sample s₁, which is simulated until time t₁ to obtain measurements ML₁. The fitness of s₁ is evaluated by comparing the measurements ML₁ with the observations OL₁. In the next sampling step t₂, the new events E₂ are appended to EL₁ to obtain EL₂, which is used to generate a new sample s₂. The samples s₁ and s₂ are simulated until time t₂ to obtain the corresponding new measurements, M₁ and M₂, which are appended to the respective historical measurements ML₁ and ML₂. At the same time, the new observations, O₂, from the real system are appended to OL₁ to obtain OL₂, which is then used to evaluate the samples s₁ and s₂ by comparing to the historical measurements ML₁ and ML₂. The same procedure is repeated for every sampling step. As there is a limit of 3 samples, at time t₄, the sample with the lowest fitness, i.e., s₁, is discarded.

The case study utilizes a scaled-down representation of a real-world wafer-fab at Harris Semiconductor described by [12]. This model was developed for a discrete-event simulation, which is an important tool for decision-making in the wafer fab [20]. It is based on an actual wafer-fab factory and captures the main complexities of the actual systems while remaining compact enough to be simulated in a reasonable amount of time.

Lots are the moving entities in the wafer fab. Each lot contains a fixed number of wafers. The lots are generated based on customer orders. In the wafer fab, the lots are processed in a series of process steps using multiple pieces of equipment.

A single wafer fab is composed of several work areas. These areas are used for wafer processing and sorting in a clean room environment. A single work area typically contains dozens of different work centers. The work centers of a work area are closely related logically or physically. A single work center is a set of machines that provide similar processing capabilities. Work centers are also called machine groups or tool groups. Some tools can process only a single wafer at a time, while batch processing tools can process a set of lots simultaneously.

The wafer fab is modeled as a set of event-driven model components whose state changes are triggered by events [25]. The time unit used in the simulation is minutes. The typical model components consist of machines and wafer-lot queues. The lot is modeled as a simulation event, and the process flow describes how a lot is routed among the machines. Each simulation event contains information regarding the lot name, the number of wafers, and process recipe, etc.

Each machine model contains input and output ports to receive simulation events and send the events to the next machine model. When the machine model receives a lot, it performs a specific process according to a recipe of the lot that determines the operation parameters such as the setup time, processing time, etc. After processing, the lot is transferred to the next machine in the process flow. Each machine has a wafer-lot queue, which is a buffer to store the incoming lots. The waiting lots in the buffer are dispatched to available machines based on pre-defined dispatching rules, such as the First-in-First-Out (FIFO) dispatching rule. The lots are also assumed to be instantaneously transferred between successive process steps. By coupling the machine models using input and output ports, the modeler can build various types of wafer-fab systems using modular model components.

There are 11 machine groups and 3 order types in this simulation model, as shown in Figure 7. The number of operations for Order Types 1, 2, and 3 is 22, 14, and 14, respectively. The bottleneck station of the fab (Station 4) is visited repeatedly, which corresponds to the photolithography process. This bottleneck is mainly due to the costly and scarce resources for the photolithography process.

The wafer-fab simulation is used as the system to generate the events and observations. The wafer lots and machines represent the orders and resources respectively in the production system. Every time a wafer lot enters a wafer-lot queue, the Arrive event is emitted. When the waiting lots are dispatched to the machine for processing, the Produce event is emitted. When the lot finishes processing at the machine, a Route event is emitted to signify the routing process to the next machine. The cycle time of a wafer lot is the total time required by the wafer-fab to process the lot according to the recipe. The cycle time is used as the performance measurement of the wafer-fab as it is one of the KPIs in wafer-fab manufacturing.

To validate the generated model, a wafer-fab simulation based on the case study is executed for 300 Order Type 1, 100 Order Type 2, and 3 respectively based on the frequency defined in the dataset. The events generated are aggregated in an event log, and the cycle time for all the orders is recorded. A process model is generated from the event log and simulated with the same orders as the original simulation model. Figure 8 shows the Petri-net model that is automatically generated from the whole event log. The cycle time of the orders in the process model is also recorded for comparison.

The observations (average cycle time) from the system are compared to the measurements of the generated model. Each point in Figure 9a shows the cycle time of the same order between the observations (x-axis) against the measurements (y-axis) for Order Type 1. A linear regression is performed between the observations and measurements, as shown in the solid line in Figure 9a . It can be seen from the regression line that there is a negative correlation between the observations and measurements. In addition, the coefficient of determination between the observations and measurements is 0.45. This means that for the individual orders, the measurements do not replicate the observations.

It is usually not feasible to determine if a model is absolutely valid over the complete domain of intended use or application [23]. Hence, the distribution of the system output is compared instead of individual order cycle time. The cycle time of the observations and measurements are sorted in ascending order before being plotted in Figure 9b . The observation for a point only matches the measurement in order, they do not represent the same order. It can be seen from the figure that there is a positive correlation between the observations and measurements. In addition, the coefficient of determination between the observations and measurements is 0.99. Hence, the overall distribution of cycle time by the generated model is similar to the output from the system.

The results are also similar for Order Types 2 and 3, which are not shown here due to space limitations. We also performed the two-sample Kolmogorov-Smirnov test between observations and measurements. The p-values obtained for each of the order types are 0.78, 0.81, and 0.21. Choosing a confidence level of 95%, the null hypothesis that the two distributions of the observations and measurements are identical cannot be rejected as the p-values are more than 0.05. These results also illustrated that standard sequential importance sampling algorithms are not applicable as the individual orders do not show the same output, only the overall distribution of the cycle time is similar.

An identical twin experiment is conducted to evaluate the feasibility of the proposed data assimilation method. The wafer-fab simulation is executed for 21 days (30,240 minutes), with a warm-up period of 5 days (7,200 minutes). After the warm-up period, the new events and observations are recorded for the data assimilation process. The data assimilation is executed every 120 minutes, i.e, the sampling period ΔT = 120. The time window w is defined as 3 days (4,320 minutes). The number of samples N_s is defined to be 10.

The fitness of the samples is evaluated by comparing the observations and measurements using the coefficient of determination. The coefficient of determination of the samples over time is plotted in Figure 10. The color represents the age of the sample, where the sample with age 0 is the newest sample generated. Darker colors indicate older samples while lighter colors indicate newer samples. As the data assimilation procedure proceeds with time, samples with lower R² are discarded and newer samples based on the windowed event logs are generated. Hence, the color of the sample becomes lighter as time progresses.

In Figure 10, the R² are very far from 1 at the beginning of the data assimilation. Due to insufficient events in the event log, the generated models are a poor estimate of the system. As the number of events aggregated in the event log increases, the R² of the newer samples improves, before converging towards 1. There are still some newer samples that are far from 1 due to different simulation states randomly initialized for each sample. However, the samples still generally converge towards 1, indicating that the newer samples have similar output compared to the system.

To evaluate the applicability of the data assimilation framework to handle the system model changes, we emulate structural changes in the original system. A time window of 7 days (10,080 minutes) is used in this experiment. As shown in Figure 7, Station 4 is one of the bottlenecks in the wafer-fab factory which will impact Order Type 1 and 2. There are two machines for the Station 4 – Machine 4.1 and Machine 4.2. To emulate changes in the system, Machine 4.1 will break down at 7^th day (10,080 minutes). After breaking down, it will no longer be able to process any more lots. The breakdown will cause a bottleneck to occur in the factory, increasing the overall cycle time. A recovery operation will be performed at 14^th day (20,160 minutes) to increase the capacity of Machine 4.2 in order to clear the backlog orders and reduce the cycle time. Figure 11a shows the cycle time for each order type over time. It can be seen from the figure that the cycle time of Order Type 1 and 2 increased from the 7^th day due to the breakdown. The cycle time of Order Type 3 is not affected as it does not require Station 4 for processing. After the recovery operation has been executed at the 14^th day, this eventually clears the backlog of orders and reduces the cycle time.

Figure 11b shows the cycle time of each order type over time for one of the samples that are generated during the breakdown period. As Machine 4.1 is already broken down, there are no events emitted from this machine; only events from Machine 4.2 are aggregated. However, since the recovery operation has not been performed, the capacity of Machine 4.2 is not increased yet. Hence, the generated model only extracted the original production capacity for Machine 4.2. The original capacity of Machine 4.2 is insufficient to service all the incoming orders for Station 4, leading to increasing cycle time for Order Types 1 and 2, as seen in Figure 11b . The cycle time of Order Type 3 is not affected as it does not require Station 4 for processing. Due to the increasing cycle time for Order Types 1 and 2, this will lead low coefficient of determination compared to the observations.

Similarly, the coefficient of determinations of the samples over time is plotted in Figure 12. Initially, the R² of the samples increases to a maximum R² = 0.77 at 12,480 minutes. This indicates that the data assimilation is eventually able to generate models that are similar to the pre-breakdown system. However, the R² eventually decreases over time after 12,480 minutes due to the large divergence of the generated models compared to the post-breakdown system.

There is a set of newer samples with increasing R² in Figure 12 after 30,240 minutes. After the recovery and a time window of 7 days (10,080 minutes), sufficient events are aggregated in the event logs such that a new set of samples that can represent the post-recovery system is generated. The coefficient determination eventually converges to the ideal value of 1 as better-estimated models of the post-recovery system are generated. Hence, the proposed data assimilation framework can adaptively handle the changes in the system over time.