Maintenance optimization for a multi-unit system with digital twin simulation

Optimization of operations and maintenance (O&M) in the industry is a topic that has been largely studied in the literature. Many authors focused on reliability-based approaches to optimize O&M, but little attention has been given to study the influence of macroeconomic variables on the long-term maintenance policy. This work aims to optimize time-based maintenance (TBM) policy in the mining industry. The mine environment is reproduced employing a virtual model that resembles a digital twin (DT) of the system. The effect of maintenance decisions is replicated by a discrete event simulation (DES), whereas a model of the financial operability of the mine is realized through System Dynamics (SD). The simultaneous use of DES and the SD allows us to reproduce the environment with high-fidelity and to minimize the cost of O&M. The selected illustrative case example demonstrates that the proposed approach is feasible. The issues of using high dimensional simulation data from DT-models in managerial decision making is identified and discussed.


Introduction
Managing large industrial plants in global competition requires a clear strategic view and a high level of control of operations. Anytime an industry relies on its physical assets, the success of operations is tightly linked to the execution of the right level of maintenance. Maintenance has both the role of keeping an asset in its best condition and to minimize unforeseen system downtime. From a managerial viewpoint, operations and maintenance (O&M) cannot be thought in isolation from the economic context within which every industry operates: a coordinated view of O&M should aim at reaching the right amount of responsiveness and throughput of a sys-tem that is required by the prevailing market conditions. In this research, the issue is investigated using an example from the metal mining industry, where efficient real-time management of operations is essential to meet the production targets, but where ultimately macro-economic variables, mainly the price of the metal, play the key role in bottom-line profitability in the long-term.
As stated by Bevilacqua and Braglia (2000) and Mobley (2002), maintenance costs can rise to 60 % of total production costs. This cost item can be affected in the short-and mediumterm by planning and optimization -unlike many other major costs of industrial operations that are fixed in nature. Despite the importance to plan operations for the impact on long-term profitability, there is only a limited amount of literature on the topic, except for Topal and Ramazan (2010) who introduced a model to estimate maintenance costs in a 10-years mine lifetime. Furthermore, considering multi-machine environments, the sheer size, and the resulting complexity due to a high number of uncertainties is a major hurdle for model development (West and Blackburn 2017). Addressing a company-wide problem-setting, like managing real-time operations and maximizing long-term profitability in a dynamic economic context, requires the help of both advanced analysis methods and control tools. We address the topic using a Digital Twin (DT) modeling concept that is used here in a meaning discussed by, e.g., Rosen et al. (2015), Grieves and Vickers (2017), to refer to "interconnected and multidisciplinary simulation models usable for operations optimization on a system level". In a recent review of Kendrik et al. (2020) five use-categories of DTs were identified of which the manufacturing stage and usage stage of a system is addressed in this paper.
A DT is a digital model of a physical entity (Negri et al. 2017;Tao et al. 2018;Redelinghuys et al. 2020) providing human-readable, semantic, data-model of reality (Negri et al. 2017;Kunath and Winkler 2018). These models reside in a high-performance, usually cloud-based, computing environment and they can be used for several types of optimization purposes (Negri et al. 2017;Tao et al. 2018;Cimino et al. 2019;Madni et al. 2019). Kendrik et al. (2020) highlight the importance of the digital counterpart of a system to optimize production performance and maximize profitability, which is the goal of the proposed model. The origins of DT can be traced back to the beginning of the 2010s in the aviation and aerospace industry. The early publications (Tuegel et al. 2011;Shafto et al. 2010) revolved around the possibilities of using ultra-high fidelity models to simulate aircrafts' maintenance under dynamic operating conditions over the equipment lifetime. In this vein, Kritzinger et al. (2018) highlighted the communication aspect between physical and virtual spaces claiming that only models transmitting data in and out from the virtual space can be regarded as DTs. The latter should not be confused with general digital models (no data connection) or digital shadows (only physical-cyber connection). The required fidelity level in DT models remains debated and, in this research, we agree with the claim by Wright and Davidson (2020) that "digital twins can use any sort of model that is a sufficiently accurate representation of the physical object being twinned". This paper focuses on the question of building and utilizing multi-domain simulation models that would integrate O&M simulation optimization with the overall profitability simulation of industrial operations in a way that could be referred to as Digital Twin. For the sake of brevity, we limit our scientific inquiry to the context of the mining industry.
To answer the research question, a two-phase methodological approach is adopted. First, the general properties of a co-simulation framework are investigated, and references to the relevant literature are provided. Second, an experimental DT model is developed on a virtual case study: a metal mine is considered due to its specific nature of O&M, and because its profitability is directly linked to the price of metal(s), which is a macroeconomic variable. To answer the question, this study presents a conceptual "digital twin" for metal mining that connects a detailed, minute-per-minute maintenance model of mobile equipment to a monthly-level profitability analysis of metal mining operations. In the model, two separate simulation modules are included: an O&M model, which replicates with high fidelity the effects of O&M decisions, and a managerial cash flow (CF) model, which is used to support decision-making at the production system level. In a co-simulation context, both models are treated as separate simulation units (SU) and when these SUs are considered as a whole, a dynamic system is created (Gomes et al. 2018).
This allows us to replicate a DT model's operational workflow and software pipelining in a controlled environment, where the O&M model optimizes some of the key system parameters before running the CF simulation for the high-level mid-term economical aspects of the system. The complexity of the system under study is a major reason to adopt a DT-inspired view: where it is not possible to express relationships analytically, a DT can help to integrate data from the field with flexible simulation tools, to achieve an overall improvement of the system's profitability. Therefore, the goal of this work is to: -Demonstrate that the DT approach in the context of metal mining operations provides a holistic method to optimize its overall operational profitability under economic uncertainty of metal prices and maintenance costs. -Point out and discuss the limitations of simulation based digital twins, when it comes to managerial decision making based on multidimensional information.
This paper continues with a brief introduction to the concepts of O&M planning in multi-equipment systems and some general considerations about system dynamics methodology in Sect. 2. In Sect. 3, a literature study on the topic of O&M simulation and DT modeling is provided to set the ground for model building. In Sect. 4, a detailed description of the models-namely the O&M module and the CF modelis provided. This is followed by the empirical application of the model, the validation of the proposed model through two experiments, and a detailed analysis of the results in Sect. 5. The paper closes with conclusions and discussion in Sect. 6, where some strategic considerations are derived from the results of numerical experiments.

