Minimum Information Variability in Linear Langevin Systems via Model Predictive Control

Controlling the time evolution of a probability distribution that describes the dynamics of a given complex system is a challenging problem. Achieving success in this endeavour will benefit multiple practical scenarios, e.g., controlling mesoscopic systems. Here, we propose a control approach blending the model predictive control technique with insights from information geometry theory. Focusing on linear Langevin systems, we use model predictive control online optimisation capabilities to determine the system inputs that minimise deviations from the geodesic of the information length over time, ensuring dynamics with minimum “geometric information variability”. We validate our methodology through numerical experimentation on the Ornstein–Uhlenbeck process and Kramers equation, demonstrating its feasibility. Furthermore, in the context of the Ornstein–Uhlenbeck process, we analyse the impact on the entropy production and entropy rate, providing a physical understanding of the effects of minimum information variability control.


Introduction
Time-varying probability density functions (PDFs) are a preferred approach for characterising the dynamics of various complex systems.PDFs commonly feature in emergent fields, such as active inference [1] or stochastic thermodynamics [2], where their value is derived either through data analysis or by solving the Fokker-Planck (FP) equation associated with an Itô/Stratonovich stochastic differential equation.
Drawing upon control theory [3], when the dynamics of the system under consideration are governed by an FP equation, we can explore control objectives such as the regulation (setting to a constant value) or tracking (following a time-varying reference) of the systems' time-varying PDFs [4].While the notion of controlling PDFs may initially appear impractical using conventional control engineering methods, advancements in technologies like optical tweezers have rendered it feasible, which is particularly evident in applications such as colloidal systems [5][6][7][8], with specific implications in biomolecule evolution control [9].
Regarding the control of FP equations, seminal works [10,11] present methodologies to control the system PDF governed by FP equations [12].Building upon this foundation, ref. [13] discusses a bilinear optimal control problem where the control function depends on time and space.In [14], the authors prove the existence of optimal controls by considering first-order necessary conditions in the optimisation problem.
In applications like Brownian motion, we can control FP equations via reverseengineering approaches such as the engineered swing equilibration (ESE) method [15,16].The ESE protocol imposes a solution to the FP equation to obtain the corresponding PDF's control parameters that provide a shortcut for time-consuming relaxations (for further details, refer to the methods section of [17]).Ref. [8] offers a comprehensive review of similar reverse-engineered approaches.
Since FP equations frequently serve as mathematical descriptions of mesoscopic systems (for further details, see [12]), i.e., systems at the nano/microscale such as molecular motors, the time evolution of the system's PDF may not only need to adhere to time constraints but also to multiple "thermodynamic" constraints to be deemed "efficient".For instance, the system may require the minimisation of the entropy production [18,19], a reduction in information variability [20], or an enhancement of self-organisation [21].The incorporation of these "thermodynamic constraints" into the optimisation process implies an extension of the current literature results on the control of FP equations.
A theory that could offer insights into addressing these optimisation challenges stems from information geometry, a field resulting from the fusion of information theory and differential geometry [22].As an evolving field, information geometry proposes novel solutions to tasks, such as maximum likelihood estimation [23], state prediction [24,25], the quantification of causality [26][27][28], or maximum work extraction [4,29].In stochastic thermodynamics [2,30], information geometry aids in obtaining time-varying descriptions of the aforementioned constraints.For instance, based on the well-known Cauchy-Schwartz inequality [31], ref. [32] presents an inequality between the Fisher divergence [33] and the information length (IL) [25,34] to quantify the disorder in irreversible decay processes.Ref. [35] introduces an inequality describing the information rate as a speed limit on the evolution of any observable.In [4], the geodesic of the IL describes the path with the least statistical variation connecting initial and final probability distributions of the system dynamics (for further details, see [20]).Consequently, information geometry can be employed in a control protocol to enforce geodesic dynamics on the system's PDF time evolution, achieving the minimum "geometric information variability" [21] and thereby establishing an optimal speed limit for the system's observables.The primary challenge addressed in this study revolves around devising and applying a technique to achieve this minimum geometric information variability.
Developing an optimal protocol for the time evolution of the system's PDFs under multiple constraints requires an exploration of existing control algorithms.The literature presents a significant amount of control procedures, spanning from classical PID control [36] to more sophisticated algorithms, like data-driven, model-free, or fractional-order controls (for instance, see [37][38][39]).However, we require an algorithm capable of handling intricate optimisation problems while remaining a feasible option for implementation in future experimental setups.One of the most popular optimisation-based control techniques is the model predictive control (MPC) scheme [40].Generally, MPC is an online optimisation algorithm designed for constrained control problems, whose advantages have been observed in applications within robotics [41], solar energy [42], or bioengineering [43].Furthermore, MPC can be easily implemented thanks to packages such as CasADi [44] or the Hybrid Toolbox [45].
Based on the presented discussion, this work offers a solution to an optimisation problem that integrates the concepts of the information length and the quadratic regulator (QR) [46] to guide the system's PDF time evolution along the path with the minimum geometric information variability (the geodesic of the information length) using MPC.Although we could implement the IL, QR, or MPC in more complex scenarios, we only study their application to stochastic processes described by linear stochastic differential equations.Hence, the system's PDF will consistently maintain a Gaussian distribution over time, assuming that the system's initial conditions follow a Gaussian distribution.Constraining the analysis to linear stochastic dynamics allows us to use a set of deterministic differential equations to describe the evolution of the mean and covariance of the Gaussian distribution within the MPC method.The applicability of such a prediction model favours MPC over other alternatives, such as reinforcement learning [47], which could have been explored to determine the geodesic of the information length within a similar context.Additionally, MPC provides a low offline complexity, mature stability-feasibility-robustness theory, and good constraint handling [47].
The algorithm is applied to the Ornstein-Uhlenbeck process [48] and the Kramers equation [25], which both describe a particle over a heath reservoir (mesoscopic stochastic dynamics).In practice, changes in temperature and optical tweezers can manipulate the dynamics of both the noise amplitude and mean value in such systems, respectively [2,7,49].Using previous findings from [50], the effects of the MPC method on the Ornstein-Uhlenbeck process are analysed by comparing IL values with the entropy production and entropy rate in the closed-loop (feedback-controlled) system.As a rigorous closed-loop stability analysis is beyond the scope of this work, we provide only a brief description of the BIBO stability conditions that are considered to constrain the control actions proposed by the MPC.Finally, we give a set of concluding remarks and discuss future work.
To help readers, in the following, we summarise our notations.R is the set of real numbers.x ∈ R n represents a column vector x of real numbers of dimension n, A ∈ R n×n represents a real matrix of dimension n × n (bold-face letters are used to represent vectors and matrices), Tr(A) corresponds to the trace of the matrix A, and A ij is the element at row i and column j of the matrix A. |A|, vec(A), A ⊤ , and A −1 are the determinant, vectorisation, transpose, and inverse of matrix A, respectively.The value I n denotes the identity matrix of order n.Newton's notation is used for the partial derivative concerning the variable t (i.e., ∂y ∂t = ẏ).In addition, the average of a random vector ζ is denoted by µ := ⟨ζ⟩, with the angular brackets representing the average.Finally, in Table 1, we provide a brief description of the variables used throughout this paper.

