A New Linear Two-State Dynamical Model for Athletic Performance Prediction in Elite-Level Soccer Players

1. Introduction

In recent decades, mathematical modeling has become a fundamental tool in sports science, allowing sports and data scientists to better understand the close relationship between training load and athletes’ performance capacity [1]. The data obtained from such modeling have provided coaches with valuable support for the prescription of increasingly effective conditioning programs aimed at maximizing athletes’ performance capacity at key moments, while avoiding illness, injury or overtraining.

The origins of these mathematical models date back to the 1970s through the studies of Banister, whose fitness-fatigue model represented the first attempt to formalize the relationship between training stimulus and performance [1]. The model proposed by Banister is based on the assumption that athletes’ performance capacity is the result of two competing factors: a first factor, called fitness, represents the positive effects of training on performance, while a second factor, called fatigue, describes the negative effects that reduce performance capacity. From a practical point of view, the model predicts that after a workload there is a brief decrease in performance due to accumulated fatigue, followed after adequate recovery by an improvement in physical capacity thanks to the adaptations induced by training [1].

Subsequent to the publication of this model, numerous studies have confirmed its effectiveness in predicting the performance capacity of athletes in several individual sports. In fact, this model has been applied in disciplines such as running [2], cycling [3], swimming [1,4], demonstrating in many of these studies a notable accuracy in predicting the actual performance of athletes [1,3,4,5].

More recently, research conducted by Busso has allowed to further develop and refine Banister’s model, successfully introducing new parameters specifically designed to capture the effects of training monotony—defined as low variability in daily workloads—and the cumulative fatigue that arises when intense sessions are performed consecutively without adequate recovery. These additions enhance the model ability to simulate real-world training responses [6]. Although these models have demonstrated their validity in many individual sports, they have proven to be less accurate when applied to athletes in team sports such as soccer, where performance is influenced by multiple physiological factors [7,8]. Indeed, a soccer match requires players to undergo numerous phases of intermittent exercise characterized by periods of high intensity alternating with phases of lower intensity or recovery [9]. This highly variable metabolic demand, to which players are subjected, represents a complex challenge for traditional models, which are often based on the assumption of a linear relationship between training load and performance.

Furthermore, soccer players are often engaged in multiple daily training sessions, each characterized by different intensities and physiological objectives (e.g. strength and power training in the morning and technical/tactical training in the afternoon). The variability of the physiological objectives of the training sessions can produce cumulative effects on the performance capacity of athletes, which are difficult to be predicted [2]. Recent studies have increasingly focused on developing predictive models of athletic performance in soccer, adopting approaches such as time-series analysis, machine learning techniques, and piecewise-linear formulations to improve the accuracy of performance forecasts [10,11,12]. Nevertheless, these approaches present relevant limitations such as reliance on subjective self-reports instead of GPS-derived external load [10], a focus on short-window extreme events with small sample sizes and acknowledged reliability constraints [11], or consideration of only the last training sessions without an explicit dynamical representation of training–performance interactions [12].

Consequently, building effective mathematical models in a sport such as soccer requires the inclusion of multiple and complex variables that reflect the athlete’s performance ability in a dynamic context typical of game actions. Recent advances in wearable technology have provided sports scientists with tools able to acquire more detailed data on athletes’ performance both during training sessions and during competitions [13,14]. Devices such as GPS systems and accelerometers enable real-time monitoring of critical performance parameters, including speed, acceleration, metabolic power and workload [15].

These technologies have generated enormous amounts of data, providing new opportunities to refine existing models or create more sophisticated ones. For example, some studies have shown that soccer performance is more accurately predicted when these parameters are combined together to calculate training loads, which take into account the intermittent nature of the sport [15,16].

Based on these considerations, the original contribution of this paper is to demonstrate that a new linear, time-varying, two-state model—linking the dynamics of a performance-related state variable with those of a lumped state variable, and whose parameters are estimated via the Differential Evolution optimization method—achieves more accurate performance predictions than the (nonlinear) Busso model. This improvement is made possible by integrating a wider range of performance metrics derived from GPS technology. The proposed model provides a more comprehensive evaluation of the interplay between training and recovery, offering a more reliable tool for predicting performance in team sports such as soccer.