Theoretical Background
Machine specific maintenance histories can be tracked with high accuracy using data series of sensory information together with maintenance reports from existing databases. A window of opportunity exists to use this accumulated maintenance information in connection with the DT model depicting the behavior of the overall system. From the reliability-theory point of view, a large-scale industrial system can be mod-

Single-item System
System boundary Unit U 1 (a) A single-item system.

Multi-item System
System boundary

Fig. 1
A schematic representation of single-and multi-item systems eled using a multi-item system made of several non-identical components, which are characterized by a common set of features, but with proprietary parameters for each feature. Two common examples of such features are the service time, and the time to failure (TTF)-distributions. In a single-item system, which can be depicted as in Fig.1a, maintenance can be optimized knowing the TTF distribution and the cost of corrective maintenance. On the other hand, multi-item systems are sets of components considered as a whole, and they can be represented as in Fig. 1b. One peculiarity of multi-item systems is that very often there is a convenience to carry out maintenance simultaneously on groups of components: since component dependencies of different natures exist -i.e. economic, stochastic, or structural (de Jonge and Scarf 2020) -they can be exploited to minimize maintenance costs and system downtime.
In multi-item systems, maintenance activities and regular operations can be organized according to a maintenance strategy, which determines the rules for scheduling of both. According to Alrabghi and Tiwari (2016), there are two broad classes of strategies: time-based maintenance (TBM) and condition-based maintenance (CBM) strategies. Both types of strategies include the possibility to perform corrective maintenance (CM) and preventive maintenance (PM) interventions, where the latter kind of activities are justified by the lower cost of stopping the system and inspect/maintain components before they fail. From the economic point of view, skipping PM can save money in the short-term, but exposes to the risk of more expensive breakdowns in the mid-and long-term.
The major difference between TBM and CBM is the principle that rules decisions: to plan maintenance activities, TBM uses only the work time, whereas CBM exploits also information on the degradation of a component. Depending on the cost and the risk generated by the fault of an item, both strategies are valuable. Concerning multi-item systems, the state of the art for both types of strategies were reviewed several times in the past (Cho and Parlar 1991;Dekker et al. 1997;Wang 2002;Nicolai and Dekker 2008;de Jonge and Scarf 2020).
A recurrent critique of many multi-item models, which is partly addressed in this paper, is the lack of integration with other fundamental parts of an industrial system-e.g., spare parts and inventory management, human resources management, or planning of operations. Tiwari (2015, 2016) confirm this by stating that the isolation of maintenance management systems is a limit to their use in practice. The experiment design used in this research resembles the one proposed by Alrabghi and Tiwari (2016) for TBM but contextualized and integrated with higher-level decision-aid tools. The DT framework offers the right testbed for simulationbased production optimization (Uhlemann et al. 2017), and for studying the integration of systems, hence to overcome system isolation.
To deal with the model integration issue, the System Dynamics methodology, originally coined by Forrester (1961), is used in this study. SD is suitable for representing the behavior of complex systems with delays and feedback loops that are constructed using intuitive graphical flowsheet diagrams (Forrester 1994). Within engineering sciences, SD has been traditionally viewed as a high-level managerial method, which is subordinated to fast-to-run, disciplinespecific computational models; however, SD has also been applied in several operations research (OR) applications, which were reviewed by Größler et al. (2008). In this paper, the role of the system dynamic model is to serve as a semantic data interface to the overall production system, where all the relevant sub-model(s) can connect.
In this paper we focus our scientific inquiry to the context of metals mining, where the role of equipment reliability is highlighted by the complexity of advanced machinery, and the pressure to meet the production targets (see discussion, e.g., Dhillon (2008)). In real mining systems, data-driven analysis of maintenance policy optimization faces the problem of the reliability behavior of equipment. As a key challenge to maintenance, Hall and Daneshmend (2003) point out that the number of (semi-)mobile equipment hinders the collection of "clean" datasets. Data collection may also be inhibited by the failure of electronic-based hardware (e.g. sensors, wiring, connectors, etc.), which is common in harsh mining environments (Dhillon 2008). The estimation of the near-future degradation state of machines and the forecasting of their most likely end of life require the use of simulation, which is recognized as a main aspect of a DT (Negri et al. 2017;Kritzinger et al. 2018;Tao et al. 2018;Cimino et al. 2019). For these reasons, we consider our model eligible to operate as a DT, although the experiments that are presented in the following do not rely on a real-world physical sys-tem, and a proper product data management system is not implemented.

