Data-Driven Framework for Uncovering Hidden Control Strategies in Evolutionary Analysis

We have devised a data-driven framework for uncovering hidden control strategies used by an evolutionary system described by an evolutionary probability distribution. This innovative framework enables deciphering of the concealed mechanisms that contribute to the progression or mitigation of such situations as the spread of COVID-19. Novel algorithms are used to estimate the optimal control in tandem with the parameters for evolution in general dynamical systems, thereby extending the concept of model predictive control. This is a significant departure from conventional control methods, which require knowledge of the system to manipulate its evolution and of the controller's strategy or parameters. We used a generalized additive model, supplemented by extensive statistical testing, to identify a set of predictor covariates closely linked to the control. Using real-world COVID-19 data, we successfully delineated the descriptive behaviors of the COVID-19 epidemics in five prefectures in Japan and nine countries. We compared these nine countries and grouped them on the basis of shared profiles, providing valuable insights into their pandemic responses. Our findings underscore the potential of our framework as a powerful tool for understanding and managing complex evolutionary processes.


Introduction
The evolution of natural dynamical systems is a complex process where the dynamics are influenced by a multitude of factors, including epidemic evolution, environmental changes, and various exogenous random events.Traditional methods of studying evolutionary processes often involve creating simplified models that may not fully capture this complexity.When it comes to influencing or controlling this evolution, optimal control theory (OCT) is a powerful and flexible tool Grüne et al. (2020); Bussell et al. (2019); Ros et al. (2004).This can provide valuable insights into the likely paths of evolution and the factors that may affect it Shen et al. (2018).The aim of this work was to investigate an evolutionary process, such as epidemic evolution, from the perspective of OCT in order to uncover the hidden mechanisms that lead to deterioration or improvement of the situation.We assumed that the evolutionary process is described by a controlled system via a hidden mechanism to which we do not have access and which we cannot change directly; we can only observe its effects through measurements related to its evolution.We also assumed that this hidden Our goal in this work was to make a simple assumption about the dynamics of COVID-19 virus transmission and then utilize data to determine the parameters governing it.An extensive initial analysis of Japanese data was performed to understand the dynamics of the epidemic in diverse prefectures (Tokyo, Osaka, Hokkaido, Fukuoka, and Okinawa).Their differing characteristics were compared and linked via a generalized additive model (GAM) to different sociodemographic structures and the related policies that were adopted in each prefecture.
After validating our data-driven model on the Japanese data, we examined whether it can be applied to other countries or geographical regions with different population structures and dynamics.We did this by applying an MPC-like procedure, as we did for the Japan data, to simultaneously obtain the dynamics and control parameters that contributed to the observed evolution of the epidemic in each country.We then performed a comparative analysis between the countries based on their management strategies, such as for vaccinations and lockdowns.
• Descriptive behavior of epidemic by country: Using the epidemic's fitted controlling parameters for each country or region, we attempted to clarify the relationship between the parameters and various socio-demographic indicators as well as the various policies that governments implemented to reduce the epidemic's spread.
• Comparative analysis and clustering: Using the contrasting characteristics and indicators of the various countries, we grouped them on the basis of the similarities of their profiles.We then looked at the key variables that led to countries being grouped together.We also conducted additional investigations for each group to determine which factors contributed most to the grouping.
Our contributions are summarized below: • We propose using a hidden controller framework in which the optimal control and the parameters for evolution in general dynamical systems including the COVID-19 epidemic dynamical system are simultaneously estimated.This transforms the linear evolution scheme into a nonlinear problem.
• We present an efficient optimization algorithm from the practical point of view that is based on a classical optimization algorithm, e.g., fmincon in MATLAB.
• We present the results of a comparative analysis between prefectures and regions of Japan with differing population structures and epidemic management policies that were obtained using our hidden controller framework.
• We show that generalizing the model using our hidden controller framework to other countries (Australia, Germany, Italy, Brazil, South Africa, and so on) can reveal the effect of different policies and population structures on the evolution of the epidemic.In our framework, users do not have access to the items in the shaded area; they know only the shapes of the optimizer, model, and system and are unable to act on that information.
We discuss in detail the methodology we used as well as related work in section 2. In section 3, we describe our experiments along with the application of our proposed model to real-world datasets.In section 4, we discuss our results and the importance of our study.We conclude in section 5 with a summary of the key points.