Preliminaries
As mentioned in Section 1, this work considers systems whose dynamics are governed by linear non-autonomous stochastic differential equations described in the form of the following set of Langevin equations: where A and B are n × n and n × 1 real time-invariant matrices, respectively; u(t) is a (bounded smooth) external input, and ξ ∈ R n is a Gaussian stochastic noise given by an n dimensional vector of δ-correlated Gaussian noises ξ i (i = 1, 2, . . .n), with the following statistical property: Readers seeking a comprehensive study of systems governed by Equation ( 1) can refer to [51].Additionally, for a concise introduction to linear time-varying and time-invariant systems, ref. [52] provides valuable insights.The Fokker-Planck equation of (1) can also be utilised to depict the system's PDF dynamics, characterised by both time t and the spatial variable vector x = [x 1 , x 2 , . . ., x n ] ⊤ .This is given as follows [53,54]: If the initial PDF is Gaussian, a solution to (3) for p(x; t) at any time t is given by [53] p Equation ( 4) describes the dynamics of the probability distribution of the random variable ζ governed by (1).The value of the mean µ(t) and covariance Σ(t) in a linear stochastic process can be computed from the solution of the following set of differential equations [51]: where D ∈ R n×n is the matrix with elements D ij (t).In this work, Equations (5a) and (5b) will be used to predict the behaviour of the probability distribution of the random variable x in our control method.