Literature Study
To clarify the connection of this work with the existing literature, a brief study on the topic of simulation-based DTs and maintenance was conducted. An overview of the models involved in manufacturing system design and operation using DES is provided by Negahban and Smith (2014), who observed that there is an on-going shift to maintenance issues and real-time control. In this vein, we used the following three combinations of keywords to conduct an inquiry on the search engine Scopus: i) "digital twin", "simulation", and "maintenance"; ii) "digital twin", "co-simulation", "maintenance"; and iii)"co-simulation" and "maintenance". Based on their relevance to this research, 30 documents were selected and listed in Table 1.
The columns "Digital Twin", "Maintenance", and "Cosimulation" are flagged if the keyword represents a relevant topic in the document. The columns "Review", "Methodology", and "Application Case" indicate if the document includes a review of the literature, a contribution to methodological aspects, or the presentation of a use case.
Based on the literature study, the number of publications concerning DTs and simulations for O&M optimization increased during the last ten years, as depicted in Figure 2. The majority of the published documents are represented by conference proceedings although the relative share of journal articles has been in a steady increase during the period of 2017-2020. This suggests that the relevance of the topic is being identified in the scientific community.
The content analysis reveals that most of the works (in Table 1) aim at developing technical models of mechanical, electrical, aerospace, and transportation systems, but only a few documents specifically addressed the combination of technical and economic aspects. There seems to be a common understanding that maintenance optimization has a central role in DT models, together with the general aim to improve operations and managerial prediction capabilities. The latter topic heavily relies on the simultaneous use of several simulation tools, but there seems to be little awareness of the co-simulation context that emerges. To verify this observation, our initial research query "i)" was tested by substituting the keyword "simulation" with "co-simulation": the low number of documents found suggested a lack of general frameworks when co-simulation models are part of a DT.
Outside the context of DTs, the principles and properties of a co-simulation model have been systematically surveyed by Gomes et al. (2018), who highlight the ability to apply separate, "black-box", simulation units as building blocks of a large (co-)simulation. This aspect is of particular impor-tance in the real world, where simulation tools for prognostic and health management (Peng et al. 2010;Kim et al. 2016) might come from different developers and they need to be integrated. Several documented industrial applications of cosimulation models within the period 2011-2016 are reported in Gomes et al. (2017).
The documents resulting from query "iii)" are similar to the references mentioned by Gomes et al. (2017) in their literature review. A closer look at these documents reveal that co-simulation models are often "stand-alone" works that do not present a connection with a physical model. In other words, although the potential of co-simulation models in maintenance optimization is clear, there is a lack of research efforts describing how these technical-economic models would be structured and how they would play out. This research work aims at contributing to close this gap by considering simulation optimization of O&M as part of a DT, and by addressing the issue according to the principles of co-simulations.

Data and Methodology
This research addresses the problem of designing a DT, which comes down to the ability to be able to simulate and optimize several models (co-simulation). Such models are not directly integrable due to their fundamental basis (such as software, modeling choice, and the level of detail), and they need to share information in an uncertain/probabilistic environment. Two models are co-simulated in this research: i) an O&M model, and ii) a managerial CF model, which operates high-level decisions based on generated CF and exogenous economic variables. The fleet capacity optimization is used as a means to achieve the economic goals: if there were to be two alternative fleets that meet a production target, the one producing the higher CF would be selected. A schematic diagram of the problem setting is illustrated in Figure 3, where the connections between separate steps are shown.
Inputs of the Digital Twin model consist of the system design and maintenance policy selection, which are marked with (i) and (ii) in Figure 3. These are used to feed the software module (iii) that replicates the operations of a metal mine modeled with a DES. The performance of O&M is evaluated by running a Monte Carlo simulation (MCS) according to the given design and policy selection. The maintenance model can be run with an optimization procedure automatically changing the system design -i.e. the number of resources in the mine, to reach a target value of ore excavated at the minimum cost. The aggregated information produced in (iii) is fed to the managerial feasibility model (v) with economic uncertainties included (iv). The aggregated system output is formed ((vi) in Figure 3), which can be, in case of operating the DT-model continuously, further looped back  The simulation ability of DT is discussed, and the concept of intelligent-DT is presented.
x x x   Table 1 divided per year. "2020" represents the year-to-date in August to the maintenance module to revise the policy until convergence of results is reached. Once an optimized maintenance policy is found, it can be used to control the physical system.
We acknowledge that the multi-disciplinary model applied in this paper is limited in nature, and we suggest that this model could be referred to as "low fidelity digital twin" (see discussion, e.g., Tuegel et al. (2011)) to distinguish them from the envisioned "full-scale" DT implementations including a wider range of simultaneously operating, high-fidelity, disciplinary sub-models.