Methods
The most commonly used data-driven control technique is based on MPC (Figure 1).It enables the controller or user to interact with a dynamical system and affect its evolution by applying a specific control variable at selected times or epochs depending on the predicted horizon.Our approach was inspired by this technique but differs in two key aspects: we do not have direct access to the system, so we cannot manipulate its evolution, and we do not have any insight into the controller's strategy or parameters.(Here we focus on a Kalman-type dynamical system, but our approach can be applied to any parametric/non-parametric dynamical system).Our knowledge is limited to the fact that the controller or the system uses a model from a specific system family (here, for example, one in the form of equation ( 2)).That is, we have access to only the observed final output data.This idea is sketched in Figure 2.
This section is organized as follows; We begin in by reviewing related work and classical optimal control techniques in Sections 2.1 and 2.2.In Sections 2.3 and 2.4 we describe our proposed hidden control framework.

Targeted control problem and its formalism
We assume a discrete-time context in which decision-makers (e.g., government agencies) observe the distribution of the disease variable (e.g., the numbers of weekly infected, recovered and deceased individuals), possibly with observational inaccuracy, for a given time period (day, week, month, etc.).For each sub-population P i , the observed distribution is generally estimated empirically from data that we denote The estimated evolutionary distributions for the whole population can be represented as π ℓ i (k), we are interested in the dynamical evolution of the discrete n-dimensional probability vector Because we estimate the optimal control by sub-population, we omit the suffix i hereafter to simplify the description (e.g., x i (k) is represented as x k , and a column of M k as M k ).This simplified description represents the state vector of the disease distribution at time k.To model the evolution, we propose using a discrete Kalman-like controlled system: • The term u k ∈ R m is an m-dimensional vector describing the unknown optimal control.This unknown control is the hidden mechanism that led to the observed shape of the evolution.
• The deterministic state matrix A = (a ij ) i,j=1,...,n , defined in R n×n , and describing the logical evolution of x k when no control is applied is also unknown.
• The unknown control matrix B = (b ij ) i=1,...,n,j=1,...,m , defined in R n×m , describes how the control is applied to the system to bring it to its targeted state.
The idea is to use only observed data M k in (1) to simultaneously estimate all the unknown objects of this control problem, namely matrices A, B and the control term u k that leads to the observed shape of the data.When A and B are known, the optimal control û * k can be obtained using classical techniques that have been widely discussed in the literature (e.g., Brogan (1991); Bryson and Ho (2018);Chen (1984).In our case, these parameters are unknown and should be estimated from the observed data.Our strategy is to use the optimal control results obtained using an iterative algorithm that minimizes a cost function.

Control algorithm and related work
For given deterministic matrices A and B, the solution of the discrete-time optimal control problem (2) of a linear system is described in Brogan (1991);Chen (1984).It consists of using a linear optimization scheme (2) and finding the optimal control u = (u 1 , . . ., u K ) that minimizes the objective or cost function: where exponent (.) ′ represents the transpose, and Q, R, and S are cost matrices defined as follows.
• Q ∈ R n×n is a positive definite weighting matrix associated with the state vector.It is linked to the relative importance assigned to state vector x k in cost function J(•).In other words, a higher value in the diagonal of this matrix indicates a bigger penalty or greater emphasis placed on minimizing related state vector x k .• R ∈ R m×m is a positive definite weighting matrix affecting the relative relevance or weight assigned to control u k in cost function J(•).It specifies how much emphasis the control algorithm should place on minimizing the control effort.A higher value in its diagonal, similar to that in the Q matrix, indicates a higher penalty or greater importance placed on reducing the control input.• S is a positive definite weighing matrix linked to the final state vector.It is a cross-weighting matrix used by cost function J(•) to penalize discrepancies between the inputs from state x k and control u k .It is aimed at identifying any cross-dependencies between the states and control inputs.We are not interested in the final state here, so we suppose that S = 0.
The optimal control problem has been widely studied (e.g., Brogan (1991); Bryson and Ho (2018)).In our case, a linear quadratic regulator variant is used to minimize (3).The solution of the state space equation (2) of a discrete time-invariant linear system has been given in a more general framework Fadali and Visioli (2012) for 2 ≤ k ≤ K: This equation expresses state vector x k in terms of the initial state vector x 1 and a sequence of control vectors, u i , i = 1, . . ., k, that are implicitly dependent: A, B and weighting matrices Q, R. In this setting, the controlled system, (2), with dimension K is controllable if the controllability matrix, ( 0 ,  0 ) ( 0 ,  0 ) For , , minimize eq.(3) as LQR solution Calculate  ,  (eq.( 5)) Calculate   * (eq.( 7)) Estimation: Find ,  to minimize ( * , ) under constraints for probabilities of  * and  and so on (eq.( 9)) Algorithm 1 Algorithm 2 Figure 3: Structure of main algorithm for estimating hidden control process has column rank K (i.e., K linearly independent columns).Under these conditions on parameters A, B, the optimal control solution Fadali and Visioli ( 2012) is given by In this setting, the closed-form expression of the gain matrix K is obtained using a Linear Quadratic Regulator (LQR).This is done iteratively in two steps using a backward approach, as detailed in Algorithm 1: • Gain matrix K is calculated from evolutionary parameter matrices A, B, Q, and R via S k+1 at step k: • S k is obtained from the solution of the Riccati equation Bittanti et al. (2012): By combining equations ( 2) and ( 6), we can write the evolutionary structure of the controlled system (2) recursively: In the general framework of optimal control, the parameters are known and fixed.In our study, we are dealing with unknown parameters.Therefore, our objective becomes the estimation of these unidentified matrices in the Kalman filter parameters, namely A, B, as well as weighting matrices Q, R.

Our hidden controller framework
Instead of the traditional objects found in control problems, our methodology relies on using discrete probability distributions as the dynamic system object, as shown in (2), to characterize the system's behavior.A further distinguishing characteristic of our methodology is its ability to estimate both the parameters of the model and their corresponding optimal control, relying exclusively on observed data.To determine unknown parameters A and B, our approach utilizes two interrelated optimization steps: 1.The first step is to find optimal control features K(A, B)) and x * k , as shown in ( 9) under controllability conditions on A and B. 2. The second step is to minimize a suitable distance that has been adapted to probability distributions between x * k and observation x k under two constraints: x k is a discrete probability distribution calculated from observed data M k and the system is controllable.This leads to the following minimization problem: The constrained minimum is taken over all A and B such that the system is controllable, as detailed in Algorithm 2.
In accordance with the concept of the probability distribution, distance d(., .) in ( 10) determines the degree of similarity between two probability distributions.A variety of metrics and distances that have been used in different circumstances and scenarios have been reported.We next discuss five of them.

Metrics and distances
• The total variation distance measures the maximum difference between the probabilities that two distributions assign to an event.Therefore, it gives a worst-case measure of the difference between the two distributions, which makes it a conservative discrepancy measure.It is given by For the discrete case, the right-hand side of ( 11) can be expressed as the ℓ 1 distance.Despite its frequent usage, the TV distance metric has certain disadvantages.It is sensitive to absolute differences between probabilities, which can lead to misleading results when applied to large event sets.Though it satisfies several metric properties like non-negativity, symmetry, and the triangle inequality, it is not a "proper" metric as it can assign a distance of zero to non-identical distributions.The TV distance also lacks the ability to capture dependencies among variables in multivariate distributions and does not provide information about the 'location' and 'spread' of differences.Furthermore, it is not a differentiable function, posing challenges for optimization algorithms.Lastly, its computational complexity increases substantially for continuous distributions, especially in high-dimensional scenarios.
• Kullback-Leibler (KL) divergence quantifies the information loss when distribution D 1 is used to approximate D 2 : .
KL divergence has several advantages and disadvantages when used as a measure of the difference between two probability distributions.Among its advantages are its high sensitivity to discrepancies in probability distributions, including those occurring in distribution tails, and its widespread application across multiple fields, such as information theory and machine learning.It also has a deep connection with maximum likelihood estimation, which enhances its value for theoretical analysis.The disadvantages include its lack of symmetry; that is, the divergence from distribution D 1 to D 2 is not the same as that from D 2 to D 1 , which can make it counterintuitive as a "distance" measure.Furthermore, it is not a metric as it does not satisfy properties like symmetry and the triangle inequality.Another critical disadvantage is its undefined nature for non-overlapping distributions.Finally, when probability distributions are estimated from data, KL divergence can be excessively sensitive to the sample size, leading to potentially unreliable estimates for smaller sample sizes.
• The Wasserstein distance, also known as the earth mover's distance, measures the minimum amount of work needed to transform one pile into another pile.For continuous or discrete probability distributions D 1 and D 2 , the p-Wasserstein distance (p ≥ 1) is defined as .
The inf bound is taken over all random variable couples (X, Y ) such that the distribution of X is D 1 and the distribution of Y is D 2 .It offers an intuitive geometric interpretation, acting as a measure of the minimum cost required to transform one distribution into another.In addition, it captures shifts in location, a property not captured by other metrics, making it advantageous in cases in which location is important.Moreover, unlike measures such as KL divergence, the Wasserstein distance remains defined even for non-overlapping distributions and exhibits robustness to noise.However, its computation can be resource-intensive due to the underlying optimization problem.It thus poses challenges, particularly in higher dimensions.The distance also depends on the choice of a cost function, which may not always have an evident best choice and can affect the distance calculation.Although sensitivity to location shifts can be beneficial, in instances where location is irrelevant, this sensitivity could capture unnecessary differences.Lastly, like many other measures, the Wasserstein metric faces scalability issues with high-dimensional data, demanding dimension reduction techniques to be effective.
• TheHellinger distance measures the total difference between the square roots of the probabilities in the two distributions: . The Hellinger distance has several advantages.It possesses metric properties including symmetry and the satisfaction of the triangle inequality, making it a true metric.It is less sensitive to outliers due to the comparison of square roots of probabilities instead of the probabilities themselves.Its fixed range, typically between 0 and 1, aids in comparing distances across various distribution pairs.Moreover, its simplicity when comparing Gaussian distributions offers computational efficiency.However, the Hellinger distance is not without drawbacks.Its structure, rooted in square roots of probabilities, makes it less sensitive to differences in distribution tail behavior.It measures global dissimilarity, often failing to capture local differences effectively.Lastly, unlike KL divergence, it does not quantify information loss when one distribution approximates another, a property that may be desirable in certain contexts.
• Theχ 2 distance is the sum of the squared differences between the two distributions, divided by the first distribution.It is not symmetric.
Each of these metrics comes with its own strengths and limitations.For instance, KL divergence (d KL ) and the chi-squared distance (d χ 2 ) are non-symmetric, making them less intuitive for certain applications.On the other hand, metrics like the TV distance (d T V ) tend to be more conservative.In this study, our primary objective was to ensure robust convergence and computational efficiency in the optimization process.To identify the most suitable metric, we conducted comprehensive experiments testing these distances on our dataset.We encountered challenges regarding convergence and stability during the optimal control and estimation processes, as shown in (10).Our experiments revealed that the Hellinger distance (d H ) exhibits superior performance in terms of convergence, stability, and efficiency.This can be attributed to its connection with the Euclidean ℓ 2 norm, which provides certain computational advantages.Therefore, it is the primary distance measure we used in this study.
arXiv Template A PREPRINT comparison of square roots of probabilities instead of the probabilities themselves.Its fixed range, typically between 0 and 1, aids in comparing distances across various distribution pairs.Moreover, its simplicity when comparing Gaussian distributions offers computational efficiency.However, the Hellinger distance is not without drawbacks.Its structure, rooted in square roots of probabilities, makes it less sensitive to differences in distribution tail behavior.It measures global dissimilarity, often failing to capture local differences effectively.Lastly, unlike KL divergence, it does not quantify information loss when one distribution approximates another, a property that may be desirable in certain contexts.
• The 2 distance is the sum of the squared differences between the two distributions, divided by the first distribution.It is not symmetric.
Each of these metrics comes with its own strengths and limitations.For instance, KL divergence (d KL ) and the chi-squared distance (d 2 ) are non-symmetric, making them less intuitive for certain applications.On the other hand, metrics like the TV distance (d T V ) tend to be more conservative.In this study, our primary objective was to ensure robust convergence and computational efficiency in the optimization process.To identify the most suitable metric, we conducted comprehensive experiments testing these distances on our dataset.We encountered challenges regarding convergence and stability during the optimal control and estimation processes, as shown in (10).Our experiments revealed that the Hellinger distance (d H ) exhibits superior performance in terms of convergence, stability, and efficiency.This can be attributed to its connection with the Euclidean `2 norm, which provides certain computational advantages.Therefore, it is the primary distance measure we used in this study.
Algorithm 1 Computing optimal control and its corresponding trajectory (see Figure 3) Require: Estimated distributions M k in (1) from observed data.
Require: A vector Y of dimension n(n + m + 2) of the system parameters.
1: Retrieve A, B, Q, R from Y .
Ensure: Q, R are symmetric positive definite.
12: Output: x ⇤ = (x ⇤ 1 , . . ., x ⇤ K ), the controlled trajectory.=0 Our parameter estimation: To prepare for the calculation of matrix parameters A, B of the hidden closed-loop equation ( 9) on the basis of observed data, we vectorize them: (a 11 , . . ., a nn , b 11 , . . ., b nm ) 0 2 R n 2 +nm .Furthermore, we assume that the hidden controller uses weighting matrices Q and R to calculate the gain matrix K in (7).There is no systematic way to

Our parameter estimation:
To prepare for the calculation of matrix parameters A, B of the hidden closed-loop equation ( 9) on the basis of observed data, we vectorize them: (a 11 , . . ., a nn , b 11 , . . ., b nm ) ′ ∈ R n 2 +nm .Furthermore, we assume that the hidden controller uses weighting matrices Q and R to calculate the gain matrix K in (7).There is no systematic way to make this choice, as it depends only on the user's preferences in terms of control over parameters A and B. Most often, the weighting matrices Q and R are chosen to be diagonal.Due to the small number of observations, we use Q = diag(α 1 , . . ., α n ) and R = diag(β 1 , . . ., β n ).We then rearrange all the free parameters in a vector: Y = (a 11 , . . ., a nn , b 11 , . . ., b nm , α 1 , . . ., α n , β 1 , . . ., β n ) ′ ∈ R n(n+m+2) .Let us now consider the objective cost function: where d(., .) is one of the distances d T V , d H , d χ 2 , d Wp , or d KL , and ′ ∈ R n is the statistically estimated distribution for population segment i at observation time k.(Index i is omitted for M k and x k to simplify the description, as mentioned in section 2.1.) • x * k is the closed-loop state of the controlled system.The computation of vector x * k for k = 1, . . ., K is shown in (9).It is implicitly dependent on the parameters in vector Y .This step is nested inside the general optimization problem (10).Vector Y should satisfy two constraints, C 1 and C 2 , on the one hand, to comply with both the controllability constraints on A and B and the positive definiteness of Q and R, and, on the other hand, to satisfy the probability structure of state vector x k .
C 1 : Vector Y should verify the constraint that A and B lead to a controllable system by verifying (5).
C 2 : For all k, vector x * k calculated using (9) should satisfy a probability structure, meaning that all its elements should be positive and their sum is less than 1.
The final optimal parameter vector, Y , is then obtained by solving the nonlinear optimization problem under constraints: From the practical point of view and to reduce the relative error in estimating the unknown coefficients of the parameters, we assume that the hidden controller resets vector x k to the observed value each time when calculating it for interval K.
To implement this assumption, we divide the range of time-ordered observations x k by a pre-selected value T (e.g., for COVID dataset, T = 3 (three days) or T = 7 (one week)) and denote the number of observation groups formed by L = [ K T ], where [x] is the integer part.We then consider the partition group {G 1 , . . ., G L } on all the observations {M 1 , . . ., M K } ordered in accordance with time: the observations in group G 1 being G 1 = {M 1 , . . ., M T }, those in group G 2 being G 2 = {M T +1 , . . ., M 2 * T }, and those in the last group being G nstep = M (L−1) * T +1 , . . ., M L * T , for example.For any 1 ≤ ℓ ≤ L, we take vector x k as follows: • x k = M k for k = (ℓ − 1) * T + 1.Thus, we reset vector x k to the first observed value in each group.
For the solution of the nonlinear optimization problem (10), we use sequential quadratic programming, which is an efficient and widely used technique de Gouvêa and Odloak (1997); Wright et al. (1999); Fletcher (2013); Schittkowski and Zillober (2003), Moreover, it is readily available in the MATLAB software used for the optimization part of this work.Call Algorithm 1 with Y , and get K(A, B) and x ⇤ = (x ⇤ 1 , . . ., x ⇤ K ).

7:
Solve quadratic programming subproblem under constraints C 1 and C 2 . 8: Update Y and calculate its corresponding F (Y ).
11: Output: Corresponding optimal control c u ⇤ .=0 These algorithms were applied to each segment of the population.The trajectories estimated using these algorithms and the observed data for different countries are plotted in Figure 5.

Interpretation of control parameters
In accordance with the preceding discussions, we computed the controlled part, denoted here as B, c u ⇤ , of the resulting dynamical system.This computation enabled us to estimate the hidden parameters, and in turn, uncovered the controller's strategy.Another key accomplishment of this study was correlating these controlling elements with the characteristics of the corresponding population.The next step was to, a posteriori, explore how these elements affected the evolution of the dynamical system.To this end, we fitted the derived controlled part to the observed population characteristics as illustrated in Figure 3.For each population segment P i , this was done using a multivariate GAM as such GAMs are well-known for their adaptability in dealing with various forms of data and their capacity to model nonlinear relationships.We conducted extensive statistical testing on the GAM output to identify the set of variables closely linked with the control part.
Multivariate GAM Let us consider a multivariate GAM in which the response is vector X = B c u ⇤ 2 R n and denote C (1) , . . ., C (J) as the dependent variables describing the characteristics of the population.As shown in Table 1, we considered six characteristic COVID-19 variables specific to Japanese prefectures.The responses and covariates (predictors) are observed at times k = 1 . . .K.
The functions S 1 (.), . . ., S J (.) are flexible smooth functions and can fit complex patterns in the data that would not be captured by a linear model.At the same time, they maintain interpretability: the fitted functions can be plotted and easily interpreted.The specifics of fitting a GAM depend on the type of response variable (whether it is binary, count data, continuous, etc.), the chosen link function, and the chosen form for the nonlinear functions S 1 (.), . . ., S J (.).This model has been well studied, and many software packages are available for fitting GAMs, such as Mixed GAM Computation Vehicle with Automatic Smoothness Estimation (mgcv) in R ??.As depicted in Figure 4, the application of the GAM yields an array of newly minted comprehensible features designated as z (1) , . . ., z (J) .These features serve to quantify the intensity of the association between the control and corresponding population characteristics, thereby offering a method to gauge the extent of a particular population characteristic's effect on the control of the obscured dynamical system.9 These algorithms were applied to each segment of the population.The trajectories estimated using these algorithms and the observed data for different countries are plotted in Figure 5.

Interpretation of control parameters
In accordance with the preceding discussions, we computed the controlled part, denoted here as B, u * , of the resulting dynamical system.This computation enabled us to estimate the hidden parameters, and in turn, uncovered the controller's strategy.Another key accomplishment of this study was correlating these controlling elements with the characteristics of the corresponding population.The next step was to, a posteriori, explore how these elements affected the evolution of the dynamical system.To this end, we fitted the derived controlled part to the observed population characteristics as illustrated in Figure 3.For each population segment P i , this was done using a multivariate GAM as such GAMs are well-known for their adaptability in dealing with various forms of data and their capacity to model nonlinear relationships.We conducted extensive statistical testing on the GAM output to identify the set of variables closely linked with the control part.
Multivariate GAM Let us consider a multivariate GAM in which the response is vector X = B u * ∈ R n and denote C (1) , . . ., C (J) as the dependent variables describing the characteristics of the population.As shown in Table 1, we considered six characteristic COVID-19 variables specific to Japanese prefectures.The responses and covariates (predictors) are observed at times k = 1 . . .K.
The functions S 1 (.), . . ., S J (.) are flexible smooth functions and can fit complex patterns in the data that would not be captured by a linear model.At the same time, they maintain interpretability: the fitted functions can be plotted and easily interpreted.The specifics of fitting a GAM depend on the type of response variable (whether it is binary, count data, continuous, etc.), the chosen link function, and the chosen form for the nonlinear functions S 1 (.), . . ., S J (.).This model has been well studied, and many software packages are available for fitting GAMs, such as Mixed GAM Computation Vehicle with Automatic Smoothness Estimation (mgcv) in R Wood (2011); Wood et al. (2016).As depicted in Figure 4, the application of the GAM yields an array of newly minted comprehensible features designated as z (1) , . . ., z (J) .These features serve to quantify the intensity of the association between the control and corresponding population characteristics, thereby offering a method to gauge the extent of a particular population characteristic's effect on the control of the obscured dynamical system.
These new features provide valuable information for investigating the dynamical system's behavior across various population segments.This facilitates the identification of similarities amongst segments for which the system exhibits For each population segment  1 ∶ ,  *  2 ∶ ,  * ⋮   ∶ ,  *

Heatmaps Dendrograms Classification
Figure 4: For each population segment, the link between control terms B, u * , and C (1) , . . ., C (J) , which describe the characteristic variables (covariates, predictors) of the population segment, was fitted to a GAM.The z-score output of the GAM was used for classification to identify how these predictors contributed to the control of the system.3 Experiments and application to real data

Data description
We used two sets of observed COVID-19 data: (i) data from five prefectures in Japan provided by JX PRESS Corporation and (ii) data from nine countries sourced from a "worldwide epidemiological database for COVID-19" Guidotti (2022) and the "COVID-19 data hub" Guidotti and Ardia (2020).Using our hidden controller framework, we estimated the unknown parameters A, B and the gain matrix K of the closed-loop system (9) to capture the dynamic behavior of COVID-19.
The prefectural data used in this study comprised the proportion of weekly infected individuals in prefecture i (π 1 i (k)), the proportion of weekly recovered individuals (π 2 i (k)), the proportion of weekly deceased individuals (π 3 i (k)), and the proportion of all other people who never had contact with the virus (π The proportions were estimated from weekly counts of numbers reported daily for 32 weeks from April 1, 2021, to November 10, 2021, coinciding with the vaccination campaign.The five prefectures (i = 1, . . ., 5) comprise Tokyo, Osaka, Hokkaido, Fukuoka, and Okinawa, respectively.These prefectures are representative of Japanese populations for which there were a sufficient number of cases for statistical analysis.For each prefecture, the data are represented as a 32 × 3 matrix, with the weekly observations in the columns.The candidate weekly predictor covariates listed in Table 1 are for the period two weeks prior.Factors such as increases in driving and public transport ("transit") in Tokyo, which are considered to affect the rate of new infections, are included.We used a two-week delay between measurement and onset in accordance with previous findings Matsui et al. (2022).
We selectively used data from nine countries in our experiments: Australia (AUS), Brazil (BRA), Chile (CHL), Colombia (COL), Czech Republic (CZE), Germany (DEU), Lithuania (LTU), South Africa (ZAF), and Japan (JPN).The proportion of weekly infected individuals in country i (π 1 i (k)), the proportion of weekly recovered individuals , and the proportion of weekly deceased individuals (π 3 i (k)) were calculated from the data on "confirmed" (cumulative number of confirmed cases), "deaths" (cumulative number of deaths), and "recovered" (cumulative number of patients released from hospitals or reported recovered) Guidotti (2022); Guidotti and Ardia (2020).The period was the same as for dataset (i): the 32 weeks from April 1, 2021, to November 10, 2021.The candidate predictor covariates for these nine countries are shown in Table 3.Four of them represent governmental policy measures: government_response_index, stringency_index, containment_health_index, and economic_support_index.Since these policy measures are correlated and the ranges of the index scores vary between countries, we first categorized the policy measures into three levels -low (L: 0-50 percentile), middle (M: 50-90 percentile), and high (H: 90-100 percentile), using the Japanese data as a reference.We then applied multiple correspondence analysis and hierarchical clustering on principal components to the categorized data and obtained three groups of policy measures, as shown in Table 2.We again used a two-week delay between measurement and onset in accordance with previous findings Matsui et al. (2022).

Experimental conditions
To calculate parameters A, B of the closed-loop equation ( 9) on the basis of the observed COVID-19 data for the five Japanese prefectures (i) and nine countries (ii), we first defined Y = (a 11 , a 12 , . . ., a 33 , b 11 , b 12 , . . ., b 33 ) T ∈ R 18=3 * 3+3 * 3 representing the unknown coefficients of parameters A, B. We then set the Hellinger distance (described in section 2.3) to d(x * k , x k ), where the term x * k represents the estimated COVID-19 dynamics, and the other term ′ ∈ R 3 represents the observed COVID-19 dynamics on the k th day.
The proportions Π(k) were calculated using the total population of Japan and the population ratios of each prefecture as of January 2021.We used the fmincon programming solver in MATLAB and selected the interior-point algorithm, for which the parameters of the termination tolerances on first-order optimality, function value, and input matrices were set to 1e − 4. The function optimized in fmincon was estimated for each window with a length of 3 using the ℓ 2 -norm.We used mgcv in R Wood (2011); Wood et al. (2016) to interpret the control parameters (section 2.4).

Results and discussion
4.1 Evaluating predictive power of hidden controller framework Figure 5 shows comparisons between our hidden controller framework's prediction (x * k ) and the observed x k for Australia, Colombia, Germany, and Japan based on the data for nine countries (data set (ii)).The predictions are reasonably accurate, suggesting that model parameters A, B were well-estimated.Similar results were observed for the other five countries, as detailed in Appendix A. These results underscore the predictive capacity of our hidden controller framework.Furthermore, control feedback u k can be effectively used in subsequent data analyses.

Predictor covariates aimed at halving the number of infected and deceased individuals
To analyze which predictor covariates effectively reduce the number of infected and deceased individuals, we applied our hidden controller framework (section 2.3) and a multivariate GAM (section 2.4) to a dataset in which the number of infected and deceased individuals was halved (dataset (i)).Figure 6 shows heatmaps depicting the z-scores of the predictor covariates associated with the numbers of (a) infected and (b) deceased individuals.Since z-scores approximately follow a Student's t-distribution and transition to a normal distribution when the degree of freedom becomes large, we calculated the z-scores at a significance level of 0.05, assuming a normal distribution.We color-coded the heatmaps blue for smaller z-scores, representing positive effects in decreasing the number of infected and deceased individuals, and red for larger z-scores, representing negative effects that increase these numbers.
We roughly categorized the predictor covariates into two groups, F1 and F2.F1 includes "week," "transit," "walking," and "num_beds," which generally had negative effects.In contrast, the predictor covariates in F2 ("driving" and "vaccination") had positive effects.These results suggest that, in Japan, traveling by automobile and getting vaccinated effectively reduced the number of infected and deceased individuals.
We can classify the prefectures into two groups: the major metropolises of Tokyo, Osaka, and Fukuoka, and the tourist cities in Okinawa and Hokkaido.The data suggest that vaccination in metropolitan areas is highly effective.
To further investigate the positive and negative effects on reducing the number of infected and deceased individuals, we analyzed the difference in z-scores for predictor covariates associated with both categories, with the goal of halving each compared with the observed numbers.Figure 7 presents heatmaps illustrating these differences in z-scores.For the infected category, vaccination was highly effective in the metropolises of Tokyo, Osaka, and Fukuoka.For the deceased category, vaccination was particularly effective in Tokyo.

Differences in predictor covariates across countries
Figure 8 presents a heatmap of z-scores for predictor covariates associated with the number of infected, deceased, and recovered individuals across the nine countries.As an illustration of the classification of characteristics based on control (section 2.4), the countries can be broadly categorized into four clusters: F1 includes South Africa, Chile, Brazil, and Colombia; F2 is comprised solely of Japan; F3 encompasses Australia and Lithuania; and F4 contains the Czech Republic and Germany.For group F4, three policy measures in Table 2 (P1 as indicated in the figure) and vaccination (P2) seem to effectively control the number of infected, deceased, and recovered individuals.However, for countries in group F1, the opposite effects were observed.

Clustering of countries on basis of policy measures
Figure 9 illustrates the clustering of countries into four groups on the basis of the three categories of policy measures outlined in Table 2.These groups align with the ones determined by the z-scores in the previous section.This underscores the critical role of policy measures in controlling the number of infected, deceased, and recovered individuals.

Conclusion
We have presented a methodology that enables the correlation of various population characteristics with the evolution of a dynamically controlled system.We assume that the observed dynamical system is the result of optimal control of a Kalman-type model, with no direct access to the various stages of this control.The choice of a Kalman-type model was made for simplicity, but this concept can be readily generalized to other models.We assume knowledge of only the model's structure, not its parameters.Our approach, therefore, is to estimate what the controller has done to arrive at the observed data; the estimates of the control parameters are made a posteriori.
A unique aspect of this work is that the evolutionary objects of the dynamical system are probability distributions.We explore how these distributions have evolved under different hidden controller settings.We have provided an algorithm   for estimating the various parameters that control the evolution of these distributions, which are based solely on a posteriori observations.
Once the control parameters were obtained, we used a multivariate generalized additive model (GAM) to explore potential relationships between the part that controls the system's evolution and the various characteristics of the studied population.We applied this methodology to a COVID-19 database for five prefectures in Japan, considering various explanatory variables.This approach was generalized to nine countries in order to explore the effects of the various policies they follow.
From the GAM, we extracted indicators to quantify the degrees of association between the different explanatory variables.These new indicators enable the analysis and visualization of groups with similar and different evolutionary profiles.The aim is to explain the reasons for these differences.This work provides a novel perspective on understanding and managing complex evolutionary processes.

Figure 1 :
Figure 1: In the MPC framework, users can access and manipulate (i.e., change or control) the system parameters.
Estimating hidden control parameters and trajectoryRequire: Estimated distributions M k in (1) from observed data.1:Define a convergence threshold ".2: Initialize A, B, Q, R.

Figure 5 :
Figure 5: Comparison between predictions made by our hidden controller framework (green dashed lines) and observed data (red solid lines) for (a) Australia, (b) Colombia, (c) Germany, and (d) Japan.

Figure 6 :
Figure 6: Heatmaps depicting z-scores of predictor covariates associated with numbers of (a) infected and (b) deceased individuals, with a goal of reducing the counts in both categories by half.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Heatmaps illustrating difference in z-scores for predictor covariates associated with numbers of (a) infected and (b) deceased individuals, aimed at reducing each category by half compared with the observed numbers.

Figure 10 :
Figure 10: Comparison between the predictions made by our hidden controller framework (green dashed lines) and observed data (red solid lines) for remaining five countries.

Table 1 :
Candidate weekly predictor covariates for five Japanese prefectures.
num_beds Number of available beds in hospitals in Tokyo (provided by Ministry of Health, Labour and Welfare ) vaccin Number of people vaccinated in Tokyo (provided by Ministry of Health, Labour and Welfare ) identical behavior.Various classification techniques have been explored for achieving this functionality.More details, illustrative figures, and discussion can be found below on the experiments performed on actual COVID-19 data.

Table 2 :
Main categorical levels consisting of three groups of policy measures.