Information Length (IL)
To obtain a minimum geometric information variability in the time evolution of the PDF of Equation ( 1), we need to investigate the concept of IL in more detail.In mathematical terms, given the joint PDF p(x; t) of the n-th order random variable ζ, we define its IL L as follows: The value 6) corresponds to the square of the Fisher information by considering the time as a parameter, and it is commonly called the information rate [21].The dimension of 1/Γ ≡ τ is time, which means that it serves as a dynamical time unit for information change.Hence, as shown in Figure 1, the integration of Γ between time 0 and t gives the total information change in that time interval; i.e., L quantifies the number of statistically different states that the system passes through in time from an initial p(x; 0) to a final p(x; t) [55].Furthermore, the IL is a model-free tool that proved to be a true metric between two PDFs in the statistical space [26].In [33], τ was shown to provide a universal bound on the timescale of transient dynamical fluctuations, independent of the physical constraints on the stochastic dynamics or their function.For Gaussian distributions, according to [25], we know that Γ can be rewritten in terms of µ and Σ as follows: Equation ( 7) simplifies the analysis of the IL when studying systems like (1).

Information-Thermodynamic Relation
However, while we recognise the utility of the IL as a tool for analysing the dynamics of time-varying PDFs, its physical significance remains unclear.Therefore, establishing a connection to a physical quantity is crucial.
In this context, consider the value of the entropy rate defined as follows: [56] Ṡ where Π is the entropy production due to irreversible processes occurring inside the system and Φ is the entropy flux from the system to the environment.Figure 2 shows a system subject to some boundary conditions to avoid matter exchange (i.e., a closed system).
The system produces entropy Π > 0 ∀t and exchanges entropy Φ with the environment (hence, it shares energy).Specifically, Φ > 0 (Φ < 0) when the entropy flows from the system (environment) to the environment (system).A system with a minimum entropy production Π produces the highest amount of free energy (useful work) [21,57].We refer the reader to [58,59] for a complete review of thermodynamics.In a first-order linear stochastic differential equation, the value of the information rate Γ can be related to Ṡ and Π as follows (for further details, see [50]): where D and Σ are the noise amplitude and the variance of the first-order Langevin equation.Equation ( 9) can be used to describe the physical effects of a minimum variability control in a dynamical system.For instance, as will be discussed later in Section 4, to obtain a minimum information variability while being out of the equilibrium (i.e., Π > 0), the control will exert entropy to the environment, creating a small but negative value for Ṡ.This means that a minimum information variability would not necessarily lead to a maximum free energy but to an optimal path where Π is limited by the value of Γ according to (9) 8).In a closed system, entropy rate Ṡ corresponds to the subtraction of the entropy produced by the system Π and the entropy exchanged with the environment Φ.