Maintenance Module
The justification to use simulation-optimization to model the mine environment in this paper emerges from two reasons: the lack of analytical expressions to model operations, and the need to adapt the configuration of resources to meet the production targets. The module aims at optimizing the maintenance policy, which is a set of heuristic rules to make maintenance decisions. The inherent complexity of the sys-tem makes it impossible to determine in advance the effects of the proposed maintenance policy, therefore, O&M of a mine's load and haul process is replicated using DESs in a Monte Carlo simulation experiment. Notwithstanding the possibility to model the environment down to tiny details, the degree of approximation was arbitrarily chosen to provide a realistic amount of complexity in a reasonable amount of time. For a further discussion on the simulation detail-level, the interested reader should refer to, e.g., Zio (2009).
System components are distinguished by type, which defines the available actions when they interact with each other. The elements used to simulate the operations of the mine present a unique behavior, and they are divided in two macro-categories: the first is server-queue components, which include shovels, dumpsites or discharge points, and workshops; in this research, server-queue components are represented as in (Law et al. (2000), pp. 12-18). The second category is represented by agents, i.e. trucks for transportation of the excavated material around the mine. According to Law et al. (2000) an "agent is an autonomous "entity" that can sense its environment, including other agents, and use this information in making decisions. Agents have attributes and a set of basic if/then rules that determine their behaviors." The agents can travel between each couple of sites in the mine, and the traveling distance between sites is described by log-normal distributions. This choice allows us to sample the travel time from one site to the other in a realistic way.
The behavior of an agent, i.e. a truck, is characterized by the parameters of the processing time distributions described in the following. A truck is unreliable in the sense that it might fail at any moment during operation, and the time to failure (TTF) is a random variable modeled using a two-parameter Weibull distribution

Fig. 3
A schematic illustration of the adopted modeling approach Server Queue . 4 A schematic representation of a server-queue entity (shovels, dump sites, and maintenance workshops) where α > 0 is the shape parameter, β > 0 is the scale parameter, and t is the time elapsed since the last maintenance intervention. The time to repair (TTR) is different in case of corrective and preventive maintenance, and it is modeled using a log-normal distribution where the parameter μ is the mean of the distribution, σ is the standard deviation, and t is the duration of the maintenance intervention. Each truck is characterized by its capacity, which varies depending on the truck model, and a cost for both preventive C T P and corrective C T C maintenance interventions. From the practical point of view, the selected TTR-approach allows us to take advantage of the cumulated maintenance data as the peculiar characteristics of each piece of equipment can be represented.
An agent can be served by server-queue objects, i.e. by shovels, dumpsites, and workshops, which are modeled according to the well-know queueing theory (Law et al. 2000). A server-queue entity presents a waiting room, the so-called queue, where the agents, or customers, wait their moment to be served by the processor, the so-called server. Figure 4 gives a schematized representation of the serverqueue object, where the agents are represented by the circles and they join the queue at an unknown arrival rate. Customers are served according to a first in first out (FIFO) logic at a serving rate that changes depending on the type of customer.
The three classes of server-queue components present subtle differences. Shovels were modeled as server-queue objects with log-normal serving time distributions, and they presented the peculiar hallmark of unreliability: as in the real world, they were subject to the aging process, hence they could unexpectedly fail, or they could be preventively maintained. Therefore, in addition to the serving time, shovels are characterized by a TTF probability density function, which is modeled using Equation 1. When a shovel becomes unavailable due to maintenance, it changes its behavior to that of an agent and it enters the maintenance workshop with maximum priority. The trucks in the queue wait for the shovel to be available again and no other trucks are assigned to that shovel until maintenance ends. Shovels are thus characterized by a cost for corrective C S C and preventive C S P maintenance in addition to TTF and TTR distributions. As soon as the  Fig. 5 Scheme of the agent's movements on a map inside the maintenance module between its sub-systems. S i represents a generic shovel site, D i represents a dump site, and W i a workshop maintenance activity is completed, the shovel is considered available again and trucks can start to join the queue and to be processed. Workshops are characterized by a FIFO logic with priority for the management of the queue (shovels with maximum priority), and they present a peculiar behavior concerning the processing time of a customer, i.e. a truck or a shovel. The service time is a function of both the type of item served (truck/shovel) and the type of maintenance intervention, i.e. corrective or preventive. Finally, dumpsite components are characterized by a log-normal service time distribution and by the presence of a stockpile; each stockpile has a limited capacity and all the stockpiles feed a single concentrator plant with a specific capacity of material per unit of time. The detailed modeling of the concentrator plant, with equipment such as crushers, conveyor belts, mills, flotation tanks, etc., is left out of the scope of this paper and it is assumed to work without interruptions.
A DES experiment was designed to replicate system operations with a high level of detail. Within the simulation procedure, all the entities interact with each other as described and illustrated in Figure 5. The mine maintenance simulation is initialized by defining the parameters of the probability distributions; the TTF and TTR distribution parameters are listed, together with the costs for maintenance, in Table 3 and Table 4 in Appendix A. Trucks are also characterized by a transportation capacity, which is a random variable sampled from the distributions reported in Table 3 in Appendix A, whereas servers are characterized by a serving time distribution, which parameters are listed for shovels and dumpsites in Table 4 and Table 5 in Appendix A respectively. The parameters listed above remain, together with the duration of the simulation horizon, un-changed for all the runs of the experiment.
Once the simulation is initialized, a truck gets assigned to a target shovel S i (see Figure 5), thus it travels to the designed site and it joins S i 's queue. After being processed at a load site, a truck can leave the site due to two reasons: it can either fail unexpectedly and thus being sent to a workshop W i in Figure 5 for CM, or it can be sent to a dumpsite D i in Figure 5 for unloading. After the unloading, a heuristic decides if the agent must be preventively maintained, or if it can continue its regular operation. The decision to submit a truck to PM is based on the age of the truck, namely a TBM policy is adopted. If the threshold value p i for the i-th agent is lower than the time elapsed since the last maintenance intervention, it undergoes PM, otherwise, it is assigned to a new load site. The maintenance policy can be represented by a list, whose components p j i are the PM thresholds for trucks T and shovels S, and it can be represented as follow: where N T is the number of trucks, and N S is the number of shovels in the system.
When a CM or PM intervention is due on a truck, a workshop W i processes the agent according to the type of maintenance needed and to the TTR distribution of the specific agent. Once the truck has been maintained, its condition is considered "as good as new" from the modeling perspective and it is ready to start a new mission. A mission is defined as a chain of actions that includes the travel to a shovel site, the waiting time in queue, the loading and unloading operations.
The shovel's mission is less detailed than a truck's mission: each shovel simply operates at its site until a failure occurs, or until it is sent to a workshop W i for PM. When a truck has been loaded, the age of the shovel is checked against the age threshold p S i and, in case the time elapsed since the last CM/PM intervention exceeds p S i , it is sent to a workshop W i with maximum priority, thus preempting each other agent in the queue.
The performance of the system was optimized based on the results of a MCS experiment. Given the stochastic nature of a DES, the problem consists in the minimization of the expected value of the cost of operations J (θ ), and it can be formalized as where θ is a vector containing the system parameters that define the number of trucks N T and shovels N S , and the maintenance thresholds P. The problem must be solved under the constraint of reaching a production target X min : Pr{X ≥ X min } ≥ 0.95.
That is, the probability that the output X satisfies the target X min must be greater than 95%; such probability can be calculated using the 95 th percentile of the output distribution from the MCS experiment.
To minimize the objective function means to act on two aspects of the model: the number of resources operating in the system and the number of unplanned downtimes. The former is minimized using an enumerative search algorithm, while the second is optimized using a more complex genetic algorithm for search over a stochastic response surface. More details about both procedures are provided in Appendix B.
The code 1 used to implement the algorithms described above is written in Python 3.7 and mostly using SimPy simulation library.

