Fractional Calculus to Analyze Efficiency Behavior in a Balancing Loop in a System Dynamics Environment

Abstract

This research project focuses on developing a mathematical model that allows us to understand the behavior of the balancing loops in system dynamics in greater detail and precision. Currently, simulations give us an understanding of the behavior of these loops, but under the premise of an ideal scenario. In practice, however, accurate models often operate with varying efficiencies due to various irregularities and particularities. This discrepancy is the primary motivation behind our research proposal, which seeks to provide a more realistic understanding of the behavior of the loops, including their different levels of efficiency. To achieve this goal, we propose the introduction of fractional calculus in system dynamics models, focusing specifically on the balancing loops. This innovative approach offers a new perspective on the state of the art, offering new possibilities for understanding and optimizing complex systems.

1. Introduction

System dynamics (SD), as a modeling and simulation discipline, has emerged as a fundamental tool in solving complex and dynamic problems across various fields from business management to environmental science. These models, grounded in systems theory and computer simulation, empower researchers and professionals to comprehend and analyze the behavior of dynamic systems over time [1]. The primary objective of this technique is to capture causal relationships and mechanisms underlying observed phenomena, providing an invaluable perspective for decision-making and addressing complex challenges in a world characterized by interconnection and interdependence [2].

System dynamics (SD) emerges as a rigorous and systematic methodology designed for the modeling and analysis of complex systems, enabling a profound understanding of how interactions among their components result in dynamic behavior. Through simulation models, this discipline facilitates the identification of long-term patterns and trends and assesses how various actions and decisions can impact systemic behavior [3]. In this context, sensitivity analysis emerges as a highly beneficial tool employed across various disciplines from economics to engineering and science. Its purpose is to assess how changes in the values of model variables influence the obtained results, providing valuable information crucial for making informed decisions [4]. Sensitivity analysis makes it possible to identify the relatively small set of parameters whose values significantly alter the model’s behavior or responses to different policies. In this way, parameters worth calculating more accurately are discovered [5].

Numerous highly relevant applications have been developed in various fields of knowledge, such as health, the environment, education, and agriculture, that leverage the sensitivity analysis technique in SD models. These applications have yielded positive results in understanding behaviors and obtaining models that meet the desired expectations [6,7,8,9,10,11,12,13].

Undoubtedly, SD models pose significant challenges. Selecting the appropriate parameters for a sensitivity analysis is intricate and often relies on subjective judgments. This situation underscores the urgent need to develop a more analytical and objective approach to identify the parameters that impact system behavior with greater clarity and simplicity. There is a significant opportunity to enhance existing mathematical models in SD. Understanding models and, therefore, reality from a feedback structure-based perspective is of utmost importance for the SD method. For years, the field has been developing tools and mathematical methods to perform automated and objective loop dominance analysis cf. [14,15,16,17,18,19,20].

The research conducted by [15] focuses on the “path participation” technique in system dynamics, systematically analyzing component interactions and assessing feedback loops’ importance. On the other hand, the paper presented by [16] explores metrics and analysis of the elasticity of eigenvalues in SD models, highlighting how different approaches converge inconsistent results. In addition, the study carried out by [18] extends the scope of methods to identify structural dominance, applying these techniques in more complex and stochastic models and providing coherent information on the system’s behavior. The three investments focus on identifying the most dominant or contributing loop in SD models. However, a common weakness is the need for a definite global mathematical model of the system, which limits analysis from a different perspective. This limitation highlights the need to approach the system from broader and more detailed approaches, as proposed in this research.

The research conducted by [17] focuses on analyzing the internal dynamics of state variables in dynamic systems, highlighting the importance of expressing time trajectories through linear combinations of eigenvectors. This approach offers a valuable perspective on the dynamics of complex systems by reducing the model to a time-invariant state space. However, there are still opportunities to address variable efficiencies in precise models, which motivates our research proposal. On the other hand, the article [19] addresses the relevance of the “Loops that Matter” (LTM) approach in understanding behaviors in SD models. Although LTM is widely applicable and easy to use, it has limitations in its original formulation, particularly regarding the impact of a flow in an accumulator on net flow. Although the article proposes reformulation to address this limitation, an integral mathematical model of the system is still needed. Our proposal complements this work by focusing on the detailed identification of balance loops, addressing the need for more precise and complete dynamic system models.

The authors of [14,20] propose innovative methodologies for transforming system dynamics simulation models into mathematical models expressed as system transfer functions. While [14] uses the Laplace transform to represent the system in the frequency domain and thus expands the understanding of models by introducing differential equations, Ref. [20] employs state space control and representation engineering to obtain a linear and time-invariant mathematical model. Both contributions strengthen simulation models in dynamic systems by offering analytical alternatives through modern control engineering, enabling more informed and appropriate decision-making. However, there is a complementary area to strengthen SD models where rolling loops are present, which is the focus of this research project that seeks to understand the behavior of equilibrium loops in greater detail and precision by introducing fractional calculus.

Current research focuses on balancing loops, also known as negative loops. These components are fundamental in the search for balance within the system, as they are designed to adjust the system’s state toward a desired objective or state. Its primary function is to counteract any disturbance that may divert the system from that objective, thus ensuring stability and optimal operation.

Negative feedback loops operate by detecting deviations from the desired state and implementing corrective actions to restore balance. This self-regulation capability allows systems to maintain stable and predictable behavior over time. All negative feedback loops share a standard structure, as illustrated in Figure 1, which facilitates their understanding and analysis in various contexts and disciplines cf. [21].

Figure 2 illustrates the typical behavior of a rolling model, where a production target for parts in the textile industry is established. The curve represents the system’s monitoring of this target, reflecting the ideal and desired behavior.