Importantly, by modeling the effects of training load on future performance, the proposed modeling approach can also support the early identification of excessive or poorly distributed training sessions, helping coaches optimize training loads and prevent maladaptive responses and reduce the risk of injury. This makes it not only a predictive model, but also a preventive tool for performance sustainability and long-term athlete’s health.

2. Materials and Methods

2.1. Subjects

Two male professional soccer players were recruited from the Italian Serie A League in order to monitor their readiness during a competitive season. Players were active members of a Serie A team and had at least 5 years of competitive experience in the top division.

2.2. Video Match Analysis

A Match analysis was conducted using the STATS SportVU® system, a validated multicamera video analysis tool from STATS LLC in Chicago, United States. The STATS SportVU® system is a sophisticated and widely used technology in sports analytics, especially soccer. It utilizes multiple high-definition cameras installed at various angles around the stadium to track players and the ball at rates up to 25 Hz, with a measurement accuracy determined by the Technical University of Munich (TUM) to have a typical error of 2.7% for total distance [17]. K-Sport provided raw data in Cartesian coordinates, and the primary data were smoothed at 5 Hz [18]. STATS SportVU® uses high-definition cameras, advanced software, and statistical algorithms to extract and process player (X, Y) and ball (X, Y, Z) coordinates [19]. Roof-level cameras captured player movements during matches, and the data were analyzed using STATS Viewer and STATS Dynamix, both provided by STATS Perform, to create a dataset of each player’s physical and technical performance. A specific software (https://www.gpexe.com) was used to download and analyze data extracted from STATS SportVU, including metabolic power (W·${\text{kg}}^{-1}$), equivalent distance index (%), and other metabolic-energetic metrics. Data were filtered according to the theoretical model based on an energetic approach where the energy cost of accelerations and decelerations plays a central role, as suggested by di Prampero et al. [20].

2.3. GPS Technology

The external load during each training session was evaluated using high-frequency GPS Pro2 devices operating at 18 Hz (GPEXE©, Exelio srl, Udine, Italy). The reliability, accuracy, and validity of the GPEXE© system have been reported in previous studies [21,22]. This level of resolution allows for a detailed assessment of players’ running patterns, directional changes, and speed profiles during both training sessions and official matches. The technology provides coaches with real-time performance metrics, supporting data-driven decision-making. The system is capable of calculating over 300 metrics in real-time and ensures rapid and accurate data transfer. Additionally, it supports three-dimensional tracking of individual players or teams over time, enabling assessments of positional workloads, training intensities, and physiological demand changes. All GPS devices were activated 15 minutes prior to data collection to ensure optimal signal quality and satellite lock. Signal quality was evaluated based on the number of connected satellites (11 ± 1.4) and the horizontal dilution of precision (0.9 ± 0.1), according to the criteria established by Witte and Wilson [23].

2.4. Parameter Estimation via Differential Evolution (DE)

To identify the optimal parameter values for both the performance prediction models (namely, the proposed two-state model and the comparative Busso model), we employed the Differential Evolution (DE) optimization algorithm. DE is particularly well-suited for solving constrained nonlinear optimization problems and has been implemented in Python to automate the parameter estimation process.

The DE algorithm operates by maintaining a population of candidate solutions and iteratively combining them through mutation, crossover, and selection mechanisms to explore the search space. This evolutionary approach ensures a global search capability while progressively minimizing the objective function, defined here as the Sum of Squared Errors (SSE) between the model-predicted performance and the actual observed performance data.

In our implementation, the SSE cost function was evaluated over the entire time series, and the algorithm was initialized with random values within the physiologically meaningful bounds subsequently reported in Table 1. The optimization was repeated multiple times (with different seeds) to avoid convergence to local minima and ensure robustness of the estimated parameters.

The DE-based estimation procedure was applied individually for each subject and for both models, allowing direct comparison of model fidelity across different formulations.

This parameter identification process ensures that each model configuration is optimally fitted to the athlete’s data, enabling reliable assessment of predictive accuracy and physiological interpretability.

2.5. Study Design

This study aimed at evaluating a novel mathematical model for predicting performance in elite-standard soccer players. Specifically, the study focuses on players participating in the Italian Premier Division (Serie A) during the 2018–2019 season. They were monitored over half a season (6 months), including the pre-season period, the first part of the regular season, and all official matches (Serie A and National League matches), which took place from July to December during the 2018–2019 season. During the pre-season, players underwent training sessions seven times a week, with friendly matches typically scheduled on Thursdays or over the weekend. These training sessions mainly focused on developing technical-tactical skills and physical fitness. During the league season, players trained six times a week, with an official match played on the weekend. Friendly and National Cup matches were usually played on Thursdays and Wednesdays, respectively.

2.6. Quantification of Training

Training load quantification was performed using the equation proposed by [16], which combines metabolic power (MP) and equivalent distance index (EDI) to assess external training load (ETL). Notably, this method demonstrated a significant correlation with internal training load (ITL), as measured using Edwards’ heart rate-based TL method. This finding supports the validity of the MP-based training load calculation in effectively reflecting internal physiological demands during soccer activities. Briefly, metabolic power metrics were developed to accurately describe the physical demands placed on soccer players during matches or training sessions [15]. Integrating metabolic power with kinematic parameters provides an integrated measurement of acceleration and velocity, making it a more precise indicator of the physiological effort exerted by soccer players during both training and match events [15]. In particular, to calculate metabolic power (MP), researchers applied the following formula [20]:

$$\mathrm{MP}=\mathrm{EC}·\mathrm{V}$$

(1)

where:

• $\mathrm{EC}$: energetic cost of running uphill at an inclination derived from the running acceleration vector;
• $\mathrm{V}$: instantaneous velocity.

To better assess the metabolic intensity of the training or game, some authors propose incorporating both the equivalent distance (ED) and the equivalent distance index (EDI) [16,24,25]. ED refers to the distance an athlete would have covered at a constant pace on grass using the total energy expended during the match, while EDI is the ratio between the ED and the total distance covered in the game. It’s important to note that the EDI provides an indirect estimate of the metabolic intensity of the match, with a high index indicating a greater contribution from anaerobic energy stores due to high-intensity acceleration efforts. To calculate ED and EDI, the following formulas were applied [15]:

$$\mathrm{ED}={\displaystyle \frac{\mathrm{W}}{{\mathrm{EC}}_c·\mathrm{KT}}},\mathrm{EDI}={\displaystyle \frac{\mathrm{ED}}{\mathrm{TD}}}$$

(2)

where:

• $\mathrm{W}$: total energy expenditure measured in joules per kilogram;
• ${\mathrm{EC}}_c$: energetic cost of running at a constant pace on flat, compact terrain (assumed to be $3.6\mathrm{J}·{\mathrm{kg}}^{-1}·{\mathrm{m}}^{-1}$);
• $\mathrm{KT}$: grassy terrain constant;
• $\mathrm{ED}$: equivalent distance measured in meters;
• $\mathrm{TD}$: total distance covered in meters.

Based on existing data, the average ratio between ED and total distance covered during soccer matches was approximately 1.20 [15], though this value varied among players (range: 1.15 to 1.35) depending on positional role [26]. In our study, we calculated the external training load (ETL) on the basis of the product of MP and EDI [16]. It defines the intensity of the external load (I) as follows:

$$\mathrm{ETL}=\int \mathrm{V}\left(t\right)dt\times \int \mathrm{I}\left(t\right)dt,\mathrm{I}=\mathrm{MP}\times \mathrm{EDI}$$

(3)

2.7. Performance Test

All performance data were obtained from official matches and collected through video match analysis. As a performance criterion, we used the ratio between the total distance covered by the athlete and the distance covered at a metabolic power output exceeding 25 W·${\text{kg}}^{-1}$ (indicative of high-intensity activity). This threshold effectively identifies efforts with elevated metabolic demands and is considered a reliable indicator for assessing high-intensity performance in soccer players [27]. For the analysis, only players who completed the entire match or were substituted no earlier than five minutes before the final whistle (i.e., ≥85 minutes of play) were included.

3. Modeling

The modeling approach employed in this study aims at originally describing and predicting athletes’ performance capacity in response to ETL derived from GPS-based metrics. We developed a novel linear, time-varying, two-state, discrete-time dynamical model that draws inspiration from the model of heart rate dynamics during variable-intensity exercises [3,4,20].

The proposed model includes two state variables: $x_1\left(k\right)$, representing a performance-related variable at time step k, and $x_2\left(k\right)$, a lumped variable at time step k whose variations directly affect the temporal dynamics of $x_1\left(k\right)$. The model is governed by the following recursive equations:

$$\left\{\begin{array}{c}{x}_1(k+1)=a_1{x}_1\left(k\right)-a_2{x}_2\left(k\right)e^{-k/{\tau }_1}+b_1{e}^{-k/{\tau }_2}u\left(k\right)\hfill \\ x_2(k+1)=a_3{x}_2\left(k\right)+b_2\left(u\left(k\right)-a_4{e}^{-k/{\tau }_3}{x}_1\left(k\right)\right)\hfill \end{array}\right.$$

(4)

where $u\left(k\right)$ is the (%-normalized) training load input, and $k\in \mathbb{N}$. The model incorporates multiple time constants ${\tau }_i$ to account for decaying memory effects and parameters $a_1$, $a_2$, $a_3$, $a_4$ to fit the ($x1\left(k\right)$ , $x_2\left(k\right)$) profiles satisfactorily. To ensure physiological plausibility and numerical stability, we introduced weak bounding constraints on all model parameters. These bounds were based on prior modeling studies [1,2,5]. Table 1 summarizes the description of the model parameters and the bounds applied during model optimization. Note that when $b_1$ is negative, the load $u\left(k\right)$ directly contributes negatively to $x_1(k+1)$, whereas it may still contribute positively to $x_1(k+2)$, in accordance with the following auto-regressive moving-average expression1:

$$\begin{array}{cc}\hfill x_1(k+2)& =\left(a_1+{\mathrm{e}}^{-1/{\tau }_1}{a}_3\right)x_1(k+1)-{\mathrm{\Sigma }}_1{x}_1\left(k\right)\hfill \\ \hfill & +b_1{\mathrm{e}}^{-(k+1)/{\tau }_2}u(k+1)-{\mathrm{\Sigma }}_{2}u\left(k\right)\hfill \end{array}$$

(5)

with

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_1& =a_3{a}_1{\mathrm{e}}^{-1/{\tau }_1}-a_2{a}_4{b}_2{\mathrm{e}}^{-(k+1)/{\tau }_1}{\mathrm{e}}^{-k/{\tau }_3}\hfill \\ \hfill {\mathrm{\Sigma }}_2& =a_3{b}_1{\mathrm{e}}^{-1/{\tau }_1}{\mathrm{e}}^{-k/{\tau }_2}+a_2{b}_2{\mathrm{e}}^{-(k+1)/{\tau }_1}\hfill \end{array}$$

(6)

This formulation allows for a more biologically plausible representation of the complex relationship between training and performance in soccer players.

On the other hand, the input term $b_1{e}^{-k/{\tau }_2}u\left(k\right)$ introduces a form of temporal “memory” in the system: past training loads are exponentially attenuated over time, meaning that, a training session at time k will progressively lose its influence on performance estimation at time k + 1 as k increases. This mechanism reflects the fading physiological impact of past stimuli and allows the model to mimic the adaptive dynamics of athletic response.

From a mathematical standpoint, the proposed formulation also marks a structural departure from prior models. Whereas performance was previously estimated via single or double summations of exponentially decaying terms related to training and fatigue, here it emerges from a linear, time-varying, two-state system, with a lumped variable. This architecture enhances both the flexibility and interpretability of the model under dynamic load conditions.

Moreover, the model generates an asynchronous output – related to the performance achieved in match events – so that the visible effects of training inputs are not reflected in the instantaneous performance. This behavior captures the delayed and distributed nature of physiological adaptations and fatigue-recovery cycles, thus aligning the model behavior more closely with real-world athletes’ response.

More specifically, the model exhibits asynchronous behavior by distributing the influence of each training input over multiple future time steps. This is governed by a set of estimated time constants that rule how long the effects of a training session persist in the athlete’s physiological state. As a result, a single high-load session may influence performance across several subsequent days—initially increasing fatigue, and only later contributing to fitness improvements. This mechanism allows the model to capture both the delayed onset of adaptation and the extended recovery period, offering a more biologically realistic depiction of training responses.

To evaluate its effectiveness and highlight the value of its novel features, the proposed two-state model was systematically compared with a well-established reference, namely the three-component (Busso) model [6], which introduces a dynamic fatigue term accounting for the monotony and cumulative load of prior sessions.

Researchers started from the one-component model, which represents the most basic approach, assuming that performance is solely driven by the accumulation of fitness as a function of prior training loads [1]. This model neglects fatigue effects and uses a single exponential decay term to reflect the residual impact of each training session. Despite its simplicity, it provides a useful baseline for evaluating the predictive power of more complex models.

In general, a fitness-fatigue model conceptualizes the performance as the sum of a positive training effect (fitness) and a negative one (fatigue). It estimates the performance through the convolution of exponentially decaying responses to training impulses, weighted by individual-specific coefficients and time constants [1,2].

The mathematical formulation of the comparative model used in this study is reported below for the sake of exhaustiveness:

• Three-component (Busso) model [28]:

  $$p\left(k\right)=p\left(0\right)+k_1\sum _{i=1}^{k-1}w\left(i\right)e^{-{\displaystyle \frac{(k-i)}{{\tau }_1}}}+\sum _{i=1}^{k-1}{k}_{2}w\left(i\right)e^{-{\displaystyle \frac{(k-i)}{{\tau }_2}}}+\sum _{i=1}^{k-1}{k}_2\left(i\right)w\left(i\right)e^{-{\displaystyle \frac{(k-i)}{{\tau }_2}}}$$

  (7)

  $$k_2\left(i\right)=k_3\sum _{j=1}^{i}w\left(j\right)e^{-{\displaystyle \frac{(i-j)}{{\tau }_3}}}$$

  (8)
  where $k_2\left(i\right)$ is a dynamic fatigue term representing the cumulative training monotony, and $k_3$ and ${\tau }_3$ model its influence.

4. Results

4.1. Athletic Performance Prediction

As previously mentioned, the predictive capacity of the implemented models was evaluated by analyzing data from two professional soccer players.

Athletic performance was assessed, for both Player 1 and Player 2, as the ratio between the total distance covered and the distance covered at a metabolic power greater than 25 W·${\text{kg}}^{-1}$, which was considered indicative of high-intensity activity. Data were obtained via video analysis of official matches and high-frequency GPS devices (GPEXE©).

The predictions were made for each player at 14 official competitive matches, held on the following dates: 2, 15, 19, and 22 September; 1, 22, and 28 October; 4 and 25 November; 1, 8, 16, 22, and 26 December 2018. Furthermore, for Player 2, performance prediction was also assessed during the pre-season phase, considering 8 friendly matches held on 14, 21, 25, and 28 July; 1, 4, 12, and 26 August 2018.

4.1.1. Model Accuracy Comparison Across Seasonal Phases

The results reported in this section show that the proposed model significantly outperforms the Busso model in predictive accuracy. Specifically, for Player 1, the proposed model achieved a Pearson correlation coefficient of $r=0.89$ and a Sum of Squared Errors ($\mathrm{SSE}=5.84$). In contrast, the Busso model yielded lower correlations and higher errors ($r=0.75$, $\mathrm{SSE}=11.95$).

The improved predictive performance of the proposed two-state dynamical model compared to the Busso model is consistent across the two athletes analyzed (Table 2), confirming that the inclusion of a lumped variable allows for a more accurate modeling of the complex relationship between training load and performance in professional soccer. Estimated model parameters for each player are reported in Table 3.

Figure 1a and Figure 1b show the performance predictions for Player 1 during the competitive season, obtained using the newly proposed two-state model and the Busso model, respectively. The corresponding predictions for Player 2 are presented in Figure 2a and Figure 2b. Furthermore, as shown in Figure 3a-b and in Table 4, the proposed model also achieves superior predictive accuracy during the pre-season (Player 2). The corresponding parameter estimates are provided in Table 5.

These findings confirm the generalizability and robustness of the proposed model across different athletes and seasonal contexts, including both the competitive and pre-season phases.

The following initial conditions at time step $k=0$ were adopted for the two-state variables $x_1\left(k\right)$ and $x_2\left(k\right)$ in the newly proposed model: $x_1\left(0\right)=31.28$, $x_2\left(0\right)=2.78$ for Player 1 (competitive season), $x_1\left(0\right)=27.98$, $x_2\left(0\right)=3.18$ for Player 2 (competitive season), and $x_1\left(0\right)=14.44$, $x_2\left(0\right)=6.75$ for Player 2 (pre-season). In all cases, $x_1\left(0\right)$ was set equal to the measured performance at $k=0$, while $x_2\left(0\right)$ was involved in the DE optimization procedure until convergence to the global minimum.

5. Discussion

This study evaluated the performance of a novel linear, time-varying, two-state discrete-time mathematical model designed to predict athletic performance in professional soccer players. The proposed approach outperformed the (nonlinear) Busso model, showing higher Pearson correlation coefficients and lower SSE values across different athletes and seasonal phases, both during the competitive and pre-season periods (Table 2, Table 3 and Table 4; Figure 1a-b, Figure 2a-b, Figure 3a-b).

Building on this structural insight, the model also introduces a predictive-preventive perspective. In fact, beyond simply describing performance retrospectively, the model enables early identification of suboptimal load trends. This capability supports timely intervention and may help mitigate injury risks. Such anticipatory modeling is of particular interest to training staff who aim to monitor and optimize load progression proactively.

To ensure consistency with biological plausibility, particular care was taken during the model-fitting phase. Specifically, parameter estimation was constrained within physiologically meaningful ranges. This approach enhances both model robustness and transferability, especially when applied to athletes with limited data availability.

On the other hand, the integration of GPS-derived external load metrics—specifically, MP and EDI—further increased the model responsiveness. These metrics allowed a more sensitive tracking of high-intensity activities and training density, which are often underrepresented by traditional indicators. Their inclusion highlights the importance of selecting input variables that reflect the complexity of the training stimulus.

The obtained results suggest that the proposed model offers a biologically informed, data-driven framework for monitoring and optimizing performance in elite soccer. Its structure allows for adaptive modeling under dynamic training regimes, and its parameters provide interpretable insights that align with established physiological mechanisms.

Additionally, implementing real-time model updating—via sliding-window mechanisms—may support dynamic performance tracking throughout entire competitive calendars.

6. Limitations

This paper presents a proof of concept, based on a relatively small sample of two professional players. As such, the observed differences in predictive accuracy between models should be interpreted with caution. Future studies will involve larger and more diverse athlete populations.

7. Conclusions

This study introduced a novel linear, time-varying, two-state discrete-time mathematical model that outperforms a classical model in predicting athletic performance among elite soccer players.

A key innovation is the inclusion of GPS-derived metrics such as metabolic power and equivalent distance index, which enhanced the model sensitivity to the physiological intensity of training sessions.

The results obtained in this paper suggest that the proposed model is a powerful tool for performance prediction while suggesting that it is also a practical framework for guiding training load optimization in high-intensity team sports. It offers a biologically grounded and data-driven approach to enhance athlete readiness and reduce the risk of overtraining. Moreover, by predicting how training load impacts future performance, the model provides valuable support in injury prevention strategies, helping to avoid maladaptive responses due to training errors or excessive load accumulation.

Future extensions of this work may include real-time applications, integration with sliding window updates, and broader validation across diverse sports and athletic populations.