Cash Flow Module
The use of system dynamics methodology allows building a compounded, close-to-reality representation of the mining operation that is still easy-to-read and modify compared to writing the model as software code. Detailed SD feasibility models of mining have been introduced in the literature by, e.g., Inthavongsa et al. (2016), Savolainen et al. (2017), who showed the flexibility of the approach and its ability to cope with complexity.
For the sake of brevity, the CF model used in this paper includes only two uncertainties: the metal price and maintenance cost. The simulation horizon is limited to one-year, and a geometric Brownian motion with and without trend is assumed to represent the uncertainty of markets adequately (for discussion see, e.g., Labys et al. (1999), Roberts (2009), Rossen (2015)). An example price simulation used in the experiments is illustrated in Figure 6 with three alternative price trend scenarios for a single random price realization. The uncertainty of maintenance costs are modeled as triangular distributions using expert estimates, which are introduced in more detail in Section 5.2.
A representation of the function block diagram of the applied CF model is provided in Figure 7. The model is divided into two sections: technical and economic models, where the inputs of the mine maintenance module are fed (blocks of the flowsheet marked with blue background). For a full list of parameters see Table 6 in Appendix C.
One of the key output variables of the CF model is the average mill utilization rate. That is, at any point in time, the mill utilization rate is either zero or one depending on the level of the ore stock that is replenished by the truckshovel system. We exclude the option to increase the size of temporary ore stock giving additional flexibility to maintenance timing, which is often used in small mines. In our case, the stock is limited to ≈ 27, 500 tons of ore, which corresponds roughly to 36 hours of production in the mill. The key added variables from the CF model include the costs of equipment leasing, fuel costs derived based on the O&M model's indicated operation hours, and other fixed costs such as buildings, and administration. The output price of metal is updated weekly.
We acknowledge that the above-described model construction, including the detailed operations & maintenance model using discrete event simulation and system dynamic cash flow modeling, could be fully implemented in a single software environment. In practice, this is usually not possible, which calls for the DT type of co-simulation approach. The reasons for this can be related to an unwillingness to share confidential financial information (from mine operator to the model owner), the effort of transferring existing pieces of core software libraries from one environment to another, and importantly, as pointed out by, e.g. West and Blackburn (2017), the uncertainty of financial return of the software product.

Model and Application
In this section, the O&M simulation optimization module's behavior is first validated with two sensitivity analyses and then used in concert with the SD model to run three experiments in a DT system setting. The mine configuration in these tests varies from one experiment to the other, and the number of components in the system is kept low to avoid over-parametrization of the model (see discussion, e.g., (Zio 2009).

