A digital twins concept model for integrated maintenance: a case study for crane operation

Abstract

The paper presents an Integrated Maintenance Decision Making Model (IMDMM) concept for cranes under operation especially into the container type terminals. The target is to improve cranes operational efficiency through minimizing the risk of the Gantry Cranes Inefficiency (GCI) results based on the implementation of the Digital Twins concept for maintenance purposes. The proposed model makes a joint transportation process and crane maintenance scheduling, relevant to assure more robust performances in stochastic environments, as well as to assess and optimize performances at different levels, from components and transport device to production systems (container terminal). The crane operation risk is estimated with a sequential Markov chain Monte Carlo simulation model and the optimization model behind of IMDMM is supported through the Particle Swarm Optimization algorithms because the objective function a non-linear stochastics problem with bounded constrains. The developed model allows the container terminal operators (management process) to obtain a maintenance schedule that minimizes the GCI (holistic indicator), as well as establishing the desired level of risk. The paper demonstrates the effectiveness of the proposed maintenance decision making concept model for cranes under operation using data from of a real container terminal (case study).

Introduction

In the last two decades, container transportation system has been faced under increasing development and it is predicted that this increase will have a rate of about 10% until 2020 according to Widarto and Handani (2015). This fact shows the importance of container transportation system as a key role of container terminals to link between sea and land. Although container terminals are growing their capacity to respond to these current demands, the rapid increase in the transportation of containerized goods has created a continuous need for the optimal use of equipment and the facilities in the port. The current target is a continuous effort to decease the operational costs in order to improve the performance of the ports.

In the container terminal operating system, gantry cranes are critical transport devices and major bottleneck restricting the working efficiency of the harbour. Container terminal managers pay more and more attentions to improve operational efficiency of gantry cranes (Ri et al. 2012; Euchi et al. 2016; Cheng-Ji et al. 2014; Dadashi et al. 2017; Seung et al. 2012; Abourraja et al. 2017; Yan and Kap 2011) and the maintenance process is one of fundamental aspects to ensure its proper functioning (Smoczek and Szpytko 2012, 2017; Liu et al. 2018; Voisin et al. 2010). The increasing of loading/unloading activities requires the proper readiness of supported infrastructure, including gantry cranes, and concepts like DT (Qi and Tao 2018; Tao et al. 2018), Closed Loop Engineering (Barari and Pop-Iliev 2009; Jie and Nan 2017), e-Maintenance (Voisin et al. 2010) are the perfect tools to improve the design, production planning, manufacturing, and maintenance in this system.

Of particular interest to this research is the relationship between maintenance process and DT framework. Almost all of today’s maintenance-oriented applications touch in some way the DT framework based on existing IT (Information-Technology) tools. Reference like (Melesse et al. 2020) explores between 2016 and 2019 papers related with DT keyword, and the second most used application domains is oriented to maintenance process. Table 1 shows a collection of papers whose main contribution, on maintenance process domains, is supported by DT frameworks and at the same time illustrates guidelines reasons of the selected approach applied in this research. The DT has been defined for more than 15 years and in the last few years the concept has been shaped, leveraging on the exponential increase in the power of the IT tools, sensors, embedded systems, IoT, able to generate and analyze plenty of data. Grieves (Grieves 2014) is considered to be the pioneer of this concept and his postulates about the DT dimensions (physical entity, virtual model, and connection) are used in almost all the papers related to the framework to guide the contribution and the proposed methodology, as shown in Table 1.

Table 1 DT maintenance-oriented approach

The papers referenced in the Table 1 are clear DT applications on maintenance process and in all cases the results achieved are accurate, but in some cases, it is not easy to distinguish the differences between the dimensions, however, the reference (Tao et al. 2018) provides the desired guidelines for the study conducted in this research: maintenance process in a container terminal (critical interaction: gantry cranes-vessel arriving frequency planned). Therefore, as an accepted definition, we adopt the one coming from (Tao et al. 2018) “A dynamic model in the virtual world that is fully consistent with its corresponding physical entity in the real world and can simulate its physical counterpart’s characteristics, behavior, life and performance in a time fashion”, as well as the five dimension framework formed by Grieves postulates (physical object, virtual counterpart, connection, data and services) introduced in the same work.

Paper contribution

Following the same definition and framework, but adapting to the requirements of this research, Fig. 1 illustrates the schema concept proposed to improve maintenance scheduling support process based on DT concepts and consequently the operational efficiency of the container terminal system considered, and Fig. 2 shows the five-dimension DT architecture adopted, which can be summed in the expression IMDMM_DT = f (Physical Object, Virtual Counterpart, Connection, Data, Services).

Fig. 1

DT concept for scheduling maintenance purposes

Fig. 2

Five-dimension DT model for maintenance management

The following sections describe each dimension according to the analyzed system and the last section adds additional remarks related to the implementation in practice.

Physical object

The current industry has higher automation levels and complexity (container terminal is only one of them), therefore, decompound the Lego system in critical pieces simplify the problem and gives us the opportunity to focus in specific process. Degradation process is inherent on technical systems; therefore, maintenance is pertinent everyday more and the human decision-making process behind is a target to be improved.

In this paper, the physical object (PO) studied is the maintenance scheduling management process in a container terminal composed by 10 gantry cranes, two power transformers, two emergency generators and two transmissions lines connecting with a Power System (see Fig. 2). The electrical Power System together with generator support system guarantee a stable operation in the terminal containers, but the relation between the frequency of arriving vessels and gantry cranes movement capacity in the container terminal is the critical relation studied (see Fig. 3).

Fig. 3

Cyber-physical system diagram

Virtual counterpart