The main objective of this research project is to develop a mathematical model that allows a more detailed and precise understanding of the behavior of equilibrium loops in SD. These loops are vital to maintain stability and meet specific objectives within a simulation model.

The behavior of these balancing loops is usually based on ideal scenarios that do not fully reflect the complexity of reality. In practice, precise models can operate with varying efficiencies due to various irregularities and particularities. This discrepancy is the main engine of our research proposal, which seeks to offer a more realistic understanding of the behavior of the loops, considering their different levels of efficiency.

To achieve this objective, we propose the introduction of fractional calculation in SD models, with a specific focus on equilibrium loops. This innovative approach opens up new perspectives on the state of the art by offering a more accurate way to model and understand complex systems. By including fractional calculation, we can better capture the nonlinear nature and temporal dynamics of equilibrium loops, allowing us to effectively address efficiency variations and other irregularities in practice.

In this research, we focus on two departments of the textile process, which we represent as two balance loops. To obtain accurate and applicable data, we collaborate with a company in the southern state of Guanajuato. By analyzing these loops, we evaluate the efficiencies of the Overall Equipment Efficiency (OEE) model, providing valuable information for business managers’ decision-making. The final objective is to contribute to developing methodologies and tools that allow more effective and efficient management of complex systems in the industrial field and, in general, to the dynamics of systems.

Fractional calculus (FC) is a branch of mathematics based on the concept of derivatives and integrals of non-integer or fractional orders. This mathematical discipline has found a wide range of applications in various fields of knowledge from physics and biology to engineering and economics. One of the most prominent features of fractional calculus is its ability to model and understand complex phenomena that exhibit long-term dependencies and nonlinear behaviors.

In the realm of modeling physical systems, fractional calculus has been used to describe phenomena such as anomalous diffusion in porous media, particle dynamics in fractal media, and wave propagation in dispersive media. For example, in a study conducted by Podlubny [22], fractional calculus was applied to model contaminant diffusion in aquifers, allowing for a more accurate capture of contaminant propagation in complex hydrogeological systems.

Fractional calculus has proven helpful in understanding biological phenomena such as population dynamics, substance diffusion in biological tissues, and neuronal activity modeling. In a landmark paper by Mainardi [23], the applications of fractional calculus in modeling population dynamics were explored, highlighting its ability to capture the complexity of evolving biological systems.

In signal processing, fractional calculus has been used to model and analyze non-stationary signals and stochastic processes. For instance, in a study by Magin [24], fractional calculus was applied to develop nonlinear time series models, leading to a better understanding of temporal variability in experimental data.

Another innovative approach in fractional calculus is used to improve the accuracy of wind energy prediction [25], which is crucial for optimizing the operation of the electrical grid and harnessing renewable resources. Researchers propose a short-term memory model to forecast missing wind speed and direction data, and a fractional-order neural network (FONN) with a fractional activation function to enhance the prediction of generated wind energy. Through two case studies, the predictive effectiveness of the FONN model in wind power prediction is demonstrated. This hybrid approach improves forecast accuracy and addresses data gaps evaluated through mean errors and R² values.

Another innovation involves presenting a new approach to numerically solving systems of fractional integro-differential equations [26]. Vieta–Fibonacci polynomials are used as essential functions, and the projection method for the fractional Caputo order is derived for the first time. This innovative approach involves an efficient transformation that reduces the problem to a system of two independent equations whose solution approximates the original problem. The efficiency and accuracy of the proposed method are validated, demonstrating the existence of a solution to the approximate problem and performing an error analysis. Numerical tests support theoretical interpretations.

Currently, FC is applied in a wide range of knowledge areas, including the analysis of viscoelastic models [27,28], fractional models of anomalous diffusion cf. [29], analysis of fractional-order electrical circuits as in [30], and medical image processing [31,32]. Additionally, the study of the role of vaccines in controlling the spread of COVID-19 involves fractional models of COVID-19 behavior when vaccines are applied [33]. These are just a few examples of applications.

Introducing fractional calculus into SD models represents a significant advancement with essential implications across various fields, including industrial engineering. This innovation enables a more precise description of complex systems, addressing nonlinear phenomena and long-term dependence. It results in a deeper understanding of constantly changing systems and more effective solutions in various disciplines.

In this study, we focused on using fractional calculus to describe the behavior of equilibrium loops in the presence of irregularities and particularities that hinder ideal behavior. However, it is essential to note that fractional models present an inherent limitation: they do not provide a precise definition of differential or geometry, which raises questions about their effectiveness. Several proposals have been made in recent years that attempt to give a physical and geometric interpretation to derivatives and integrals of non-integer order. Although there is no unified criterion, as in the case of integer operators, some phenomena, such as viscoelasticity, could not be explained in terms of whole operators. Moreover, violating the Leibniz rule allows fractional operators to be used in practical applications. In our work, fractional models seek to interpret the model’s efficiency (balancing loop) where the ideal case would correspond to an order of derivative equal to 1. Although there is no unified criterion regarding fractional calculus’s physical and geometric interpretation, some proposals have been made, such as [34,35,36]. On the other hand, Tarasov [37] proposes that one of the rules that fractional operators must meet is violating the Leibniz rule; this allows fractional calculation to be used in various physical applications.

Significant studies have focused on studying systems of integral-differential fractional equations in fractional calculus. These investigations stand out for their resolution through innovative numerical methods, simplifying the search for solutions [38,39]. However, the main strength of our research project lies in the fact that we offer an analytical solution to this type of equation, which represents a remarkable advance over previous approaches mainly focused on numerical strategies.

It is a versatile tool with multidisciplinary applications ranging from modeling physical and biological systems to process optimization. It has already been successfully applied in areas such as viscoelastic model analysis, anomalous diffusion, electrical circuits, and the study of COVID-19 propagation when vaccines are applied. This approach not only enriches the field of dynamic systems but also significantly impacts various areas of science and technology.