Minimum Information Variability Problem
Now that we understand the meaning of the IL, let us explain in more detail how the IL can be used to minimise deviations from the geodesic of the system's PDF time evolution.In [32], the authors use the inequality J (t F ) ≥ L(t F ) 2 where J = τ (Fisher divergence) with τ = t f − t 0 and L given by (6).As mentioned in Section 1, such inequality follows from the Cauchy-Schwartz inequality Γ 2 dt u 2 dt ≥ ( Γu dt) 2 with u = 1.But, most importantly, the equality holds for the minimum path where Γ is constant (see, e.g., [19,32]), and the deviation from this equality is said to quantify the amount of the disorder in an irreversible process [32].
From [20], such a statement can be clarified by the following procedure.Let us define the time-average for a function f (t Then, we can define the time-averaged variance Equation ( 10) describes an accumulative deviation from the geodesic connecting the initial and final distributions of the system dynamics.Thus, we can conclude that by setting Γ as a constant, we obtain a geodesic that defines a path where the process has the minimum geometric information variability.

BIBO Stability of the Linear Stochastic Process
As we address a control problem, it is essential to examine the bounded-input, bounded-output (BIBO) stability of (1).For this purpose, we revisit the BIBO stability conditions for Equations (5a) and (5b), which will be instrumental in our discussion of MPC stability.
Theorem 1.The mean (5a) and covariance (5b) dynamics of (1) are BIBO-stable if and only if the eigenvalues λ i of the matrix A satisfy the following inequality: where ℜ[s] stands for the real part of the complex value s ∈ C.
Proof.For a detailed proof of this result, please refer to [60,61].
Remark 1. Theorem 1 is considered to be satisfied throughout our examples; i.e., the control method is applied to stable systems only.Furthermore, the control actions are constrained to finite values.

Main Results
To guide the system's PDF time evolution along the geodesic of the IL while achieving a desired set point at time t = t f , we propose the following cost function: where In this work, we call Equation ( 12) the Information Length Quadratic Regulator (IL-QR).As will be discussed in Section 3, MPC enables us to obtain the solution of (12) via a numerical scheme, circumventing analytic complexities while being applicable to practical scenarios.Analytically finding the geodesic dynamics involves using the solution of the Euler-Lagrange equations of the IL.The steps of such an approach are discussed and successfully applied in [4] for a first-order stochastic differential equation.In Appendix A, we give the details of the procedure when considering a more generalised scenario.From (12), we observe that the first term on the right-hand side imposes a constant Γ 2 (needed to minimise the deviations from the geodesic [4]).The term involving Q guides the system towards a specified PDF defined by Y d .The third term on the right-hand side of (12) regularises the control action c to prevent abrupt changes in the inputs.Furthermore, Q and R are matrices that penalise the error ϵ = Y − Y d and the control input u, respectively.Additionally, I L penalises the square of the error (Γ 2 (t) − Γ 2 (0)), aiming to maintain Γ 2 (t) close to Γ 2 (0) over time.
In our approach, the control of the dynamics for µ is governed by u(t), while the dynamics of Σ(t) are adjusted by controlling the noise width using a time-dependent vector w(t), where its elements replace the non-zero constant values of the matrix D in (5b).As discussed in the numerical examples, the noise width can be altered by modifying a macroscopic observable such as temperature (for further details, refer to the Brownian motion models presented in [12]).

Model Predictive Control
As discussed in Section 1, the solution of our optimisation problem will be computed through the MPC method.Hence, the following discrete optimal control problem encoding the MPC formulation is required: In Equation (13), the symbol over Y and Γ refers to their predicted values over the influence of the control c throughout the optimisation process in the finite horizon of length N. Note that the value of Y is constrained by the discretised version of the set of Equations (5a) and (5b) given by where A d = I + T s A and B d = T s B if a first-order approximation of the time derivative considering the sample period T s is applied (we apply a fourth-order Runge Kutta instead of a first-order approximation to compute Y in our simulations).Note that we have added the argument c in Equation ( 13) when describing Equations ( 14) and ( 15) to emphasise the application of the control during the optimisation procedure.The initial conditions µ[0], Σ[0] of ( 14) and ( 15) change every time an optimal control c solution has been computed and they are subject to the measurements m, S of the real/simulated process.Every element c i of the control vector c is constrained by a lower and an upper limit denoted as c l,i and c u,i , respectively.Finally, f is the function describing the predicted value Γ 2 , which is defined as follows: For a clearer grasp of the MPC method's application in solving the IL-QR problem in real-world scenarios, Figure 3a,b show the control diagram and the functionality of the MPC's optimiser when considering a second-order stochastic system, respectively.Figure 3a illustrates the real-time operation of the MPC algorithm.As the process evolves, the MPC algorithm receives input parameters, including a given set point Γ 2 [0], Y d , the prediction values Ỹ of µ and Σ from a prediction model, a set of constrains, the cost function J N , and the current system dynamics Y, to solve the optimisation problem given in (13).Afterwards, the optimal control c solution of Equation ( 13) is applied to the system.The MPC method determines the optimal control c by utilising the differential equations of μ and Σ as a prediction model in a finite horizon of length N. Figure 3b provides a brief overview of the working principle of the MPC's optimiser block.
In the case of a stochastic process described by a bivariate time-varying PDF p with random variables x 1 and x 2 , the MPC optimiser method initiates the optimisation process using the measured system's PDF output (depicted in black).The optimisation process involves extrapolating the values of the PDF p within a finite horizon of length N and comparing them with the reference trajectory described by the PDF p d .The optimisation problem is solved using the interior point method with CasADi [44].In this study, the use of deterministic descriptions of the first two statistical moments over time has facilitated the control, prediction, and simulation of the PDF, owing to the nature of the Langevin equations under consideration.However, in scenarios involving pure data or more complex stochastic differential equations, estimating time-varying PDFs may require inference methods [62] or stochastic simulations [63].

c[k]
Gaussian Process x 1 x 2 x ∼ N(⟨x⟩, Σ) Reference Trajectory kT k using the interior point method with CasADi [44].(b) A diagram illustrating a discrete MPC scheme applied to a second-order stochastic process.The algorithm predicts the behaviour of the dynamical system within a finite horizon of length N and compares it with the reference trajectory described by the PDF p d .

Simulation Results
To demonstrate the numerical implementation of MPC for solving the IL-QR cost function (13), the subsequent subsections delve into applying the minimum information variability to both the Ornstein-Uhlenbeck process and the Kramers equation.Throughout the simulations, we utilise various parameters, which are conveniently summarised in Table 2.Each system undergoes a pair of simulations or "experiments", as detailed in Table 2.  First, let us consider the Ornstein-Uhlenbeck (O-U) process (see Figure 4) defined by the following Langevin equation: where ζ(t) is a random variable, u(t) is a deterministic force, γ is a damping constant, and ξ(t) is a short correlated random forcing such that ⟨ξ(t)ξ(t 1 )⟩ = 2Dδ(t − t 1 ) with D ≥ 0 and ⟨ξ(t)⟩ = 0.The results of the MPC implementation are shown in Figures 5-8.In the following, it is important to note that, in the figures, black colour is used to represent values labelled on the y-left axis, while blue colour is used for values labelled on the y-right axis, unless otherwise specified.Figure 5 illustrates the scenario where the desired state

It also shows the time evolution of the states Y(t) = [µ, β(t)]
and controls c(t) = [u(t), D(t)] ⊤ .From the results, we see that the method finds a geodesic motion (solution to the IL-QR) from the initial to the final state in less than 0.4 seconds.The geodesic motion is corroborated by the constant value of Γ 2 (t) ≈ Γ 2 (0) = 2.4 and the plot of the information length L whose shape is a line with a slope of 1.5526.Γ 2 (0) is computed by considering that u(t = 0) = D(t = 0) = 0.It is noteworthy that the value of Σ is observed to vary slightly over time compared to the hyperbolic analytical solution in [4], which is provided for a non-constant Σ (refer to Appendix A).
When analysing the stochastic thermodynamics of the O-U system controlled by the MPC technique, Figure 6 presents the plot of the entropy rate Ṡ in comparison with the entropy production Π, along with a plot of the value of Γ 2 using Equation ( 9).The analytical expressions for Ṡ and Π, along with their derivations, are given in Appendix C. In the closed-loop system, we can see that the MPC method induces slight changes in both the entropy production Π and the entropy rate Ṡ during the process.Since the values of Σ and µ in the desired state Y d are close to, and far from, their initial conditions at state Y(0), respectively, the balance between Ṡ, Π, and Γ 2 as given by ( 9) is kept by maintaining an almost constant D(t) and a u(t) with a nearly constant velocity.The second plot demonstrates that the balance expressed via Equation ( 9) is maintained as expected.The control generates a small entropy rate Ṡ converging to a negative value, indicating that most of the entropy flows from the system to the environment.The control induces entropy production to maintain the system out of equilibrium with respect to Y d .This departure from equilibrium is supported by the presence of entropy production Π > 0.

−→ t
Under conditions almost similar to those in Figure 5, Figure 7 shows the behaviour of the closed-loop system PDF, the states µ and Σ, and the controls D and u, as well as the behaviour of Γ 2 for Y d (t) = [1/30, 1/(2 × 3)] ⊤ .Here, Γ 2 (0) is computed by considering that u(t = 0) = 0 and D(t = 0) = 1/(2 × 0.3).The final state Y d is achieved at approximately t = 2.8.Once again, the geodesic behaviour is supported by the constant value of Γ 2 (t) ≈ Γ 2 (0) = 0.41 and the plot of the information length L, which depicts a line with a slope of 0.64759.9) demonstrates that to maintain a constant information rate Γ while being out of equilibrium, we increase the value of entropy production Π, and most of the entropy flows out to the environment, as indicated by the negative sign of Ṡ.