Maintenance Model Validation
A simple sensitivity analysis was performed to validate that the lower the age threshold for PM, the higher the possibility to avoid unplanned downtimes. On the other hand, the higher the age thresholds, the less effective PM should be in reducing costs. To validate this hypothesis, the maintenance thresholds of all the items were parameterized as follows: where i identifies the item, j is the class T for trucks or S for shovels, and MT B F i is the mean time between failure of the i-th item. The parameter a ∈ (0, 3] is a scale factor that allows to vary the age threshold p j i of all the equipment included in the experiment in question. By parametrizing the age thresholds, it was possible to estimate both the cost of CM and PM for the whole system by changing only the parameter a. In all the other experiments, the cost of CM/PM depends also on C C , C P , and on the TTF distributions, but here these parameters are fixed. The sum of the cost of PM and CM at different values of a are plotted in Figure 8. The fleet used to realize the sensitivity analysis included two trucks and one shovel. When the age thresholds are very low ( MT B F) the cost of PM is high because PM events are carried out extremely often. However, the cost of PM decreases sharply when a increases and, with maintenance thresholds p j i equal to 0.5 times the MT B F values, the cost of CM starts to be higher than the cost of PM, thus making it inconvenient to perform PM more rarely. As it is depicted in Figure 8, the total cost of maintenance presents a minimum cost as a function of PM and CM, which makes clear the need to optimize the PM age threshold of all the items before running the whole simulation procedure. A second maintenance model validation experiment was carried out to test the performance of the system with varying maintenance resources; in particular, the difference between two-and three-workshops configurations were analyzed. A total of sixty configurations were tested, namely all the possible combinations of 2 or 3 workshops, 1 to 10 trucks, and 1 to 3 shovels. The statistics used to present the results are the average throughput and the average cost of maintenance obtained from 50 simulations over a 2-year time horizon. For each configuration, the maintenance thresholds p j i are optimized and then the DESs are run.
Since the dumpsites present limited capacity, i.e. material excavated cannot exceed the mill production rate, the configurations with two and three maintenance workshops produce different results. As shown in Figure 9, many system configurations deliver the maximum possible amount of material, but at different costs. Interestingly, highly different configurations lead to similar results: for instance, the 2-workshops 3-shovels and the 3-workshops 1-shovel configurations deliver almost the same throughput at the same cost using a similar number of trucks. The two solutions are however very different from a managerial point of view: the investments required to purchase or to rent the equipment, the skilled personnel needed to operate the facilities, and the resilience of the resulting system are meaningful aspects to be considered.
The above-mentioned issues go beyond the reasonable modeling scope of the DES, but these are the issues that can be easily integrated into the managerial profitability model to produce further insights to support operational decisionmaking. In a dynamic economic environment, provided by the SD, the proposed analysis can be repeated with better implicit knowledge of the production process, such as production targets and planned maintenance, thus producing a probabilistic evaluation of the future scenarios.

Digital Twin Testing and Validation
The first experiment aimed to verify and validate the overall DT approach, whereas the second experiment consisted of a more detailed optimization of the system under the assumption of uncertain maintenance costs. The parameters used in the CF model are illustrative, and they were chosen in a way that approximately 5-6 trucks (max. 10) with 1-2 shovels (max. 3) would satisfy the mill requirements for material tonnage.
The fleet design study was carried out to screen all the possible system configurations that can be produced by ten trucks, three shovels, and the maintenance policies provided in Table 2. In Table 2, the item 'design optimization' indicates if the number of trucks and shovels is set already in Fig. 9 A multi-criteria comparison of different system configurations. Each point corresponds to a system configuration (number of workshops, shovels, and trucks), and it is characterized by the average cost of maintenance and the average throughput of metal obtained over the simulation horizon the discrete event simulation. The configurations that have a smaller number of trucks than shovels were discarded manually as irrelevant, when the optimization option was turned off. Three different maintenance policies were used: a "max-corrective" (or "run-to-failure"), a balanced, and a "max-preventive" policy. The first and last policy represent the theoretical endpoints of the available scale of the simulation space: according to the "max-preventive", a PM action is performed after every mission, whereas according to the "max-corrective" the maintenance threshold is set (de-facto) to as infinite. The balanced policy foresaw one PM event per week of simulation. In a more advanced setting, the balanced policy could be defined by the simulation-optimization algorithm, which searches for the PM thresholds p j i that return the minimum expected cost of maintenance for a given configuration. In Experiment 1, maintenance costs were assumed to be fixed (known ex-ante), and the price trends were those displayed in Figure 6. The number of simulations displayed in the last row of Table 2 was determined by the number of combinations, e.g., in experiment one with 'pre-optimized' fleet design, the number of combinations to be simulated was nine as only the number of maintenance policies multiplied by the number of price trends.