2. Materials and Methods

Fractional calculus (FC) is a discipline that studies fractional order integration and differentiation operators, meaning operators of non-integer order. Its origin dates back to the beginnings of conventional calculus when the first notations for derivatives were proposed [40]. Over the years, different definitions of fractional derivative have emerged, some based on the generalization of the limit definition of differentiation, such as the definition by Lacroix, the fractional derivative [41], and in recent years, the conformable derivative [42]. Others are based on the generalization of multiple integrations and the properties of semigroup, namely, the Riemann–Liouville fractional integral (IRL) and derivatives: Riemann–Liouville (DRL), Caputo derivative [43], Caputo–Fabrizio derivative (DCF) proposed in 2015 [44], and the Atangana–Baleanu derivative (AB) from 2016 [45].

2.1. Gamma Function

The factorial of an integer number $n$ is defined as the product of all consecutive integers from 1 to $n$ and is expressed:

$$n!=1 · 2 · 3···\left(n-1\right)· n,$$

(1)

where zero factorial is defined as 0! = 1. For an extension for values other than integers, use the Gamma function defined by the integral [46]:

$$\Gamma (z)=\underset{0}{\overset{\infty }{\int }}{t}^{z-1}{e}^{-t}dt$$

(2)

for any $z$∈ $C$, as to $\mathfrak{R}\left(z\right)>0$, when integrating by parts, you have an attractive property called recurrence property:

$$\Gamma \left(z+1\right)=z\Gamma \left(z\right) ,$$

(3)

When applied to integer real values, we have the relation:

$$n!=\Gamma \left(n+1\right),$$

(4)

in this way, the factorial concept can be generalized for real or complex values.

2.2. Mittag-Leffler Function

Just as the exponential function plays a fundamental role in solving many ordinary differential equations, the Mittag-Leffler function (ML) [47], defined in 1903, consistently appears in the solution of fractional differential equations.

$$E_{\alpha }\left(z\right) = \sum _{k=0}^{\infty } \frac{z^k}{\Gamma \left(\alpha k + 1\right)},$$

(5)

where $\alpha$, $z$ ∈ $C$, α is an arbitrary parameter; subsequently, a generalization was proposed in 1905, and other updates in [48,49]. The function is

$$E_{\alpha ,\mathsf{\beta }}\left(z\right) = \sum _{k=0}^{\infty } \frac{z^k}{\Gamma \left(\alpha k + \mathsf{\beta }\right)},$$

(6)

where $e$ α, $\beta$∈ $C$ are also arbitrary parameters, and the expression is called a two-parameter ML function.

$${E^{\mathsf{\gamma }}}_{\alpha ,\mathsf{\beta },}\left(z\right) = \sum _{k=0}^{\infty } \frac{\left(\mathsf{\gamma }\right)\mathrm{k}}{\Gamma \left(\alpha k + \mathsf{\beta }\right)} \frac{z^k}{k!},$$

(7)

where again, α, β, and γ C are arbitrary parameters, and the expression is called the Mittag-Leffler function of three parameters. When γ = 1, the equation becomes Equation (6). In addition, it can be observed that Equation (5) is a special case when β = 1 in the Equation (5). Therefore, the exponential function is considered a special case of ML when α = 1, that is:

$${E^1}_{\mathrm{1,1}}\left(z\right) = \sum _{k=0}^{\infty } \frac{z^k}{\Gamma \left(k + 1\right)}=\sum _{k=0}^{\infty } \frac{z^k}{k!}=e^z,$$

(8)

The aforementioned is why the ML function and its generalizations are called generalized Mittag-Leffler functions. Research into the properties of the ML function has continued [50,51]. Let us consider the Laplace transforms of the following functions: Equations (9) and (10) in terms of the Mittag-Leffler function [52]:

$$\mathcal{L}\left\{t^{\alpha }\right\}=\frac{\mathsf{\Gamma }\left[\alpha +1\right]}{s^{\alpha +1}},$$

(9)

$$\mathcal{L}\left\{\left[1-E_{\alpha }\left(\omega t^{\alpha }\right)\right]\right\}=\frac{\omega }{s\left(s^{\alpha }+\omega \right)},$$

(10)

where $\omega \in \mathbb{R}$ y $\alpha \in \mathbb{C}$, with $\mathfrak{R}\left(\alpha \right)>0$.

2.3. Caputo Derivative

The fractional Caputo derivative was proposed in 1971 [43] to avoid the fractional order initial conditions of the Riemann–Liouville derivative. It is defined as follows: see Equation (11).

Definition 1.

Let  $f\left(t\right):t\to t$ be a continuous function  $(a,b)$