←− t
In comparison to the stochastic thermodynamics shown in Figure 6, Figure 8 exhibits minor changes in the entropy production Π and significant variations in the entropy rate Ṡ of the closed-loop system, given that the values of Σ and µ in the desired state Y d differ from its initial condition Y(0).This disparity also leads to slight fluctuations in both the D(t) and u(t).Additionally, Figure 8 demonstrates the validity of the balance in Equation (9).

Kramers Equation
To study the solution of the IL-QR problem in a higher-order system via the MPC method, let us consider the so-called Kramers equation.The Kramers equation describes the Brownian motion of particles in an external field [12].The non-autonomous version of the Kramers equation is given by the following set of Langevin equations: Here, ξ is a short correlated Gaussian noise with a zero mean ⟨ξ⟩ = 0 and a strength D ≥ 0 with the following property: In practice, as shown in Figure 9, the Kramers Equation ( 18) is also a good first approximation to describe the dynamics in one-dimension of a particle in an optical trap [5].The Kramers equation control c and state Y vectors are defined by  18) as a mass-spring-damper system.In this system, the external force u(t) is deterministic, while ξ(t) represents a stochastic perturbation on the process, which varies due to the temperature of the fluid [12].In addition, they include the time evolution of the bivariate PDF p(x; t) with the spatial vector x = [x 1 , x 2 ] ⊤ , with x 1 and x 2 representing the position and velocity of the system dynamics, respectively.The remaining parameters used in the simulations are listed in Table 2. Again, in the figures, black colour is used to represent values labelled on the y-left axis, while blue colour is used for values labelled on the y-right axis, unless otherwise specified.
For the first numerical experiment, Figure 10 demonstrates that the MPC method is capable of maintaining Γ 2 as constant through time while reaching the desired state Y d .Here, the value of Γ(0) in ( 13) is obtained as follows: Tr Σ −1 (0)(AΣ(0) + Σ(0)A ⊤ + 2D(0)) 2 = 6.16667, (22) where u(0) = 0 and vec(D(0)) = [0, 0, 0, 1/(2 × 0.3)] ⊤ , while µ(0), Σ(0) and A are taken from the corresponding Y(0) and the mathematical model (18), respectively.The geodesic dynamics give a behaviour with slow oscillations in the state Y.The controls u and D present high oscillations as the system reaches the desired state Y d .The system reaches Y d at t ≈ 7 with an error of 1 × 10 −3 .The geodesic behaviour is supported by the linear behaviour of the information length L compared to the fitted equation y = 24.8332t.).The control effectively adjusts the values of the covariance matrix Σ and the mean vector µ by dynamically modifying u and D, transitioning the system's PDF from one state to another out of equilibrium.It is noteworthy that, to maintain the system on the IL's geodesic, the MPC method maintains a constant Γ 2 , resulting in abrupt changes in u and D as the system approaches the desired state Y d .Physically, this leads to a high entropy production due to the control action.
Figure 10 shows a second numerical experiment where Y d is even farther from the system's equilibrium.Yet, the MPC method can maintain Γ 2 as constant through time while reaching Y d .Like the example of Figure 10, in this case Γ 2 (0) = 6.16667.Small oscillations persists in the time evolution of µ 1 , µ 2 , Σ 11 , cΣ 12 , and Σ 22 .The system reaches the desired state Y d at t ≈ 8.5.Thus, the farther the desired state Y d is from the initial state Y, the longer the time it takes to reach it while following the geodesic path.The geodesic behaviour is evidenced by the plot of the information length L, whose behaviour is compared to the fitted equation y = 24.8336t.Similar to the example in Figure 10, the controls exhibit highly oscillatory behaviour as the system reaches Y d .