Experiment 1 -Fleet design
Results of Experiment 1 without fleet optimization are provided in Figure 10 which shows that the policies "maxcorrective" and "balanced", with 6-10 trucks and 1-2 shovels, would be the most profitable ones. It is noticeable that with the given parameters of fixed costs, and duration of PM-and CM-events, the "max-corrective" policy was favored over the balanced option. The "max-preventive" policy produced negative profits in all cases within the selected set of parameters, which highlighted the need for further maintenance threshold optimization. That is, in this case, to maximize the amount of preventive maintenance leads to the lost of overall cost efficiency due to excess queuing times to the workshop, whereas increasing the number of workshops in the initial design would also have a bloating effect on operation's costs.
To inspect the results of Experiment 1 in more detail, the total number of maintenance events in the case of full PM are plotted in Figure 11. Figure 11a shows that the "maxpreventive" policy, with a two-workshops design, is possible only in the case of one truck and one shovel. As the number of trucks increases, the relative share of CM actions goes up since there is not enough capacity in the workshops ( Figure  11b) for adequate equipment intake.
Such effect is further highlighted in the case of one shovel and nine or ten trucks: the queuing time spent by trucks (either at loading, unloading, or maintenance) increases, thus making it more probable for them to fail before the next scheduled PM event.
In Experiment 1, the simulation-optimization algorithm in the O&M model was also tested to screen out the infeasible fleet designs already at the beginning of the simulation. The maintenance optimization, as designed, favored the high mill utilization rate options that were gained with eight to ten trucks and one or three shovels (Table 10).
It is clear that the results of maintenance optimization efforts are uncertain ex-ante even with the assumption of fixed maintenance costs. To take a further step, the role of cost uncertainty was included in the analysis by replacing the fixed maintenance costs with triangular probability distributions in the CF model. These distributions are depicted in Figure 12 using box-plot diagrams (for numerical values see Appendix C), and they represent expert knowledge. In a real case, these distributions could be derived from proprietary maintenance data that are available from the organization's historical records.
The applied cost distributions of PM and CM differ in shape. To reflect the risk of CM, the distributions of costs have long tails that can produce up to five times the fixed cost, whereas the positive risk of CM is limited to 1% respectively. The PM cost distribution is weighted more in the center of the distribution, thus giving less uncertain results on costs; it is capped to a maximum of 1.1 times the original assumption, and it can go below -10%. From the modeling point of view, this setting creates a strong incentive towards accepting preventive policies over the "max-corrective" given in Fig. 11 The limited capacity of workshops illustrated using the simulation statistics. Scenarios are taken from results with full preventive maintenance policy settings using an increasing-price array   the first experiment. Another option for the inclusion of cost uncertainty could be to include them in the O&M optimization module that is run before the CF simulation. However, this creates additional complexity to the genetic algorithm used for O&M-model's simulation optimization, and it complicates the user's abilities to interact with the CF simulation experiments within the SD-flowsheet that steers the CF model. A trade-off between model choice is made by using averages of distributions based on 10,000 draws in the optimization model (see Appendix B) and a single random draw in the CF model, which makes it more volatile in terms of results. Running Experiment 1 with uncertain costs returned the previously suggested outcome {9, 1, 98%} with no PM; this was due to the limited workshop capacity as previously discussed (see Figure 11). Therefore, the question of optimal maintenance policy boils down to finding out whether and what is the optimal time between maintenance events that would keep the amount of CM within reasonable limits.

Experiment 2 -Optimal Timing for Preventive Maintenance
The issue of optimizing preventive maintenance thresholds is addressed in this last experiment. The age threshold, marked as n, for the PM event timing was set as a ratio versus one round of simulation of the maintenance module. That is, a value of n = 1 means that there is approximately one preventive maintenance event per week of simulation in the O&M model, and n = 0.5 indicates one PM event every two simulated weeks respectively. In this experiment, n varies from 0.125 to 1.000 with a step size of 0.125.
The visual insights are provided in Figure 13, which shows that the {9, 1, 91.6%} combination with n = 1 yields a profit of approximately 30 million. On the other hand, the simulation indicates another option with two trucks less, namely {7, 1, 81.7%} with n = 0.875, that has only some 0.5 million profit less than the "best option". From the decisionmaking point of view, we can observe that the interpretation and efficient utilization of simulation results becomes difficult because of the high number of co-existing solutions. In this vein, as previously discussed by, e.g., Negahban and Smith (2014), Min et al. (2019), there is an increasing demand to implement meta-model based solutions to simplify simulation data of DT-models into implementable managerial insights. This issue is, however, beyond the scope of this paper.
As depicted in Figure 13, maintenance policies with n ≥ 0.5 seem to deliver the best performances in terms of profitability and utilization and, furthermore, a positive relationship between the mill utilization and a high num- ber of trucks exists. On the other hand, based on Figure  13, a reduced number of shovels accounted for an increase of profitability. To verify when adding further equipment is no more beneficial from the profits point of view, Figure  14 is laid out, which depicts the "the profit per truck" as a function of the number of trucks in the system with varying maintenance policy thresholds and shovel counts. Figure 14 suggests that selecting a reduced fleet size could increase the overall profitability of operations (see the middle subplot). The unprofitable configurations, i.e. those showing an unbalanced amount of trucks and shovels, locate themselves into the lower-left quadrants of the plots, where profits are negative. Moreover, Figure 14 highlights that, when investing in three shovels, one should prefer maintenance policies with n ≈ 1.0, whereas a smaller fleet size could perform at its highest with a n ∈ [0.8, 0.9].

Result Summary and Analysis
The limited set of results shows that in the mining industry the co-simulation of detailed O&M models for fleet operations in connection with the managerial CF models may provide a way to improve the overall profitability. The potential for economic improvement in this study was brought by the possibility to combine fleet design and maintenance policy, which, based on the detailed simulation, could be able to produce a comparable high revenue per invested unit of money without simply striving to maximize either the rate of excavation or mill utilization that are usually viewed as the most important key performance indicators. In the model demonstration, a total time frame of one year with weekly simulations in the O&M model was used. For real-world applications with constantly updating operational data, a more frequent optimization might be a more interesting choice, where the managerial considerations could involve the choice of alternative mine plans with less focus on the current market price.