$${}_0{}^C{D}_t^{\alpha }f\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(n-\alpha \right)}{\int }_a^t\frac{f^{\left(n\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha }}d\tau ,$$

(11)

where $a\in \mathbb{R},n-1<\alpha <n$, y $f\in H^n\left(a,b\right), b>a$. And his Laplace transform is given by:

$$\mathcal{L}\left\{{}_0{}^C{D}_t^{\alpha }f\left(t\right)\right\}=s^{\alpha }F\left(s\right)-\sum _{k=0}^{n-1}{s}^{\alpha -k-1}{f}^{\left(k\right)}\left(0\right),$$

(12)

where $\mathcal{L}\left\{f\left(t\right)\right\}=F\left(s\right)$, $f^{\left(m\right)}\left(0\right)={\left.\frac{d^{m}f\left(t\right)}{d{t}^m}\right|}_{t=0}$ and $m\in [0,n]$. for n = 1 the Laplace transform is given by:

$$\mathcal{L}\left\{{}_0{}^C{D}_t^{\alpha }f\left(t\right)\right\}=s^{\alpha }F\left(s\right)-s^{\alpha -1}f\left(0\right).$$

(13)

3. Methodology

In this investigation, the process begins with creating the system dynamics model, which is based on identifying and analyzing equilibrium loops and sets a specific goal to achieve. Then, we formulate the system of differential equations that govern the system’s behavior. Ordinary differential equations initially solve this system to obtain an initial solution. The OEE method is proposed, and the efficiencies of each department, which relate to alphas, are obtained. Later, fractional calculus is implemented using the Caputo method to solve fractional differential equations, allowing a more accurate description of the system’s behavior. Simulations are performed using these mathematical models for both the Tissue and Threading departments, providing a thorough evaluation of the system’s performance in different operating scenarios and conditions. This comprehensive methodological approach ensures a thorough and detailed understanding of the dynamics of the industrial processes studied, as shown in Figure 3.

3.1. System Dynamics Model

The Forrester diagram, as shown in Figure 4, represents the dynamics of the textile production process through an SD model. This model shows the two fundamental departments of the textile process: Weaving and Basting. It comprises two state variables, x and y, representing these departments, respectively. These state variables are the system’s central axis, reflecting each department’s cumulative production. The diagram also incorporates six auxiliary variables and two flow variables, forming the system, interactions, and feedback mechanisms.

The auxiliary variables, Xd1 and Xd2, are particularly prominent as they mean the system objectives or the target number of parts (for example, parts of a specific sweater model) that each department seeks to produce. These goals are not static; they guide operational dynamics, shaping workflow and production efforts within each department. The diagram also introduces adders 1 and 2, which calculate the discrepancy between the actual production and the targets set for each department, highlighting areas where adjustments may be needed to align with production targets.

The constant k, assigned to each department, represents the ideal rate at which the desired production level can be reached, expressed in units of (1/min). These constants are essential to calibrate the system to optimal performance, ensuring that production targets are met efficiently.

Table 1. Variables of the Forrester diagram.

State Variable	Notation of the Variable	Xd (Goals)	Adder	Differential Equation
Weaving	$X$	$Xd1=72$	$Xd-x$	$Ed1=(Xd1-x)k1$
Basting	$Y$	$Xd2=72$	$Xd2-y$	Ed$2=(\mathrm{Xd}2 -$ y) k2 − Ed1

describes these variables’ specific values and roles in the Forrester diagram, providing a clear framework for understanding the dynamics at play.

3.2. System of Differential Equations

It is important to note that a system dynamics model has been developed to represent the production process in two departments. From this model, the two fundamental differential equations that describe the behavior of these departments, presented as Equations (14) and (15), have been identified and formulated. These equations form the analysis’s backbone, providing a detailed understanding of production dynamics in these specific contexts.

The proposed system, with initial conditions $x\left(0\right)=0, y\left(0\right)=0$,

$$\begin{array}{c} \frac{dx}{dt}=k_1\left(x_{d1}-x\right), \\ \end{array}$$

(14)

$$\frac{dy}{dt}=k_2\left(x_{d2}-y\right)-k_1\left(x_{d1}-x\right).$$

(15)

3.3. Solution Ordinary Differential Equation

Equations (16) and (17) are obtained by solving the whole system of differential equations.

$$x=x_{d1}\left(1-e^{-k_{1}t}\right),$$

(16)

$$y=x_{d2\left(1-e^{-k_{2}t}\right)}+\frac{x_{d1}{k}_1}{k_2-k_1}\left(e^{-k_{2}t}-e^{-k_{1}t}\right)$$

(17)

To obtain the values of $k$, the solution of each of the state variables must be equal; it is equal to the operating time of each department, where they are 480 and 100 min, respectively; the value of k is sought, where 98% of the goal is obtained, obtaining a value of k1 = 0.009541 for the Weaving department and k2 = 0.12 for the Basting department.

3.4. OEE Methodology

The methodological approach of this study focuses on the Weaving and Basting departments, which are fundamental segments of the textile production chain. Working with a company located in the southern part of Guanajuato, real-world data were collected to assess these departments’ efficiency accurately. The evaluation leverages the Overall Equipment Effectiveness (OEE) methodology, a comprehensive measure that evaluates a department’s efficiency based on equipment availability, performance speed, and product quality. For this study, product quality is specifically related to the production of a particular model of sweater, underlining the practical application and relevance of this research for the operational challenges of the textile industry.

In this research project, fractional differential equations and their alpha value, which will be addressed, are intrinsically linked to the efficiency derived from the OEE methodology. Therefore, obtaining the absolute efficiency values of the company located south of Guanajuato is essential. These efficiency values are directly related to the fractional order of the derivative, underlining the importance of their accuracy and validity for proper modeling and understanding of the dynamics of the textile production process.

3.4.1. Availability Factor

Availability is calculated by dividing the actual operating time by the planned production time. This factor, expressed as a percentage, provides a key measure of operational efficiency by considering the time the equipment operates according to the established plan. The availability of the Weaving department for the week of October 25 to 30 can be observed in

Table 2. Availability of the Weaving department.

Day	25	26	27	28	29	30
Total time (min)	690	690	690	690	690	630
Breaks (min)	120	120	120	120	120	120
Maintenance stoppage (min)	10	10	10	10	10	45
Shutdowns of machine records (min)	20.43	29.21	36.23	20.78	32.45	30.27
(A) Planned time available (min)	560	560	560	560	560	465
(B) Productive time (min)	539.5	530.7	523.7	539.2	527.5	434.73
(B/A) Availability (%)	96.35	94.78	93.53	96.29	94.21	93.49

.

Table 3. Availability of the Basting department.

Day	25	26	27	28	29	30
Total Time (min)	100	100	100	100	100	90
Breaks (min)	20	20	20	20	20	20
(A) Planned time available (min)	80	80	80	80	80	70
(B) Productive time (min)	79.26	79.35	76.55	78.41	79.95	69.73
(B/A) Availability (%)	99.10	99.19	95.70	98.03	99.95	99.63

shows the availability of the Basting department for the week of 25–30 October.

3.4.2. Quality Factor

Quality in Overall Equipment Effectiveness (OEE) is measured by considering the units that meet quality standards against the total units produced. The quality factor of the Weaving department for the week of 25–30 October can be observed in

Table 4. Quality of the Weaving department.

Day	25	26	27	28	29	30
(A) Actual production- total parts	80	79	72	72	69	55
Contaminated canvases	0	1	0	2	0	0
Overlay	5	6	5	10	7	10
Non-functional canvases	1	3	3	5	3	7
(B) Good pieces	74	69	64	55	59	38
(B/A) Quality (%)	92.50	87.34	88.89	76.39	85.51	69.09

.

Table 5. Quality of the Basting department.

Day	25	26	27	28	29	30
(A) Actual production-total parts	75	64	72	79	88	81
(B) Good pieces	75	64	72	78	88	80
(B/A) Quality (%)	100	100	100	98.73	100	98.77

shows the quality factor of the Threading department of the week from 25–30 October.

3.4.3. Performance Factor

Efficiency is measured by the planned production factor divided by the actual production.

Table 6. Performance of the Weaving department.

Day	25	26	27	28	29	30
(A) Planned production (total pieces)	82.7	79.3	72.98	73.44	69.24	55.21
(B) Real production	80	79	72	72	69	55
(B/A) Performance (%)	96.70	99.59	98.66	98.04	99.65	99.63

shows the efficiency factor of the Basting department for the week of 25–30 October.

The performance factor of the Basting department for the week of 25–30 October can be observed in

Table 7. Performance of the Basting department.

Day	25	26	27	28	29	30
(A) Planned production (total pieces)	78.0	78.1	75.3	83.7	96.5	84.2
(B) Real production	75	64	72	79	88	81
(B/A) Performance (%)	96.10	81.93	95.54	94.38	91.16	96.20

.

The OEE factor is determined by multiplying the three key factors in the methodology.

Table 8. Calculation of the OEE methodology in the departments.

Day	25	26	27	28	29	30	OEE Weekly
OEE (Weaving) (%)	86.1	82.4	82.0	72.1	80.2	64.3	77.90
OEE (Basting) (%)	95.2	81.2	91.4	91.3	91.1	94.6	90.84

details the efficiency calculation for the two departments.

Overall Equipment Effectiveness (OEE) was calculated for both departments, yielding an alpha value of 0.7790 for the tissue process and an alpha value of 0.9084 for the Threading department.

4. Results

In this work, the results obtained through the use of the derivative of Caputo are presented, accompanied by their demonstration and the corresponding graphic representations in function of time. These graphs show variations in system behavior under different levels of efficiency.

4.1. Solution of The Fractional Order Model Per Caputo

It is proposed to replace the first-order derivative operators of the system of differential equations with fractional-order operators of order α, where 0 < α < 1. When α = 1, it returns to the first-order operators and, therefore, to the original system of equations, see Equations (18) and (19).

$$\frac{d^{\alpha }x}{d{t}^{\alpha }}=k_1\left(x_{d1}-x\right),$$

(18)

$$\frac{d^{\alpha }y}{d{t}^{\alpha }}=k_2\left(x_{d2}-y\right)-k_1\left(x_{d1}-x\right).$$

(19)

When using the Caputo derivative operator in the system of fractional equations, we have:

$${}_0{}^C{D}_t^{\alpha }x=k_1\left(x_{d1}-x\right),$$

(20)

$${}_0{}^C{D}_t^{\alpha }y=k_2\left(x_{d2}-y\right)-k_1\left(x_{d1}-x\right).$$

(21)

To solve the system of fractional order differential equations, we apply the Laplace transform:

$$\left[s^{\alpha }X\left(s\right)-s^{\alpha -1}x\left(0\right)\right]=\frac{k_1{x}_{d1}}{s}-k_{1}X\left(s\right),$$

(22)

$$\left[s^{\alpha }Y\left(s\right)-s^{\alpha -1}y\left(0\right)\right]=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{s}-k_{2}Y\left(s\right)+k_{1}X\left(s\right).$$

(23)

By simplifying both equations and applying the initial conditions, we obtain:

$$\begin{gathered} \begin{array}{c}X\left(s\right)\left[s^{\alpha }+k_1\right]=\frac{k_1{x}_{d1}}{s}\\ \end{array}, \\ Y\left(s\right)\left[s^{\alpha }+k_2\right]=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{s}+k_{1}X\left(s\right). \end{gathered}$$

(24)

When resolving to $X\left(s\right)$:

$$X\left(s\right)=x_{d1}\frac{k_1}{s\left(s^{\alpha }+k_1\right)},$$

(25)

we can find the solution by $x\left(t\right)$ applying the inverse Laplace transform using the formula of Equation (10), so we have Equation (26).

$$x\left(t;\alpha \right)=x_{d1}\left[1-E_{\alpha }\left({-k}_1{t}^{\alpha }\right)\right],$$

(26)

To find the solution of $y\left(t\right)$, we can use $X\left(s\right)$ of Equation (25):

$$Y\left(s\right)\left[s^{\alpha }+k_2\right]=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{s}+x_{d1}\frac{k_1^2}{s\left(s^{\alpha }+k_1\right)},$$

(27)

Clearing for Y(s), we have:

$$Y\left(s\right)=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{s\left(s^{\alpha }+k_2\right)}+x_{d1}\frac{k_1^2}{s\left(s^{\alpha }+k_1\right)\left(s^{\alpha }+k_2\right) }.$$

(28)

To solve Y(s), we separate in two terms $Y_1\left(s\right)$ and $Y_2\left(s\right)$, where the first has the same structure of X(s), therefore we can use the same formula of inverse transform and solve the first term:

$$Y_{1\left(s\right)}=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{k_2}\frac{k_2}{s\left(s^{\alpha }+k_2\right)},$$

(29)

applying Laplace transform to $Y_{1\left(s\right)}$:

$$y_1\left(t\right)=\frac{k_2{x}_{d2}-k_1{x}_{d1}}{k_2}\left[1-E_{\alpha }\left(-k_2{t}^{\mathsf{\alpha }}\right)\right],$$

(30)

and for the second term $Y_2\left(s\right),$ we have:

$$Y_{2\left(s\right)}=x_{d1}\frac{k_1^2}{s\left(s^{\alpha }+k_1\right)\left(s^{\alpha }+k_2\right)},$$

(31)

where we multiply the denominator of $Y_{2\left(s\right)}$ by $s^{2\mathsf{\alpha }}/s^{2\mathsf{\alpha }}$:

$$Y_{2\left(s\right)}=x_{d1}\frac{k_1^2}{\frac{s^{2\alpha +1}}{s^{2\alpha }}\left(s^{\alpha }+k_1\right)\left(s^{\alpha }+k_2\right)},$$

(32)

and by rearranging the denominator products, we obtain:

$$Y_{2\left(s\right)}=x_{d1}\frac{k_1^2}{s^{2\alpha +1}\left(1+k_1{s}^{-\alpha }\right)\left(1+k_2{s}^{-\alpha }\right)},$$

(33)

We apply the properties of the geometric series to both products: $1/\left(s^{\alpha }+k_1\right)$ and $1/\left(s^{\alpha }+k_2\right)$:

$$Y_{2\left(s\right)}=x_{d1}{k}_1^2 s^{-2\alpha -1}\sum _{m=0}^{\infty }{\left(-k_1{s}^{-\alpha }\right)}^m\sum _{n=0}^{\infty }{\left(-k_2{s}^{-\alpha }\right)}^n,$$

(34)

Applying the properties of the series product we obtain a more compact format:

$$Y_{2\left(s\right)}=x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n{s}^{-\left(m+n+2\right)\alpha -1}.$$

(35)

We express the argument of the series as a fraction and multiply it by $\mathsf{\Gamma }\left[(m+n+2)\alpha +1\right]/\mathsf{\Gamma }\left[(m+n+2)\alpha +1\right]$:

$$Y_{2\left(s\right)}=x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n}{s^{\left(m+n+2\right)\alpha +1}}\frac{\Gamma \left[\left(m+n+2\right)\alpha +1\right] }{\Gamma \left[\left(m+n+2\right)\alpha +1\right]},$$