Conclusions
In this work, we demonstrated the application of the MPC method, an online optimisation algorithm for constrained control problems, to achieve the minimum information variability in systems governed by linear stochastic differential equations.The linearity of the system results in time-varying Gaussian PDFs, with statistical moments determined by a set of deterministic differential equations.Our simulations demonstrate that the MPC method is a practical approach for determining the geodesic of the information length between the initial and desired probabilistic states through the solution of the proposed IL-QR loss function.
From a thermodynamic standpoint, simulations of the MPC in the O-U process reveal that the MPC directly influences the amount of entropy production generated by the system to meet all optimisation requirements.Future work will involve extending this approach to nonlinear stochastic dynamics, such as toy models [64], systems with uncertain physical parameters [65], Brownian motion [66], or diffusion [67,68].This extension aims to maximise the free energy [21] by minimising the entropy production and to analyse closedloop stochastic thermodynamics for higher-order systems.It is worth noting that extending the method to nonlinear dynamics can be achieved through the Laplace assumption [69] and/or by employing Kalman/particle filter methods [70,71].
where Q = sinh −1 ( ϕ 2 ) using ϕ = β(0)(µ(0) − µ(t F )). Clearly, through Equations (A5) and (A6) and the dynamical model of the O-U process (17), we can construct the optimal control f (t) of the input u(t) and the optimal noise amplitude D I (t) of D(t).From [4], given that u(0) = 0, such optimal controls are given by f (t) = µ(t) − β(0)µ(0) β(t) , (A8) To compare the analytical and MPC solutions of the geodesic of the IL, Figure A1 shows the behaviour of the O-U process when controlled through the analytical solutions (A8) and (A9) or the IL-MPC method.The figure contains different subplots that show the time evolution of µ, β 1/2 , Γ 2 , and L and the optimal controls D I and f .In the simulation, the desired state and the damping are Y d (t) = [1/30, 1/(2 × 0.3)] ⊤ and γ = 1, respectively.Additionally, we set a fixed final time t F = 2A = 0.9304 (one cycle of the hyperbolic geodesic motion (A5) and (A6)) by considering the initial state Y(0) = [5/6, 1/(2 × 0.3)] ⊤ .Figure A1 uses dashed and non-dashed lines for the MPC and the analytical response, respectively.Additionally, recall that black colour is used to represent values labelled on the y-left axis, while blue colour is used for values labelled on the y-right axis.From the comparison, a major conclusion is that the time evolution of β is no longer hyperbolic when using the MPC method.This means that the MPC method finds an almost constant Σ solution but not the hyperbolic solution shown in [4].The MPC allows us to reach the final state Y d at t F with an error of 6.6 × 10 −4 .As a second example, Figure A2 shows the dynamics of the controlled O-U process when the initial state is Y(0) = [5/6, 1/(2 × 3)] ⊤ (fixing t F = 2A = 0.7367), the desired state Y d (t) = [1/30, 1/(2 × 3)] ⊤ , and the damping γ = 1.Again, the MPC method recovers a geodesic solution where the β time evolution is constant.In this scenario, the MPC method reaches the desired state Y d with an error of 9.8 × 10 −4 in a time t > t F , demonstrating that the numerical optimisation scheme may not recover an optimal time.As a final remark, note that if the n-th order case is considered, Equations (A3) and (A4) form a set of non-linear differential equations whose solution may be obtained by a numerical procedure.But, even for the case of a second-order stochastic process, this becomes a challenging problem (we have a boundary value problem of 12 non-linear differential equations).Hence, the MPC method provides an alternative solution to this problem while being an experimentally feasible approach, as demonstrated by the application to the Kramers equation in Section 4.2.