Discussion and Conclusions
This paper presented a DT modeling approach aimed to dynamically optimize the O&M in the mining industry with respect to the uncertain price of the end-product. A twostage simulation was adopted: firstly, a selected set of mining equipment with their unique failure distributions was chosen, and its capabilities to meet the production demand in the concentrator plant was simulated using discrete event simulation. In the second phase, the aggregated information provided by the DES was fed to the managerial CF model with costs fully accounted to evaluate the overall economic optimum from the operations perspective.
The DES experiment produced useful information on the availability and utilization of mining equipment, and we were able to connect the aggregated output data with the overall profitability of the business via the managerial CF model to form a "digital twin of manufacturing" in mining. Within the limits of the selected simulation values, it seems that in most cases maximizing mill utilization and production throughput would yield the highest expected revenues regardless of the realized price array, but some opportunities for downsizing may exist. The results highlight the role of maintenance as a necessary evil with only little potential to the economic upside. However, these results are highly dependent on the selected values and should be re-evaluated using a more extensive sensitivity analysis of the DT presented or a case example of a real mine.
Some implications to managerial decision making of the mining industry can be suggested. First, simulation-based DT modeling can perform overall operational optimization while considering the stochasticity and high dimensionality of operational data, which, as the second point, allows simultaneous consideration of operational and investment decisions. From the methodological point of view, the whole concept of DT is still taking its shape, and often it is used just to bring extra buzz to promote one's simulation efforts. We used the term specifically to refer to a holistic system optimization model that consists of several, data-based, selfreliant, and discipline-specific, simulation models. These models run parallelly and they are connected in the virtual space, which could have a (near) real-time connection with the physical process. Whether in the future the approach would be labeled DT or not, the models evidently have novelty value beyond simpler simulations applied today. Insights gained from our modeling efforts corroborate the earlier findings and discussion around the problems in the implementation of DTs. These issues include, but are not limited to, i) the dimensionality and the degrees of freedom of data in large models, ii) a suitable level of data aggregation when transferring it between discipline models, and iii) the validation and verification (V&V) of the results. A central limitation related to the V&V of this paper is the absence of a true physical mining system; on the other hand, fully virtual modeling enabled testing and running experiments rapidly for scientific research. It should be also acknowledged that the TTF distributions are not able to fully capture the reality of maintenance: additional details could be added in the O&M model by making the failure probabilities dynamic or component-based to reflect the agents' actions better than the static ones used here. As future avenues of research, a case-specifically tailored DES-model could be taken into action, which is connected to the relevant datagathering systems for up-to-date information as well as the existing CF models. On the other hand, rather than increasing the details of modeling, there is an emerging need to develop meta-model based methods that can interpret the results of DT-simulation in a fast and managerially digestible manner.
Funding Open Access funding provided by LUT University (previously Lappeenranta University of Technology (LUT)).

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Experiment Parameters
The parameters used to characterize the components used in the simulation experiments are reported in the following. The TTF of trucks and shovels are the same used in (Mena et al. 2013). The letters μ and σ represent the mean and the standard deviation of a Gaussian and a log-normal distribution. C T C and C T P are respectively the cost of corrective and preventive maintenance of a truck The letters μ and σ represent the mean and the standard deviation of a log-normal distribution. C S C and C S P are respectively the cost of corrective and preventive maintenance of a truck

B Simulation Optimization Details
To optimize the configuration of the system means to identify the number of trucks and shovels that guarantees to achieve the production target X at the minimum cost. The response surface E[J (θ )] is obtained through an MCS experiment and the effect of constrained resources -i.e. the number of maintenance workshops and the maximum capacity of stockpiles -is not easily predictable. We observed that the number of available workshops makes the increase in throughput non-monotonic with respect to an increase in the number of operating trucks and shovels. Therefore, the possibility to guide the search according to some heuristic must be dropped, and given the relatively low number of possible combinations, an enumerative search procedure was implemented. Once the average values of throughput produced by the system were available, the solutions were ranked in ascending order concerning the expected cost of maintenance E[J (θ )].
Maintenance thresholds are continuous real variables, which are constrained to be positive. Since there are no Legend: letters μ and σ represent the mean and the standard deviation of a log-normal distribution other constraints, an effective tool to find a minimum of the response surface E[J (θ )] is the genetic algorithm (GA). Each individual of the population was represented by a vector P, which was introduced in Eq. (3), and the size of the initial population was set to 50 individuals. The selection of parents occurred according to the "fitness proportionate selection" mechanism, whereas single-point crossover operations and mutation operations were performed with a preference for the latter (a ratio of 4 mutations to 1 crossover was used). Crossover operations consisted in the selection of a random element of two vectors P, at which the vectors are separated in a head and a tail, and subsequently the tails are swapped. Mutation operations occurred with a probability of 0.2 for each element of P, and when the mutation occurred a quantity sampled from a normal distribution ∼ N (0, 0.1) was added to the element. The stopping criterion used in the GA was the limit on the number of generations, which was set equal to 25. The GA was indeed tested using a different number of generations up to 100 and the algorithm showed no meaningful improvement of the best solution found after 30 generations. The number of DESs required to obtain a reliable estimate of the response surface value returned by an individual was estimated in 50 repetitions. Such requirement makes the execution of the algorithm computationally demanding, therefore a 64-bit Windows Server 2016, Intel Xeon Platinum 8160 CPU 2.10 GHz, with 768 GB of RAM workstation was used to run the algorithm; a single optimization required on average 30 min.