(36)

and the products are rearranged in the denominator:

$$Y_{2\left(s\right)}=x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n}{\Gamma \left[\left(m+n+2\right)\alpha +1\right]}\frac{\Gamma \left[\left(m+n+2\right)\alpha +1\right]}{s^{\left(m+n+2\right)\alpha +1}},$$

(37)

Now you can find the inverse transform:

$$y_{2\left(t\right)}=x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n}{\Gamma \left[\left(m+n+2\right)\alpha +1\right]}{t}^{(m+n+2)\alpha },$$

(38)

and the solution of is given by:

$$y\left(t;\alpha \right)=\left(x_{d2}-\frac{x_{d1}{k}_1}{k_2}\right)\left[1-E_{\alpha }\left(-k_2{t}^{\alpha }\right)\right]+x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n}{\Gamma \left[\left(m+n+2\right)\alpha +1\right]}{t}^{\left(m+n+2\right)\alpha }.$$

(39)

4.2. Mathematical Demonstration of Fractional Functions

Demonstrating fractional functions is fundamental to validate the mathematical solidity and formality of the models developed. When considering the case where the alpha value equals 1, the fractional solution is expected to match the solution obtained by traditional methods. This coincidence demonstrates the consistency of the fractional approach with conventional methods and confirms the mathematical validity of the models. Therefore, the accuracy and coherence of the fractional approach in the representation of dynamic systems can be verified by comparing the solutions obtained with and without the use of fractional functions. This comparison is crucial to ensure that the mathematical tools are appropriate and effective in modeling complex phenomena.