Figure 1 .
Figure 1.Graphical description of the value of L. From information geometry, L quantifies the number of statistically different states that the system passes through in time from an initial p(x; 0) to a final p(x; t). The complete derivation of how to compute the values of Ṡ, Π, and Φ for a first-order Langevin equation is shown in Appendix C.

Figure 2 .
Figure 2. Graphical description of Equation (8).In a closed system, entropy rate Ṡ corresponds to the subtraction of the entropy produced by the system Π and the entropy exchanged with the environment Φ.
2 is the vector of states µ and vec(Σ) that define the time evolution of p(x; t);Y d = [µ d , vec(Σ d )] ⊤ isthe desired position of the n + n 2 states defined by µ d ; and Σ d at time t, and c ∈ R m (such that m ≤ 1 + n 2 ) is the vector of controls defined by c = [c 1 , c 2 , . . ., c m ]

Figure 3 .
(a) A control diagram describing the main components of the implemented MPC methodology.The algorithm comprises a prediction model utilised to determine the optimal control c[k]

Figure 4 .
Figure 4.The O-U process equation is commonly used to describe a prototype of a noisy relaxation process, for instance, the movement of a particle confined to a harmonic potential V(ζ) = 1 2 γ(ζ − u(t)) 2 and thermal fluctuations with temperature T (D = k B Tγ, and k B is the Boltzmann constant) such that ζ(t) fluctuates stochastically.