The virtual counterpart (VCo) is an IMDMM based on the PO computational modelling, in fact, is a DT tool that simulate the same process performed in the time-real, but in this case, an interactive (because the model consider freedom variables, planned variables) optimization algorithm chooses the best maintenance scheduling given the vessel frequency planned in the coming year and the historical data from the previous process as a result of machine learning analysis, and then gives the feedback to the container terminal entity manager (planner) as a closed-loop engineering system. In our case, VCo depends on two main variables: system capacity and system load, VCo = f (VC, TC) in our model. The implementation of VCo or IMDMM in practice is aligned with the field of robots. References like (Neil et al. 2008; Foumani et al. 2015, 2018; Baniasadi et al. 2020) are examples of similar studies in different environments. They all describe the same function studied in this research, IT tools, robots in this case, that help managers to face scheduling challenges, it means, once the sequence of the process to be conducted is well known (PO), software and robots together (VCo) perform all the tasks. Depending on the real-world application, restrictions are needed. For instance, typical restriction to be considered in cooperative gantry crane system is the safety distance between them, but, in our case, we assume a well-designed system (each gantry crane works in a well-defined set of space), reason why, the restriction is not considered and as a consequence the system capacity (TC variable in the model) is a parallel composition of gantry cranes.

Figure 4 shows the conceptual relation between VC, and TC. In Fig. 4, two main process are conducted: operation and maintenance. Knowing the operating system capacity and planned system load in the container terminal (modelling and simulation processes), allow us to assess the convolution between them and measure the loss capacity, defined forward as (GCI). Consequently, we can also evaluate the maintenance scheduling impact with the same indicator GCI in a closed-loop engineering system using an optimization model because the loading capacity of the gantry cranes consider both process (maintenance and degradation), VC = f (LC^D & LC^M) in our model. Therefore, the maintenance strategies (maintenance as a support process of the operation) change depending on operation variables (see red arrows in Fig. 4).

Fig. 4

Maintenance management process

Connection

The formal connection (C) between the PO and the VCo is a database structure created for maintenance purposes (see Fig. 2). Basically, the database stores the collected, filtered, analyzed data by unique component ID considered. This intermediate structure guarantees the practical implementation of the optimization model because allow us to merge the historical data and the planned process, C = f (Planned Parameters Set, Historical Parameters Set), in our case as follows:

• 1. Planned Parameters Set = {λ_VF, VF_n, VN, VC}

• 2. Historical Parameters Set = { $\overline{{LC_{i,j} }}$, µ_GCi, σ_GCi, TTF_i,s, TTR_i,s, TTM_i,k, TDM_i,k}

Data

The connection structure is linked with data, our case the data used is historical degradation information of all the components in the system, previous planned process, system structure, and so on, collected by SCADA (Supervisory Control And Data Acquisition) and SAP (Systems, Applications & Products in Data Processing) systems, defined as D = f (SCADA system, SAP system), in fact, historical information available related with the maintenance scheduling management process:

• 1. Nominal parameters of the components involved.

• 2. Historical degradation, time to failure (TTF), time to repair (TTR).

• 3. Maintenance sequences according with operation times, life cycle, manufacturer recommendations.

• 4. Planned vessel frequency.

• 5. System structure.

• 6. Components localization.

• 7. Components efficiency.

Services

The DT tool finality in our case is give services to the entity manager of the container terminal ensuring fasters and optimal decision-making process. Once the maintenance scheduling management process is optimal based on the model freedom variables, consequently, the holistic operational efficiency in the container terminal increase. In our case, S = f (GCI indicator, Planned Parameters Set).

Additional remarks

Generator support systems, including power transformers and transmission lines, are examples of complex system which have a dynamic characteristic time by time according to (Szpytko 2004). Stochastic type methods can simulate the dynamic behavior of a system time by time and in some conditions and assumptions. In our case, the physical object studied include the generator support system in the maintenance scheduling process because the coordination is crucial to guarantee a proper container terminal functionality (reliability standards).

The paper is focusing on developing a joint transportation process and crane maintenance scheduling relevant to assure more robust performances in stochastic environments, as well as to assess and optimize performances at different levels, from components and transport device to production systems.

The human decision-making process behind of the maintenance scheduling management process is the coordination of components maintenance and/or replacement that make up the system but maintaining holistic and/or clustering objectives defined by the decision makers. Coordination itself can be an complex problem and humanly dreadful to find a faster optimal solution, mathematically speaking is a nondeterministic polynomial time NP-complete problem, therefore, listing all possible operational scenarios to make a coherent coordination in order to find the best scenario is difficult, reason why, computational modeling supporting by DT concepts are pertinent in this case. In addition, maintenance scheduling is an open problem because engineering systems increase their complexity and variety every single day. Below, we discuss the selected approach used to implement and solve the current problem identified briefly justifying the reasons.

Scheduling as a general problem, can be decomposed essentially by hierarchical levels: holistic objective (referring to global or integrated strategies) or multi-objectives (referring to decentralized strategies), and optimization criteria: cost, reliability, or hybrid approach. Examples are cost-holistic (Zhi et al. 2020), cost-multi-objectives sequential (Briskorn and Zey 2019) and holistic-reliability approach (Luo et al. 2019). For any of the cases, the problem is to find the best scheduled sequence of actions for each component considered in the system. Generally, the objectives and restrictions are not well defined because depends on the individual system requirements (available information, relevant goals, application in practice). However, as a consensus, the optimization problem is defined as a multi-criteria combinatorial problem of non-linear objective functions with constrains (adding in our case stochastics), and the problem objective is to determine the timing and sequence of the tasks periods of each component analyzed. Therefore, the variables x in a scheduling problem is represented by the start time of the tasks for all the component considered. Referring to the maintenance, the conclusion is the same, only in this case, the task to coordinate are maintenance tasks.

The reference (Briskorn and Zey 2019) introduces in the literature overview, the contributions of different authors since 2009, almost all the contributions, including the reference (Briskorn and Zey 2019), focus on the interaction between cranes and do not consider the interaction frequency arriving vessel-gantry crane, however, all the proposals are viable practical applications regardless of strategies to solve the problem. On the other hand, the authors (Castilla-Rodríguez et al. 2020) in the literature review highlight the potentialities of the simulation approach to fill the gaps in the well-established historical models proposed since 1989, but all applications are oriented to the operation, not maintenance oriented.