The demonstration is performed for $x\left(t;\alpha \right)$, starting from Equation (26) when $\alpha =1$:

Now, we replace it in the Mittag-Leffler function:

$${ E}_{\alpha }\left(z\right) = \sum _{k=0}^{\infty } \frac{z^k}{\Gamma \left(\left(1\right)k + 1\right)},$$

(40)

resulting in the denominator equal to $\Gamma \left((k + 1\right))$=k!. Then, we have the Euler function:

$${ E}_{\alpha }\left(z\right) = \sum _{k=0}^{\infty } \frac{z^k}{k!}=e^k,$$

(41)

replace the value of k, by $k1$ Equation (16) is obtained by performing the demonstration.

Resulting in being equal to the ordinary solution by testing the correct solution by the Caputo method for state variable 1. The demonstration is performed for $y\left(t;\alpha \right)$, when $\alpha =1$:

$$y\left(t;1\right)=\left(x_{d2}-\frac{x_{d1}{k}_1}{k_2}\right)\left[1-E_{\alpha }\left(-k_2{t}^1\right)\right]+x_{d1}{k}_1^2\sum _{m,n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}{k}_1^m{k}_2^n}{\Gamma \left[\left(m+n+2\right)+1\right]}{t}^{\left(m+n+2\right)}.$$

(42)

To demonstrate the function found with the fractional calculation when $\alpha =1$, we apply the Laplace transform, which divides the transformation into two parts:

$$Y_{1\left(s\right)}=\left(x_{d2}-\frac{x_{d1}{k}_1}{k_2}\right)\left(\frac{k_2}{s\left(s+k_2\right)}\right),$$

(43)

The partial fractions solve the second term left in the equation as follows:

$${ Y}_{1\left(s\right)}=\left(x_{d2}-\frac{x_{d1}{k}_1}{k_2}\right)\left(\frac{1}{s}-\frac{1}{s+k_2}\right),$$

(44)

Later, the transform of the place is applied to the second part:

$$Y_{2\left(s\right)}=x_{d1}{k_1}^2\left(\sum _{m=0}^{\infty } \frac{{\left(-1\right)}^m{k_1}^m}{s^{m+1}}\right)\left(\sum _{n=0}^{\infty } \frac{{\left(-1\right)}^m{k_2}^m}{s^{m+1}}\right)s^{-1},$$

(45)

where the denominator passes to multiply the numerator in each of the summations.

$$Y_{2\left(s\right)}=x_{d1}{k_1}^2\left(\sum _{m=0}^{\infty } {\left(-1\right)}^m{k_1}^m{s}^{-m}\right)\left(\sum _{n=0}^{\infty } {\left(-1\right)}^m{k_2}^m{s}^{-m}\right)s^{-3},$$

(46)

The terms are arranged to obtain the geometric series in the two terms:

$${ Y}_{2\left(s\right)}=x_{d1}{k_1}^2\left(\sum _{m=0}^{\infty } {\left({-k}_1{s}^{-1}\right)}^m\right)\left(\sum _{n=0}^{\infty } {\left({-k}_2{s}^{-1}\right)}^m\right)s^{-3},$$

(47)

The approximation of the geometric series obtained is represented as follows:

$${ Y}_{2\left(s\right)}=x_{d1}{k_1}^2\left(\frac{1}{1+k_1{s}^{-1}}\right)\left(\frac{1}{1+k_2{s}^{-1}}\right)s^{-3},$$

(48)

The terms are arranged as follows:

$${ Y}_{2\left(s\right)}=x_{d1}{k_1}^2.\frac{1}{s^3}\left(\frac{1}{1+k_1{s}^{-1}}\right)\left(\frac{1}{1+k_2{s}^{-1}}\right),$$

(49)

Multiply and rearrange the denominator:

$${ Y}_{2\left(s\right)}=x_{d1}{k_1}^2\frac{1}{s(s+k_1)(s+k_2)},$$

(50)

The linear partial fractions are made:

$$Y_{2\left(s\right)}=x_{d1}{k_1}^2\frac{1}{s(s+k_1)(s+k_2)}=\frac{A}{s}+\frac{B}{s+k_1}+\frac{D}{s+k_2},$$

(51)

Resulting in the partial fractions:

$$Y_{2\left(s\right)}=x_{d1}{k_1}^2\left(\frac{\frac{1}{k_1{k}_2}}{s}-\frac{\frac{1}{k_1{(k}_2-k_1)}}{s+k_1}+\frac{\frac{1}{k_2{(k}_2-k_1)}}{s+k_2}\right),$$

(52)

multiplying and simplifying the equation.

$$Y_{2\left(s\right)}=\frac{\frac{x_{d1}{k}_1}{k_1{k}_2}}{s}-\frac{\frac{x_{d1}{k}_1}{k_1{(k}_2-k_1)}}{s+k_1}+\frac{\frac{x_{d1}{k_1}^2}{k_2{(k}_2-k_1)}}{s+k_2},$$

(53)

are added $Y_{1\left(s\right)}+Y_{2\left(s\right)}$

$$Y\left(s\right)=\frac{x_{d2}}{s}-\frac{x_{d2}}{s+k_2}-\frac{x_{d1}{k}_1}{k_{2}s}+\frac{x_{d1}\frac{k_1}{k_2}}{s+k_2}+\frac{\frac{x_{d1}{k}_1}{k_1{k}_2}}{s}-\frac{\frac{x_{d1}{k}_1}{k_1{(k}_2-k_1)}}{s+k_1}+\frac{\frac{x_{d1}{k_1}^2}{k_2{(k}_2-k_1)}}{s+k_2},$$

(54)

reducing similar terms and drawing common ground $x_{d2}$.

$$Y\left(s\right)=x_{d2}\left(\frac{1}{s}-\frac{1}{s+k_2}\right)+\frac{\frac{x_{d1}{k_1}^2}{k_2{(k}_2-k_1)}}{s+k_2}+\frac{x_{d1}\frac{k_1}{k_2}}{s+k_2}-\frac{\frac{x_{d1}{k}_1}{k_1{(k}_2-k_1)}}{s+k_1},$$

(55)

$$Y\left(s\right)=x_{d2}\left(\frac{1}{s}-\frac{1}{s+k_2}\right)+\frac{\frac{x_{d1}{k_1}^2}{k_2{(k}_2-k_1)}}{s+k_2}+\frac{x_{d1}\frac{k_1{(k}_2-k_1)}{k_2{(k}_2-k_1)}}{s+k_2}-\frac{\frac{x_{d1}{k}_1}{k_1{(k}_2-k_1)}}{s+k_1},$$

(56)

reducing similar terms in $Y\left(s\right)$.

$$Y\left(s\right)=x_{d2}\left(\frac{1}{s}-\frac{1}{s+k_2}\right)+\frac{\frac{x_{d1}{k}_1}{{(k}_2-k_1)}}{s+k_2}-\frac{\frac{x_{d1}{k}_1}{{(k}_2-k_1)}}{s+k_1},$$

(57)

the common factor is removed in $\frac{x_{d1}{k}_1}{{(k}_2-k_1)}$ the equation.