Figure 5 .
Figure 5. Minimum information variability control of the O-U process (Experiment 1).In this experiment, the initial and desired PDFs have similar variances but different mean values.In the simulation, Γ 2 is maintained constant at a value of Γ 2 = 2.4.To follow the geodesic, the input force u(t) needs to linearly decrease, while the noise amplitude D(t) decreases slightly exponentially.Conversely, the mean value undergoes a linear change, while the variance exhibits hyperbolic behaviour.

Figure 6 .
Figure 6.Stochastic thermodynamics under minimum information variability for the O-U process (Experiment 1).The second plot demonstrates that the balance expressed via Equation (9) is maintained as expected.The control generates a small entropy rate Ṡ converging to a negative value, indicating that most of the entropy flows from the system to the environment.The control induces entropy production to maintain the system out of equilibrium with respect to Y d .This departure from equilibrium is supported by the presence of entropy production Π > 0.

Figure 7 .Figure 8 .
Figure 7. Minimum information variability control of the O-U process (Experiment 2).In this experiment, the initial and desired PDFs have different variances and mean values.Γ 2 is maintained constant at a value of Γ 2 = 0.41, while both the input force u(t) and the noise amplitude D(t) decrease exponentially.

Figure 9 .
Figure 9. Mechanical representation of the Kramers Equation (18) as a mass-spring-damper system.In this system, the external force u(t) is deterministic, while ξ(t) represents a stochastic perturbation on the process, which varies due to the temperature of the fluid[12].

Figures 10 and 11
Figures 10 and 11 show the simulation results of the closed-loop Kramers equation when considering the desired states Y d(t) = [0, 0, 1/(2 × 3), 0, 1/(2 × 3)] ⊤ and Y d (t) = [−5/6, 0, 1/(2 × 3), 0, 1/(2 × 3)] ⊤ ,respectively.These figures include the time evolution plots of the mean values µ 1 and µ 2 and the covariance matrix values Σ 11 , Σ 12 , and Σ 22 of the Kramers equation random variables ζ 1 and ζ 2 .In addition, they include the time evolution of the bivariate PDF p(x; t) with the spatial vector x = [x 1 , x 2 ] ⊤ , with x 1 and x 2 representing the position and velocity of the system dynamics, respectively.The remaining parameters used in the simulations are listed in Table2.Again, in the figures, black colour is used to represent values labelled on the y-left axis, while blue colour is used for values labelled on the y-right axis, unless otherwise specified.For the first numerical experiment, Figure10demonstrates that the MPC method is capable of maintaining Γ 2 as constant through time while reaching the desired state Y d .Here, the value of Γ(0) in (13) is obtained as follows:

Figure 10 .
Figure10.Minimum information variability control of the Kramers equation (Experiment 1).The control effectively adjusts the values of the covariance matrix Σ and the mean vector µ by dynamically modifying u and D, transitioning the system's PDF from one state to another out of equilibrium.It is noteworthy that, to maintain the system on the IL's geodesic, the MPC method maintains a constant Γ 2 , resulting in abrupt changes in u and D as the system approaches the desired state Y d .Physically,

Figure 11 .
Figure 11.Minimum information variability control of the Kramers equation (Experiment 2).In this scenario, the value of Y d results in slower dynamics compared to Experiment 1, as MPC now requires more time to reach Y d .Once again, significant changes are observed in u and D to maintain a constant information rate Γ when the system approaches the desired state Y d .

Table 1 .
Description of the different variables used throughout this work.

Table 2 .
Parameters used in the simulation results for the O-U and Kramers systems.The table includes figure numbers showing the simulation where the set of parameters was used.