Reference (Liu et al. 2018) is a maintenance-oriented contribution that highlights how maintenance activities have been ignored in the literature on terminal container systems. The contribution is aligned with the operation and considers the interaction between vessel capacity-gantry cranes because minimize the turn-around time of the vessel, but the interaction is fixed, it does not analyze how the vessel arriving frequency can condition the maintenance schedule. Furthermore, the proposal is mathematically complex and is only applicable in specific scenarios with necessary simplifications. The model proposed in this paper attempts to find the missing interactions suggesting a simple model mathematical formulation with a minimum set of input parameters.

Depending on the approach selected to resolve the maintenance scheduling problem, the volume of initial information required, the number of models constrains, and the complex mathematical formulations are issues to consider. These three dimensions above sometimes define the approach to be used. Multi-objectives approach introduces constrains and modeling challenges but is closer to the system needs. Holistic-objective approach sometimes do not represent all the system requirements. Cost-objective approach required high information levels and assumptions to standardize the selection criteria. Reliability-objective approach is not well accepted and understood by the industry. In conclusion, the approach selection depends largely on the characteristics of the system under analysis and implementation requirements.

For us, maintenance scheduling solution implemented is a holistic-reliability approach. Holistic, for an easy understanding of the maintenance impact on the system, understood as decision-making oriented to analyze the process as a whole (unique indicator) by the entity manager, and reliability, because a container terminal is an isolated and critical system working in continuous process, by construction, the system must be reliable. In addition, the selected approach requires generally available information and the probabilistic mathematical formulation is translated into technical terms by the GCI indicator. Therefore, the approach selection is guided by individual system requirements.

Especially, the paper is focus on Integrated Maintenance Decision Making Model (IMDMM) concept definition for cranes under operation into the container type terminals, and in a specific scenario description as an application example. The DT type model target is to improve cranes operational efficiency through minimizing the risk of the Gantry Cranes Inefficiency (GCI) indicator as a result of implement the maintenance scheduling process described in Fig. 4.

The DT definition and description is extensive and has several layers. The current paper assumes a certain DT connection and resolves a potential scenario with the virtual counterpart.

The document is organized as follows; firstly, the mathematical formulation of the optimization problem is presented with the constraint set and the flow diagram, as well as all the equations and assumptions of the proposed model organized in subsections. The optimization problem formulation is the virtual counterpart shadow, and the flow diagram describe the connection between the entity manager (services) and maintenance scheduling management process (physical object). Following the same idea, in next sections the scenario (parametrization and calibration) and methods (solution) are described, discussed and validated with a potential real system linking with the database structure for maintenance purposes. Further several comparative and possible scenarios are presented with a real container terminal (case study). Finally, some important conclusions are drawn to highlight potential outcomes in this research area and the connection with another investigation pieces.

IMDMM: Mathematical modeling

The proposed model objective is to minimize the expected value of the convolution function defined as E[R] (GCI indicator), between the gantry cranes capacity distribution function of the container terminal affected for the generator support system influence (container terminal capacity) defined as X and the containers number distribution function of the arriving vessels (container terminal load) defined as Y. The model is defined below:

$$ \begin{aligned} & \hbox{min} \quad E[R] = \frac{1}{Ns}\sum\limits_{n = 1}^{Ns} {r_{n} } \hfill \\ & \quad = \frac{1}{Ns}\sum\limits_{n = 1}^{Ns} {\left( {\left\{ {\begin{array}{*{20}c} {\frac{{\sum\limits_{t = 1}^{T} {Y_{t} - X_{t} } }}{{\sum\limits_{t = 1}^{T} {Y_{t} } }} \times 100} & {\text{if}} & {X_{t} < Y_{t} } \\ 0 & {\text{if}} & {X_{t} \ge Y_{t} } \\ \end{array} } \right.} \right)_{n} } \hfill \\ &\quad {\text{subject}}\;{\text{to:}}\quad 0 < x_{i} < T - (\sum\limits_{k = 1}^{{NM_{i} }} {TTM_{i,k} } + \sum\limits_{k = 1}^{{NM_{i} }} {TDM_{i,k} } ) \hfill \\ &\quad \;\;X,\;Y \in {\mathbb{R}}^{ + } \to R \in {\mathbb{R}}^{ + } \hfill \\ \end{aligned} $$

(1)

where

• E[R] is the expected value of the risk function (GCI).

• Ns number of simulations, R risk function, R = {r₁, r₂, …, r_Ns} values of risk, is the convolution product of the distribution functions X and Y. Therefore, R → f (X, Y).

• X gantry cranes capacity distribution function of the container terminal affected for the generator support system influence, depends on the total loading capacity LC_i of the container terminal and other variables. Therefore X → f (LC_i, …).

• LC_i loading capacity of the crane i, depends on TTM_i,k (start time to maintenance), TDM_i,k (time duration maintenance) and other variables. Therefore LC_i → f (TTM_i,k, TDM_i,k, …).

• x_i independent variable of the objective function, x_i = TTM_i,1, the start times for the first maintenance of each gantry crane or component considered in the system.

• T is the simulation window.

• N_C number of components considered in the system.

• NM_i is the number of maintenances of each gantry crane or component considered in the system.

• i = 1, 2, …, N_C.

• t = 1, 2, …, T.

• k = 1, 2, …, NM_i.

The stochastic non-linear optimization model with bounded constrains proposed for the gantry cranes, generator support system, power transformer and transmission lines maintenance scheduling problem solution in the container terminal present only continuous variables x_i = TTM_i,1 and is defined in the model constraint intervals (user definition, usually the intervals are defined by the container terminal manager in the database structure, periods depending on financial and/or human resources). The independent variable of the objective function to be optimized x_i = x₁, x₂, …, x_NMi depends on the quantity of maintenance tasks NM_i to be coordinated for each gantry cranes, generator, power transformer and transmission line (database structure stores cycle life information for each devices).

The optimization variables are only the start times for the first maintenance of each unit TTM_i,1. Once TTM_i,1 is established, the remaining TTM_i,k where k ≠ 1, are calculated adding the corresponding maintenance intervals, which are invariable and depends on the operation time between two consecutive maintenance tasks (SCADA system monitors devices operation times and SAP system establishes cycle life information, therefore, the link between both process converge in the database structure defining the maintenance cycle for each devices considered).

Figure 5 shows the computational algorithm flow diagram implemented to solve the online maintenance scheduling given the scenario inputs. The input parameters and the model variables are defined at the beginning of the paper. The computational algorithm flow is linear and only has two conditioning moments, first one to guarantee the simulations error criterion, and second one to guarantee the best solution of all the IMDMM proposals evaluated.

Fig. 5

Flow optimization model

The problem is solved as follows: the optimization algorithm proposes a set of TTM_i,1 and N_S Monte Carlo simulations are performed to determine N_S values of risk (r₁, r₂, …, r_Ns). The Risk mean E[R] and variance V[R] are determined from (r₁, r₂, …, r_Ns) and the error criterion is checked. If the desired error is not achieved, N_S is augmented and the Monte Carlo simulations are repeated for the same set of TTM_i,1. When the desired error is achieved the process is repeated for another set of TTM_i,1. This is done several times (1, 2, 3, …, N) determining decreasing values of Risk mean E[R₁], E[R₂], …, E[R_N]. The set of TTM_i,1 leading to the lowest Risk mean (E[R_N]) is the solution. The output model is the maintenance scheduling for all the devices considered in the system. By construction, the optimization model variables are associated with the maintenance process, but other stochastics variables considered (resulting from machine learning analysis of historical degradation data) in the model to estimate the risk are described in the next sections.

The optimization model construction formalized in this section is structured in two steps. First one, the vessels demand, gantry cranes, generator support system, power transformers and transmission lines stochastic mathematical model for a container terminal are defined. The Markov model used for the generator support system has four reliability states, considers random faults intrinsic to these systems, included the faults repair times and the probability of starting failure for the generators. The Markov model used for the gantry cranes, power transformers and transmission lines have two reliability states. The vessels demand considers the frequency stochastic behaviour hourly with which they arrive at the port and the vessels container capacity. Second one, the Monte Carlo Markov Chain (MCMC) model used to estimate the GCI risk indicator and the block diagram structure used to relate gantry cranes/generator support system is formalized.

Vessels demand modelling

The period of time between the arrival of two consecutive vessels in this paper is considered that follows an exponential distribution rounded to the nearest integer with parameter λ_VF, which we will denote in this investigation as VF_n ~ E(λ_VF), where n = 1, 2, …, VN—number of vessels that arrive at the container terminal in the simulation window T. The n-th independent random numbers generated from the distribution function guarantees the follows restriction $\sum\limits_{n = 1}^{VN} {VF_{n} } \ge T$.

One of the vessel features is its length and capacity. Each vessel carries several containers to the port for unloading, and each vessel loads a specific number of containers and leaves the port. The containers number $\overline{{VC_{n} }}$ of the n-th vessel is chosen according to empirical distribution shown in Table 2, generating u uniformly distributed random numbers [0, 1].

Table 2 Empirical distribution

The compound distribution function that models the containers number behavior VC of the n-th vessel that arrives at the port, at the time instant t is defined in Eq. (2) below:

$$ VC = f\left( {t,\theta } \right) = \left( {\left( {\overline{{VC_{n,t} }} } \right)_{t = 1}^{{VF_{n} }} } \right)_{n = 1}^{VN} $$

(2)

where θ is a set of parameters, θ = { $\overline{{VC_{n} }}$, VF_n, VN}, t = 1, 2, …, T (simulation window) and n = 1, 2, …, VN ensuring the restriction $\sum\limits_{n = 1}^{VN} {VF_{n} } \ge T$.

Vessel demand or vessel frequency is a crucial variable for the container terminal manager, in fact, is the freedom planned variable in the model (see Fig. 1) and defines the starting point in the DT model proposed. Given the planned vessel frequency (DT input), the model return the maintenance scheduling and the GCI indicator measured (DT output).

Gantry cranes modelling

The gantry cranes operation is continuous, eventually fails and is repairable. This random behavior can be described from Markov processes (Soszynska 2012). Considering the operation effectiveness of the container gantry crane, in this paper, we fix that the system and its components have two states z = 0, 1 and between them transition rates are defined. If the probability function of the transition rates from one state to another is exponential, as an example, they are denoted as λ (failure rate) and μ (repair rate) of the gantry crane, but as we know, the probability distribution can be change depending of the analyzed historical degradation data. Figure 6 shows a Markov process with two states: available and unavailable, and its transition rates λ and μ when exponential is used. The probability of moving from one state to another depends on the failure or repair rate of each gantry crane. For the two-state system represented in Fig. 6, the system of differential equations with initial conditions P₀(t) + P₁(t) = 1, P₀(0) = 1 and P₁(0) = 0 that models the Markov process is shown in Eq. (3).

Fig. 6

Two-state model for a gantry crane

$$ \begin{aligned} \frac{{dP_{0} }}{dt} = - \lambda P_{0} \left( t \right) + \mu P_{1} \left( t \right) \hfill \\ \frac{{dP_{1} }}{dt} = - \mu P_{1} \left( t \right) + \lambda P_{0} \left( t \right) \hfill \\ \end{aligned} $$

(3)

Equations (4) and (5) are the stationary solution of the system of differential Eqs. (3).

$$ P_{0} = \frac{\mu }{\lambda + \mu } = A $$

(4)

$$ P_{1} = \frac{\lambda }{\lambda + \mu } = U $$

(5)

In the two-state model, the gantry cranes are considered fully available (z = 1) or totally unavailable (z = 0). According to the standard systems, each gantry cranes should carry out 25 moves per hour which is equal to 144 s for every movement (Azimi and Ghanbari 2011). In order to simulate a real operation, the j-th number of movements $\overline{{LC_{i,j} }}$ ~ N($\mu_{{GC_{i} }}$,$\sigma_{{GC_{i} }}$), where $\mu_{{GC_{i} }}$ is the average move/hour and $\sigma_{{GC_{i} }}$ the standard deviation assumed. The stochastic capacity $LC_{i,j}$ at the time instant t of a gantry cranes i is determined by the TTF_i,s, TTR_i,s and $\overline{{LC_{i,j} }}$. The loading capacity LC is a result coming from gantry cranes standards, but the variables TTF_i,s and TTR_i,s are results of historical degradation data monitoring by SCADA system and consequently, filtered and organized in the database for maintenance purposes (DT intermediary connection). The maintenance scheduling is a planned process, therefore, in order to simulate failures in the system for risk calculation purposes (prediction), theoretical distribution functions are fitted to the historical degradation data of each gantry crane using machine learning (maximum likelihood estimation MLE framework), then generating TTF and TTR random samples from the fitted distributions is possible simulate the chain of TTF_i,s and TTR_i,s. The parameters TTF_i,s, TTR_i,s and $\overline{{LC_{i,j} }}$ allow to simulate with (6) the behavior of $LC_{i,j}$ generating s-th independent random numbers. Based on the previous experience (Azimi and Ghanbari 2011), we assume TTF_i,s ~ W(α_i, β_i) and TTR_i,s ~ N(μ_i, σ_i), where α_i and β_i are the shape and scale parameters of the Weibull distribution function respectively, µ_i is the average time repair and σ_i the standard deviation assumed of the Normal distribution function respectively for each i-th gantry crane. The model proposed to simulate gantry cranes capacity is defined below:

$$ \begin{aligned} & LC_{i,j}^{D} = f(t\left| \theta \right.) \\ &\quad = \left\{ {\begin{array}{*{20}l} {\overline{{LC_{i,j} }} } & {\text{if}} & {t < \sum\limits_{s = 1}^{m} {TTF_{i,s} } + \sum\limits_{s = 1}^{m - 1} {TTR_{i,s} } } \\ 0 & {\text{if}} & {\sum\limits_{s = 1}^{m} {TTF_{i,s} } + \sum\limits_{s = 1}^{m - 1} {TTR_{i,s} } \le t < \sum\limits_{s = 1}^{m} {TTF_{i,s} } + \sum\limits_{s = 1}^{m} {TTR_{i,s} } } \\ \end{array} } \right. \end{aligned} $$

(6)

where θ is a set of parameters, θ = { $\overline{{LC_{i,j} }}$, TTF_i,s, TTR_i,s}; i = 1, 2, …, N_GC; j = 1, 2, …, T (simulation window); s = 1, 2, …, NK_i and m = 1, 2, 3, …, NK_i knowing that NK_i depends on the simulation window used in the optimization model. The s-th independent random numbers generated from distribution functions guarantees the follows restriction $\sum\limits_{s = 1}^{m} {TTF_{i,s} } + \sum\limits_{s = 1}^{m} {TTR_{i,s} } \ge {\text{simulation window}}$.

The container terminal capacity compound distribution function TC defined in Eq. (7) is determined by the j-th gantry cranes total capacity in the container terminal (all the gantry cranes are in parallel) and the time of each vessel in the container terminal VF_n:

$$ TC = f(t\left| \theta \right.) = \left( {\left( {\sum\limits_{j = 1}^{{T_{{}} }} {\sum\limits_{i = 1}^{{N_{GC} }} {LC_{i,j,t,n} } } } \right)_{t = 1}^{{VF_{n} }} } \right)_{n = 1}^{VN} $$

(7)

where θ is a set of parameters, θ = {LC_i,j, N_GC, T, VF_n, VN}, i = 1, 2, …, N_GC, t = 1, 2, …, T (simulation window), j = 1, 2, …, T and n = 1, 2, …, VN ensuring the restriction $\sum\limits_{n = 1}^{VN} {VF_{n} } \ge T$. where n = 1, 2, …, VN.

On the other hand, one of the factors that affects the capacity of container terminals, is not stochastic and is not considered a random phenomenon, it is type maintenance strategy used into gantry cranes. The maintenance is contemplated within the strategies of a container terminal because it guarantees planned work time for the cranes. Maintenance is the activity designed to prevent failures in the production process and in this way reduce the risks of unexpected stops due to system failures (see Fig. 4). In a container terminal, to perform some maintenance tasks it is necessary that the gantry crane does not work, and this causes loss of capacity in the terminal. Due to this reason, it is advisable that some maintenance task be carried out at the time of the year where the least frequency of vessels exists, so that equilibrium and adequate flow are guaranteed in the container terminal. To consider this effect, in this work the parameters TTM_i,k and TDM_i,k (time duration maintenance) are introduced in the Eq. (8), then we combine the Eq. (6) and (8) using junction symbol & representing the AND logic, as we show below:

$$ \begin{aligned} & LC_{i,j}^{M} = f(t\left| \theta \right.) \\ &\quad = \left\{ {\begin{array}{*{20}l} {\overline{{LC_{i,j} }} } & {\text{if}} & {t < \sum\limits_{k = 1}^{r} {TTM_{i,k} } + \sum\limits_{k = 1}^{r - 1} {TDM_{i,k} } } \\ 0 & {\text{if}} & {\sum\limits_{k = 1}^{r} {TTM_{i,k} } + \sum\limits_{k = 1}^{r - 1} {TDM_{i,k} } \le t < \sum\limits_{k = 1}^{r} {TTM_{i,k} } + \sum\limits_{k = 1}^{r} {TDM_{i,k} } } \\ \end{array} } \right. \end{aligned} $$

(8)

where i = 1, 2, …, N_GC; j = 1, 2, …, T; k = 1, 2, …, NM_i; r = 1, 2, …, NM_i; and θ = { $\overline{{LC_{i,j} }}$, TTM_i,k, TDM_i,k}. Therefore, and as a combination consequence of both process LC = LC^D & LC^M knowing that the junction symbol & is used for the AND logic. As a modelling result, the LC variable defined in Eq. (7), consider both process, degradation (D) and maintenance (M).

Usually, SAP system store cycle life information for each device defined previously by manufactures and consequently, the information is filtered and organized in the database for maintenance purposes (DT connection).

Generators, transformers, and transmission lines modelling

The generating unit, transformer or transmission line operation is continuous, eventually fails and is repairable. This stochastic behavior can be described with Markov processes. The two-state model is used to represent the generators that operate as base load (Salgado Duarte et al. 2020), knowing this result, we used two-state Markov model and the simulation function given by (Salgado Duarte et al. 2020) for transformer and transmission line operation models case. However, in the power systems there are peak load units, intermittent operating units, or emergency units. This difference is because they are turned on and off frequently and/or it is important consider the probability of starting operation. This behavior is necessary to consider it in the generating unit model. The Sub-Committee on Application of Probabilistic Methods of the IEEE (Calsetta et al. 1972) proposed a four-state model for these generating units. In this paper, the procedure and the final simulation model used for intermittent operating units or emergency units is taken from the paper (Salgado Duarte et al. 2020). The generator, transformer, and transmission lines information treatment into the database structure for maintenance purposes (DT intermediary connection), follows the same gantry cranes description above.

Gantry cranes-generator support system modelling

Once we know the LC_i for each gantry crane i, the relation between the gantry cranes system and the generator support system is defined in Fig. 7. Two electrical systems in parallel supply the electricity for the container terminal. Power system supply the electricity under normal operating conditions, but if the service is interrupted because the transformer T2 or the transmission lines LT1 and LT2 fail or have programed maintenance, two emergency generators supply the electricity of the container terminal. We assume that the power system does not fail, this means that it is always available and the service is only interrupted by transformer T2 or the transmission lines LT1 and LT2, and the generator G1 and G2 are equal and they can, together, supply all the container terminal demand, therefore the possible states of the gantry cranes-generator support system are determinate for the following rules, where the junction symbol & is used for the AND logic while the symbol || is used for the OR logic:

Fig. 7

Gantry cranes/generator support system diagram

• State 1: Generator support system are totally unavailable; the operational capacity of the container terminal is 0% because fail or have programed maintenance, therefore (T1 & T2) || (T1 & TL1 & TL2) || (T2 & G1 & G2) || (G1 & G2 & TL1 & TL2) || (T1 & T2 & G1 & G2 & TL1 & TL2).

• State 2: Generator support system are available 50%; the operational capacity of the container terminal is 50% because fail or have programed maintenance, therefore (T2 & G1) || (T2 & G2) || (TL1 & TL2 & G1) || (TL1 & TL2 & G2).

• State 3: Generator support system are fully available; the operational capacity of the container terminal is 100%.

As a conclusion, during the simulation process, we simulate independently the LC_i for each gantry crane i, as well as all the components in the generator support system, then using the Eq. (7), we estimate the system stochastic capacity (terminal capacity), then by construction, the block diagram in Fig. 7 join the generator support system in the terminal capacity estimation.

The defined states for the generator support system allow analyzed the influence of the electricity supply reliability system in the container terminal and establishes the connection to coordinate the maintenance process of all these components at the same time using the models defined before.

Risk indicator modelling

The risk function denoted as R can be generated with the sum of X + Y random, independent, and non-negative variables. By definition, the product of R(s) = P(s)Q(s) is defined with the generating function $P\left( s \right) = \sum\nolimits_{j = 0}^{\infty } {p_{j} s^{j} }$ of X and the generating function $Q\left( s \right) = \sum\nolimits_{j = 0}^{\infty } {q_{j} s^{j} }$ of Y. Consequently, the generating function of R(s) is generally defined by the convolution formula (14):

$$ r_{k} = \sum\limits_{j = 1}^{k} {p_{j} q_{k - j} } $$

(9)

where p_j and q_j are the generated sequence from P(s) and Q(s) respectively.

In this investigation, X defined in Eq. (10) is the gantry cranes capacity distribution function of the container terminal affected for the generator support system influence, and Y defined in Eq. (11) is the containers number distribution function of the arriving vessels:

$$ X = TC $$

(10)

$$ Y = VC $$

(11)

The risk function is denoted in this investigation as R. This function is the convolution product of Eqs. (10) and (11), defined in Eq. (12):

$$ R = \left\{ {\begin{array}{*{20}l} {\frac{{\sum\nolimits_{t = 1}^{T} {Y_{t} - X_{t} } }}{{\sum\nolimits_{t = 1}^{T} {Y_{t} } }} \times 100} & {\text{if}} & {X_{t} < Y_{t} } \\ 0 & {\text{if}} & {X_{t} \ge Y_{t} } \\ \end{array} } \right. $$

(12)

The expected value of the risk function E[R] is defined in this paper as Gantry Cranes Inefficiency (GCI), when t = 1, 2, …, T considers the 8 760 h of the year (simulation window). In this work, to estimate E[R] the Monte Carlo simulation method is used. The convergence process is fluctuating in this method. However, the error level decreases when the number of samples increases, according to the law of large numbers. In this method it is not practical to run a simulation with many samples, because more calculation time is required. Therefore, it is necessary to balance the required precision and the calculation time with a stop criterion. This criterion guarantees that the simulation continues, until the risk indicator has the precision specified for the simulation. The parameter used as stopping criterion in the method is the coefficient of variation β defined in (Salgado Duarte et al. 2020).

Summarizing, the GCI indicator estimated by Monte Carlo simulation is defined below:

$$ \begin{aligned} {\text{GCI}} & = E[R] \pm \beta E[R] = E[R] \pm \frac{\sigma [R]}{E[R] \cdot \sqrt N }E[R] \\ & = E[R] \pm \frac{\sigma [R]}{\sqrt N } \end{aligned} $$

(13)

GCI indicator assess the system risk level given a maintenance scheduling proposed (DT outputs). Based on the result, the container terminal entity manager can take some decision in time (DT services), preventing container terminal capacity loss related with inefficient planification.

IMDMM: Mathematical parameterizing and calibrating

As a starting point for parametrization and calibration, it should be emphasized that the model is fully implemented in MATLAB. (2019). version 9.7.0.1,319,299 (R2019b). Natick, Massachusetts: The MathWorks Inc., and can be run on a personal computer. The implementation was tested in a system Core i5/CPU 1.6 GHz/RAM 8 GB and all results presented below were performed in the same system.

Conceptually, all the IMDMM parameters information are into the database structure for maintenance purposes (DT connection). To analyse the DT virtual counterpart applicability (optimization model), a potential scenario parameterization is described below.

The container terminal parameterization data used into the paper is taken from the work (Azimi and Ghanbari 2011). The container terminal analyzed is a system with 10 gantry cranes and each one has a number-of-movements variable following a Normal distribution functions with the average of 21 move/hour and the standard deviation of 5.56. The system analyzed assumes 10 gantry cranes available with the TTR for each gantry cranes following a normal distribution with the average of 23.6520 h and the standard deviation of 12.2640 (Soszynska 2012), the related index of TTF for all gantry cranes following Weibull distribution with scale parameter of 1219 h and shape parameter of 0.84 (Azimi and Ghanbari 2011). The parameterization proposed for the generators, power transformer and transmission lines are shown in the Table 3 below.

Table 3 Components data

In the paper (Azimi and Ghanbari 2011), the period between the arrival of two consecutive vessels follows an Exponential distribution with the average of 9.41 h, but we change the approach variating the average value depending on container terminal operation year seasons. In order to evidence the proposed approach consistency, we describe visually all the modelling steps, and, in the meantime, we analyze the results as a calibration. Figure 8 shows the simulated containers demand in a year according to the vessel’s frequency in the Table 4 that arrive at the container terminal and show how many containers we must move each hour with the gantry cranes system. It is visible in Fig. 8, year moments with arriving high vessel’s frequency.

Fig. 8

Containers demand

Table 4 Frequency according to the time of year

Figure 9 shows the simulated terminal stochastic capacity with 10 gantry cranes according to the vessel’s frequency, that is, the capacity conjunction of all the gantry cranes per hour, giving us the container terminal capacity must move per hour, according to the gantry cranes availability affected by the generator support system. Figure 9 is consistent because when the vessel frequency is higher, the less time the vessel stays in the port, as a result the simulated operating terminal capacity has higher frequency (gantry cranes with higher operating frequency) and less holistic capacity available, therefore, more time the vessel in the port, more time available to move the load, then more capacity available.

Fig. 9

Container terminal stochastic capacity

Figure 10 shows the convolution between the simulated containers demand in a year according to the arriving vessels frequency and the simulated stochastic capacity of the container terminal with 10 gantry cranes. In addition, Fig. 11 shows the relation between X, Y and R in the first 200 h of the year for the study case discussed. The two figures above are general views of the simulation approach because measure the capacity loss in the container terminal when the vessel frequency overload the movement availability of the gantry cranes.

Fig. 10

Convolution curve for a simulated year

Fig. 11

Convolution curve of 200 simulated hours

One of the paper objectives is to highlight the influence of maintenance on the system. To analyze how the maintenance process affects the system, a scenario is proposed where, for each gantry crane, generator support system, power transformer and transmission lines, two weeks of continuous preventive maintenance are scheduled for the year, and it is assumed that the components will not be available for the terminal during the two weeks of maintenance.

Given the above parametrization of the model, is possible visualize the relation between all the components in the system and merge all the steps described in the “IMDMM: Mathematical parameterizing and calibrating” section with the parameters values in the current section, therefore, the maintenance scheduling behind of the Integrated Maintenance Decision Making Model is ready to be solve in the next section.

Before we jump to the next section, Fig. 12 shows how the first condition in the flow diagram is guaranteed given a maintenance scheduling of all the components in the system. In this case, 62 iterations are necessary to obtain the desired criterial. In addition, the minimum value that the optimization algorithm can reach, given the described parametrization is GCI = 10.54%. The target value is known because the best possible scenario is no maintenance at all.

Fig. 12

Expected and simulated value behavior

In the case of the simulation parameters, the simulation window is one year (8760 h) and we assume robust expected value estimation (GCI indicator) by Monte Carlo simulation when ε = 0.01 (Salgado Duarte et al. 2020).

In order to estimate the expected value of the risk function, we must know the system distribution function, but in this case the construction is only possible by bootstrapping the individual’s distribution functions by Monte Carlo simulations (gantry cranes individual’s distribution) and then based on the rules described before, the system distribution function is created (see Fig. 7), this technique is known as Parametric Bootstrap. When the complexity and dimension of the system is large, simulation time is an issue. Equation (14) declare the criterial for robust expected value estimation (optimization target) and the clear dependency on number of simulations and expected value standard deviation. Therefore, the number of simulations needed to guarantee a robust estimation depends on the scenario evaluated, but for all the cases the same rule is considered (ε = 0.01), as Fig. 12 shows.

IMDMM: Mathematical solving

As we describe above, the model has a stochastic non-linear objective function with bounded constraints. In order to solve this specific problem, PSO algorithmic was used because GCI risk indicator (objective function value given the maintenance scheduling (DT services), is the results of a convolution by Monte Carlo simulation, therefore we do not know the objective function derivate and Newton’s, Lagrange, quasi-Newton or Sequential Quadratic Programming traditional methods cannot be used.

PSO is fully applicable to this problem and a brief description below clarify the reason. The next subsections we describe the PSO algorithmic implemented and the solution proposed for the previous parametrization.

Particle swarm optimization (PSO) algorithm

The PSO algorithm used in this investigation is based on the algorithm described in (Kennedy and Eberhart 1995), using modifications suggested in (Mezura-Montes and Coello Coello 2011) and in (Pedersen 2010). In addition, the algorithm was successfully tested with the Rosenbrock’s function or banana function.

The particle swarm algorithm begins by creating the initial particles and assigning them initial velocities. It evaluates the objective function at each particle location and determines the best (lowest) function value and the best location. It chooses new velocities, based on the current velocity, the particles’ individual best locations, and the best locations of their neighbors. It then iteratively updates the particle locations (the new location is the old one plus the velocity, modified to keep particles within bounds), velocities, and neighbors. Iterations proceed until the algorithm reaches a stopping criterion.

By default, the algorithm creates particles at random uniformly within bounds (model constrains, restriction intervals for TTM_i,1). If there is an unbounded component, the algorithm creates particles with a random uniform distribution from –1000 to 1000. If you have only one bound, the algorithm shifts the creation to have the bound as an endpoint, and a creation interval 2000 wide. Particle j has position x(j), which is a row vector with n_vars elements (number of TTM_i,1 to be coordinated). Control the span of the initial swarm using the TTM ⁰_i,1 (initial swarm span).

Similarly, the algorithm creates initial particle velocities v at random uniformly within the range [-r, r], where r is the vector of initial ranges. The range of component i is min(upper-bound(i)—lower-bound(i), initial-swarm-span(i)), in our case, min((Simulation Window − (TTM_i,1 + TTM_i,k≠1 + TDM_i,k)) − 0, TTM ⁰_i,1 ).

The algorithm implemented evaluates the objective function at all particles. It records the current position p(j) of each particle j. In subsequent iterations, p(j) will be the location of the best objective function that particle j has found. And b is the best overall particles: b = min (objective-function(p(j))) and d is the location such that b = objective-function(d), in our case, the objective-function E[R] = f (TTM ^j_i,1 ).

The algorithm initializes the neighborhood-size N to N_min = max (round (swarm-size × min-neighbors-fraction)). The inertia W = max(inertia-range), or if inertia-range is negative, it sets W = min(inertia-range) and the stall counter c = 0. For convenience of notation, set the variable y_local = self-adjustment-weight, and y_global = social-adjustment-weight, where self-adjustment-weight and social-adjustment-weight are options.

The algorithm implemented updates the swarm as follows by iteration steps. For particle j, which is at position x(j):

• 1. Choose a random subset S of N particles other than j.

• 2. Find f_best(S), the best objective function among the neighbors, and g(S), the position of the neighbor with the best objective function.

• 3. For u₁ and u₂ uniformly (0,1) distributed random vectors of length n_vars, update the velocity v = W·v + y_local·u₁·(p-x) + y_global·u₂·(g-x). This update uses a weighted sum of:

  • a. The previous velocity v

  • b. The difference between the current position and the best position the particle has seen p-x

  • c. The difference between the current position and the best position in the current neighborhood g-x

• 4. Update the position x = x + v.

• 5. Enforce the bounds. If any component of x is outside a bound, set it equal to that bound. For those components that were just set to a bound, if the velocity v of that component points outside the bound, set that velocity component to zero.

• 6. Evaluate the objective function f = function(x).

• 7. If f < function(p), then set p = x. This step ensures p has the best position the particle has seen.

• 8. The next steps of the algorithm apply to parameters of the entire swarm, not the individual particles. Consider the smallest f = min(f(j)) among the particles j in the swarm. If f < b, then set b = f and d = x. This step ensures b has the best objective function in the swarm, and d has the best location.

• 9. If, in the previous step, the best function value was lowered, then set flag = true. Otherwise, flag = false. The value of flag is used in the next step.

• 10. Update the neighborhood.

  If flag = true:

  • a. Set c = max (0, c−1).

  • b. Set N to N_min.

  • c. If c < 2, then set W = 2·W.

  • d. If c > 5, then set W = W/2.

  • e. Ensure that W is in the bounds of the inertia-range.

    If flag = false:

    • a. Set c = c + 1.

    • b. Set N = min (N + N_min, swarm-size).

The algorithm iterates until it reaches a stopping criterion, in this case, when the relative change in the best objective function value g over the last max-stall-iterations is less than 0.001 (Function Tolerance described in the diagram flow). The algorithm settings used in the solution is described in the following Table 5.

Table 5 Optimization algorithm properties

Solution proposed

Once the optimization algorithm used on the solution and the full parameterization are described, the results of the proposed optimization model are shown in Fig. 13 and Table 6.

Fig. 13

Convergence process of the optimization algorithm

Table 6 Solution for the scenario

Figure 13 is the convergence process of the optimization model and shows how the GCI indicator varies when gantry cranes, generator support system, power transformers and transmission lines maintenance scheduling changes during the year. The value of the objective function achieved (GCI = 11.48%) shows the proposed approach robustness because the result is close to GCI = 10.54% (no maintenance at all).

Proper planned maintenance scheduling process improve the operational efficiency in the container terminal, and the GCI indicator measure the impact. Table 6 throws the expected results, the support system component scheduling maintenance must be carried out during the periods of the year where the vessels have the lowest frequency in the container terminal, therefore, the GCI risk indicator is valuable for the container terminal managers (DT services) because they can decide according to the risk level of the decision-making process, what would be the best moment in the year to carry out the maintenance process in the system (DT counterpart).

Conclusions

The paper describes the proposed IMDMM concept (DT counterpart) applicability based on risk assessment, and how the model can efficiently find the maintenance scheduling for all the devices under operation considered in practice. The above results confirm, in the presented scenario, the relation between the estimation of GCI indicator (DT services) used into a container terminal by the entity manager and the vessels arriving frequency behavior to the harbor (free variable in the DT physical object). Experimental results show that presented model can help to organize the maintenance scheduling strategy in the container terminal. The paper solves the assessing risks problem of transportation process in the container terminals through the simulation-based approach (MCMC) which considers the relationship between random factors (historical degradation data fitted by machine learning framework) during the container reloading and maintenance scheduling process efficient used in the container terminal capacity estimation. The presented DT model opens the way to extensive simulations under various scenarios and conditions, with the possibility to be updated in real-time, to detect anomalies and to conduct accurate diagnostics and prognostics of cranes into selected scenarios.