$$Y\left(s\right)=x_{d2}\left(\frac{1}{s}-\frac{1}{s+k_2}\right)+\frac{x_{d1}{k}_1}{{(k}_2-k_1)}\left(\frac{1}{s+k_2}-\frac{1}{s+k_1}\right),$$

(58)

Now, you can find the inverse transform in Equation (17):

Therefore, by checking that the solution obtained by the Caputo method of the fractional calculation, by substituting α = 1, which coincides with the function in ordinary time, we validate the fractional solution. This result confirms the consistency and precision of the fractional approach in the representation of dynamic systems. By evaluating both fractional functions in 1 and obtaining the function over time for the two state variables, the validity and reliability of the models developed are corroborated. This finding reinforces confidence in using fractional calculus as a practical mathematical tool to address complex dynamic phenomena.

4.3. Simulating Diverse Scenarios

In this section, we will explore several scenarios, varying both production goals and alpha values in each department, to better understand how these variables affect the performance of the Weaving and Basting process. This will provide valuable information for decision-making in the management of textile production.

4.3.1. Simulation of Weaving Department Efficiencies

It is important to note that the alpha parameter plays a fundamental role as an indicator of operational efficiency, and its range of variation allows a thorough assessment of the efficiency of the department in question.

Table 9. Parameter variance for different simulations of the Weaving department.

Scenarios	Category	Target (Pieces)	Efficiency
Scenario 1	Optimal	100	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$
Scenario 2	Optimal	70	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$
Scenario 3	Optimal	50	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$

shows the different scenarios for the Weaving department.

Figure 5 shows simulation 1, which compares the six efficiency categories with 100% performance.

Figure 6 shows simulation 2, which compares the six efficiency categories with 100% performance.

Figure 7 shows simulation 3, which compares the six efficiency categories with 100% performance.

4.3.2. Simulation of Basting Department Efficiencies

Table 10. Parameter variance for different simulations of the Basting department.

Scenarios	Category	Target (Pieces)	Efficiency
Scenario 1	Optimal	100	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$
Scenario 2	Optimal	70	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$
Scenario 3	Optimal	50	$\alpha =0.9$
	Very good		$\alpha =0.8$
	Good		$\alpha =0.7$
	Fair		$\alpha =0.6$
	Poor		$\alpha =0.5$
	Very poor		$\alpha =0.4$

details the different scenarios considered for the Basting department. Each scenario represents a unique combination of variables, including the department’s production goal and the efficiency values associated with the processes, expressed in alpha terms. These scenarios provide a comprehensive view of the possible configurations under which the performance of the Basting department in the manufacturing process is evaluated.

Figure 8 shows simulation 1, which compares the six efficiency categories with 100% performance.

Figure 9 shows simulation 2, which compares the six efficiency categories with 100% performance.

Figure 10 shows the third scenario for the Basting department, where the six comparative labels are shown.

Simulations provide us with an invaluable tool for decision-making. By exploring a variety of scenarios that encompass different efficiency values for both departments, we can gain a deeper understanding of how these factors affect the production process. This ability to simulate and compare multiple scenarios gives us valuable insight to optimize our operations and make informed decisions that drive efficiency and productivity across the organization.

4.4. Field Validation

In this section, we contrast the efficiency of 100% with the efficiencies observed in each of the two departments validated in the field. This meticulous analysis allows us to see the operational reality of both departments. In addition, it allows us to identify significant discrepancies and areas of potential improvement, which contributes to more informed decision-making and continuous optimization of our processes.

4.4.1. Simulation of Weaving Department Efficiencies

The results of field validation using the OEE methodology are presented in Figure 11 and detailed in the third section. In this context, the efficiency obtained for the Weaving department was 77.90%.

4.4.2. Simulation of Basting Department Efficiencies

The results of field validation using the OEE methodology are presented in Figure 12 and detailed in the third section. In this context, the efficiency obtained for the Basting department was 90.84%. It can be concluded that the Weaving process is in the category “Good” since it exceeds 70% of estimated efficiency and has a wide margin for improvements in the two categories. It is estimated that it could achieve approximately 50 pieces. On the other hand, in the Basting process, it is in the “Optimum” category, exceeding 90% efficiency, which indicates a very good performance in this department, obtaining 67 pieces during the given period.

The results presented in Section 4.4 are specifically designed and analyzed for domain 0 < α < 1. This range was selected intentionally because it makes sense for the application. Since the value of α does not measure absolute production but production efficiency, it is limited to the range from 0 to 1. In this context, a value of α greater than 1 could indicate an efficiency greater than 100%, which would not be practical or relevant to interpreting efficiency in the process. Therefore, we focused on the range 0 < α < 1, where the results are most significant and applicable for this study.

5. Conclusions

This research project presents an innovative proposal by introducing fractional calculus into system dynamics models, explicitly focusing on equilibrium loops. By addressing discrepancies between ideal simulations and actual models, this initiative seeks to provide a more realistic understanding of systems’ behavior, considering their different levels of efficiency. The integration of fractional calculus significantly strengthens existing models, providing more accurate and helpful information for decision-making in various contexts. This contribution represents a valuable advance in the field of system dynamics in line with the legacy of engineer Jay Forrester.

The results obtained from the simulations in Figure 3 and Figure 4 reveal a series of significant findings. First, the validity of the solutions obtained by comparing ordinary differential equations with fractional solutions for different alpha values has been confirmed. This finding reinforces the solidity of the approach used in the study. In addition, the analysis of efficiency curves in the different areas of the company determines the completion time of orders based on the percentage OEE indicator. For example, in the Tissue department, with 91% efficiency. The proposed approach determines the respective duration time. This efficiency level highlights the importance of implementing continuous improvement strategies in production management to meet deliveries promptly. In short, the results provide a solid foundation for informed decision-making to diagnose company processes and improve overall performance. Future work will analyze the efficiency behavior of an industrial process using the theory of conformable calculus.