Optimal Control in Stochastic Thermodynamics

We review recent progress in optimal control in stochastic thermodynamics. Theoretical advances provide in-depth insight into minimum-dissipation control with either full or limited (parametric) control, and spanning the limits from slow to fast driving and from weak to strong driving. Known exact solutions give a window into the properties of minimum-dissipation control, which are reproduced by approximate methods in the relevant limits. Connections between optimal-transport theory and minimum-dissipation protocols under full control give deep insight into the properties of optimal control and place bounds on the dissipation of thermodynamic processes. Since minimum-dissipation protocols are relatively well understood and advanced approximation methods and numerical techniques for estimating minimum-dissipation protocols have been developed, now is an opportune time for application to chemical and biological systems.


I. INTRODUCTION
In this article we review recent progress in optimal control of stochastic thermodynamic systems. We focus on classical isothermal stochastic thermodynamics describing control by linear-response, thermodynamic geometry, and optimal-transport theory.
Historically, modern thermodynamic control began with the study of finite-time thermodynamics of macroscopic systems, 1-3 the natural extension beyond quasistatic (infinitely slow) processes. Fundamentally, any finite-time thermodynamic control will induce some degree of irreversibility, manifesting as energy dissipated into the environment. A goal of finite-time thermodynamics is to quantify and minimize this dissipation through the use of designed control strategies. For example, Ref. 4 studied the optimal cycle for finite-time operation of a heat engine and found that instantaneous jumps in control parameters are necessary to minimize dissipation.
In parallel, a thermodynamic-geometry framework was developed to provide a novel means to describe thermodynamic processes on a smooth (generally Riemannian) manifold. [5][6][7] Ref. 8 showed the connections between thermodynamic geometry and minimum-dissipation protocols, opening the door for the development of a geometric description of minimum-dissipation protocols. Although theoretically compelling, the utility of the framework was not fully realized until the development of stochastic thermodynamics.
The aforementioned descriptions focused on macroscopic systems that equilibrate rapidly and whose fluctuations are relatively small. The advent of modern experimental techniques, including single-molecule biophysical experiments, created demand for a theoretical description of the energetics of microscopic systems. Optical tweezers and magnetic traps 9-12 allow for the precise manipulation of individual polymer strands [13][14][15][16][17][18] and nanoscale molecular machines 19 (e.g., ATP synthase, [20][21][22] kinesin, [23][24][25][26] and myosin [27][28][29][30] ). Due to their small scale, these systems are not accurately described by macroscopic thermodynamics since the fluctuations (of order k B T ) are comparable to the systems' internal energy scales, and the operating speeds of single-molecule experiments and molecular machines are comparable to those systems' natural relaxation times.
To describe these small-scale systems, the field of stochastic thermodynamics was developed, which aims to describe the nonequilibrium energetics of stochastic (fluctuating) microscopic systems 31,32 . Just like its macroscopic counterpart, a central goal of stochastic thermodynamics is the description of optimal control strategies: methods for performing a given task at minimum energetic cost 33,34 .
There are two distinct but related types of control we consider in this review: full control (section IV)and parametric control (section V). Full control assumes we have complete control of the probability distribution ( Fig. 1, top). Parametric control adjusts a finite number of control parameters (Fig. 1, bottom), and in doing so drives the probability distribution.
For either full or parametric control, the exact minimum-dissipation protocol is known if the probability distribution is Gaussian [35][36][37] . These exact solutions provide a glimpse into the properties of optimal control processes. For example, just like the finite-time thermodynamic control described previously, the minimum-dissipation protocol has discontinuous changes in control parameters at the start and end of the protocol but remains continuous between these endpoints 35 . These discontinuities are present even for underdamped dynamics 38 . The control-parameter jumps have been observed in a number of different systems [38][39][40] and are now well understood and have been shown to be a general feature 41 .
For more general solutions under full control, the study of minimum-dissipation protocols can be mapped onto a problem of optimal-transport theory, a well developed branch of mathematics for which there exist numerous algorithms and methods for determining the optimal-transport map 42,43 . The connection between minimum-dissipation arXiv:2212.00706v2 [cond-mat.stat-mech] 11 Apr 2023 Full control (top) assumes complete control of the probability distribution p(r, t) (shaded) which can be optimally driven between the endpoints by a potential V (r, t) (red dashed curves). Parametric control (bottom) adjusts a finite number of control parameters λ(t) according to a protocol Λ between specified endpoints, thereby driving the probability distribution p(r, Λ).
protocols and optimal-transport theory was first shown in Ref. 44 for overdamped dynamics: the protocol that minimizes dissipation when driving a system obeying overdamped Fokker-Planck dynamics between specified initial and final distributions is governed by the Wasserstein distance [45][46][47][48] and the Benamou-Brenier formula 49 . This technique eventually led to new fundamental lower bounds on the average work required for finite-time information erasure 50,51 . Initially only applicable to overdamped dynamics, the connections between optimal-transport theory and minimumdissipation protocols have recently been shown for discrete-state and quantum systems 47,52-56 . General solutions for parametric control are typically difficult to determine, although recent progress has been made towards exact solutions for general systems building off of optimal-transport 54 or advanced numerical techniques 39,57,58 . Although exact solutions are convenient where possible, the determination of minimum-dissipation protocols can be considerably simplified through approximate methods. Inspired by a diagram presented in Ref. 59, we schematically show in Fig. 2 the limits where minimum-dissipation protocols are known.
Linear-response theory can be used to determine the minimum-dissipation protocol for weak perturbations and performs relatively well at any driving speed and beyond its strict range of validity 59,60 . For slow control, the thermodynamic-geometry framework has been generalized to stochastic thermodynamic systems 7,61 and has been used to explore a diverse set of model systems [61][62][63][64][65][66][67][68][69][70][71][72][73][74] , including DNA-pulling experiments 75 and free-energy estimation 76 . In the opposite limit of fast control, minimum-dissipation protocols are describe by short-time efficient protocols 41 , which can be combined with the thermodynamic-geometry framework to design interpolated protocols that perform well at any driving speed 77 . Leveraging known solutions from optimal-transport theory, strong control can be described by the strong-trap approximation, yielding explicit solutions for minimum-dissipation protocols 37 .
This review is organized as follows: we begin with examples of both experimental and theoretical model systems in section II, followed by a brief introduction to stochastic thermodynamics of heat, work, and entropy production in The space of thermodynamic control. Horizontal axis is the driving speed from slow to fast, and the vertical axis is the strength of driving from weak to strong. Linear-response theory is applicable to weak and slow driving (blue), and can be simplified to a thermodynamic-geometry framework for slow driving. Short-time efficient protocols are valid for fast driving (purple) and can be combined with thermodynamic geometry to bridge the space between slow and fast with interpolated protocols (green). The strong-trap approximation (red) is only valid for overdamped dynamics, with region of applicability schematically indicated by a distinct dotted line. Exact solutions for Gaussian distributions serve as a window into the properties of minimum-dissipation protocols and are valid at any driving speed or strength of driving. section III. Section IV reviews recent progress exploiting optimal-transport theory to determine minimum-dissipation protocols under full control, yielding explicit solutions for the minimum-dissipation protocol under a strong-trap approximation in section IV B and allowing for constrained final control parameters in section IV C. Section V reviews parametric control, focusing on approximation methods in the fast (section V A), weak (section V B), and slow (section V C) limits. Applications to free-energy estimation are discussed in section VI before comparing the performance of designed protocols in section VII and finally concluding in section VIII with a perspective and outlook for the study of optimal control in stochastic thermodynamics.

II. MODEL SYSTEMS
In this section we provide a brief introduction to a few paradigmatic model systems that motivate and guide the study of optimal control in stochastic thermodynamics. As discussed in section I, the growth of stochastic thermodynamics coincides with the advent of new experimental techniques used to manipulate and measure single-molecule biophysical systems [20][21][22][23][24][25][26][27][28][29][30] . First and foremost among these techniques are laser optical tweezers 10-12 which can be used to trap microscopic Brownian systems.
The simplest experimental apparatus for studying stochastic thermodynamics is that of a microscopic bead trapped in an optical potential. From a theoretical perspective, this system is well approximated by continuous overdamped Brownian motion in a quadratic constraining potential. In these experiments, the center and stiffness of the trapping potential can be dynamically controlled to  apparatus can be augmented to realize a virtual constraining potential of any form 102-104 and can, for example, be used to study fundamental bounds on information processing through bit erasure 105,106 . Microscopic beads trapped by laser optical tweezers can be attached to biopolymers to probe their properties. For example, dual-trap optical tweezers can be used to fold and unfold DNA or RNA hairpins by modulating the separation between the trapping potentials ( Fig. 3 a). [13][14][15][16][17][18] Monitoring the position of the probe beads provides insight into the properties of the indirectly observed biopolymers. The simplest model representing this process is that of a driven barrier crossing 107 , where a Brownian system is dynamically driven over an energy barrier by a time-varying quadratic trapping potential (Fig. 3 d). 61,77 Magnetic traps can be used to probe the F1 component of the rotary motor ATP synthase ( Fig. 3 b), [20][21][22] which is driven periodically to synthesize ATP, an essential and portable energy source for the cell. Once again, microscopic beads are used to probe the properties of the molecular machine and can be dynamically driven (Fig. 3 e); however, the control differs from driven barrier crossing in that the driving is periodic.
As a final example, consider a nanomagnetic bit characterized by its spin state or average magnetization (Fig. 3 c). By applying an external magnetic field, the system state can be driven from all spin-down to all spin-up, reversing the magnetization and resulting in a bit flip 108 . This type of system is typically modeled with a discrete state space, e.g., the Ising model ( Fig. 3 f), and optimal control of this system has been investigated [70][71][72] . Due to the discrete state space, the properties of optimal control can differ from those for a system with a continuous state space.
Throughout this review we will use the model system of overdamped driven barrier crossing previously studied in Ref. 77 in order to give some intuition and examples for optimal control. The model consists of an overdamped Brownian particle in a double-well potential (symmetric for simplicity) constrained and driven by a quadratic trapping potential (Fig. 3 d). This system represents a simplified model of hairpin pulling experiments and Landauer erasure 50,51,77 . The two-state nature of the system is also representative of activated processes such as chemical reactions, and the barrier crossing mimics the main features of experiments performed on ATP synthase, whose dynamics can be approximated as a series of barrier crossings 22,69,109 . This model is also typical of steered molecular-dynamics simulations, which use a time-dependent quadratic potential to drive reactions [110][111][112] .
is the sum of the static hairpin potential V land [x] and time-dependent trap potential V trap [x, x c (t), k(t)] (shown schematically in Fig. 3d). The hairpin potential is modeled as a static double well (symmetric for simplicity) with the two minima at x = 0 and x = ∆x m representing the folded and unfolded states 17,18,61,107 , for barrier height E B , distance x m from the minimum to barrier, and distance ∆x m = 2x m between the minima. The system is driven by a quadratic trap with time-dependent stiffness k(t) and center x c (t).

III. THERMODYNAMICS
In this review we focus on thermodynamics at the distribution level, which is generally described by dynamics of the form governing the time evolution of a nonequilibrium probability distribution p(r, t) over position vector r at time t according to the time-evolution operator L[r, t].
Continuous-space stochastic systems are described by the Fokker-Planck equation, which for overdamped dynamics has the time-evolution operator for mean local velocity 46 v(r, t) ≡ −D∇ [βV tot (r, t) + ln p(r, t)] , total potential V tot , diffusivity D, and β ≡ (k B T ) −1 for temperature T and Boltzmann's constant k B . For a discretestate system, L[r, t] is the transition rate matrix. For the example of driven barrier crossing, the total potential includes both the hairpin and trapping potential, and the time dependence arises from dynamic changes in the trap center and stiffness which drive the system over the barrier.
The average system energy is and the rate of change in energy is The first term quantifies work W done on the system by an external agent controlling the potential V (r, t) (e.g., an experimentalist dynamically driving a trapping potential) and the second quantifies heat flow Q into the system from changes in the system distribution p(r, t) (e.g., the system responding and relaxing towards a new equilibrium distribution in response to movement of the center and stiffness of the trapping potential): Throughout, a dot above a variable denotes the rate of change with respect to time. Substituting (8) and (9) into (7) gives the first law of thermodynamics: any change in system energy equals work and heat flows into the system. Optimal control in thermodynamics is often discussed in terms of minimizing either entropy production or excess work incurred during the protocol. The entropy production as defined in this review will be used for periodic systems (e.g., ATP synthase), and when we have full control over the distribution (section IV). The excess work is used for parametric control (section V) and model systems like DNA hairpins and nanomagnetic bits which are driven between control-parameter endpoints rather than periodically.
To understand these two concepts, consider the nonequilibrium free energy for dimensionless entropy If the probabilities in (10b) are equilibrium distributions, then the nonequilibrium free energy reduces to the equilibrium free energy F eq . For an isothermal process the rate of change in free energy is The dimensionless entropy production isṠ for environmental entropy S env . Equation (13c) follows from the second law of thermodynamics. Importantly, this definition of entropy production only accounts for dissipation during the protocol, so if the system is not in equilibrium at the end of the protocol, additional energy may be dissipated into the environment during subsequent relaxation. A measure of energy dissipation which accounts for relaxation to equilibrium even after the protocol terminates is the excess work W ex ≡ W − ∆F eq , the amount of work done in excess of the equilibrium free-energy difference. The excess work and entropy production are related by for net entropy production ∆S prod ≡ ∆t 0 dtṠ prod , nonequilibrium free-energy difference ∆F neq between the initial p(r, 0) and final p(r, ∆t) distributions, and equilibrium free-energy difference ∆F eq between the initial and final equilibrium distributions.
Equation (14) clarifies the distinction between entropy production and excess work: if both the initial and final states of the system are at equilibrium, the two quantities are equal, otherwise the difference between the two equals the difference between the nonequilibrium and equilibrium free-energy changes. If the system is allowed to relax to equilibrium, then this excess energy is dissipated into the environment, resulting in additional entropy production not accounted for in the present definition (13a), so that the total entropy production for such a process is (13a) plus the entropy production from the subsequent relaxation. In contrast, the excess work always includes the energy dissipated into the environment from the system relaxing towards equilibrium after the protocol terminates. Essentially, quantifying dissipation by the entropy production in (13a) assumes one can harness the nonequilibrium free energy at the conclusion of the protocol to perform a useful task (generally true for periodically driven systems like ATP synthase), while excess work quantifies dissipation when all the excess free energy is dissipated into the environment after the protocol terminates (generally true for two-state barrier crossings like hairpin experiments).

IV. FULL CONTROL
For continuous-state systems, complete control over the shape of the potential V (r, t) grants full control over the probability distribution p(r, t), which can considerably simplify the optimization process and allows us to exploit known results from optimal-transport theory 44,46,113 . Optimal transport describes the most efficient methods to move mass (e.g., a pile of sand) from one location to another; this is useful for describing methods that minimize dissipation in transporting probability from an initial to a final distribution 42 .
Since the final distribution is constrained, the energy dissipated into the environment throughout the protocol is determined by the average entropy produced in driving from initial probability distribution p(r, 0) to final probability distribution p(r, ∆t) as 46 Angle brackets · · · denote an average over p(r, t).
Expressing the entropy production in this form makes precise the connections between optimal-transport theory and minimum-dissipation protocols. In optimal-transport theory a common measure of the distance between two distributions is the L 2 -Wasserstein distance, defined in the Benamou and Brenier dual representation as 49 Therefore, the entropy production is bounded by the squared L 2 -Wasserstein distance between initial (p 0 ) and final (p ∆t ) probability distributions 46 : This allows us to exploit existing procedures from optimal-transport theory to determine protocols that minimize entropy production 42,43 . Notably, the exact solution is known in two situations: one-dimensional systems and Gaussian probability distributions (section IV A). Extending minimum-dissipation full control and the connections to optimal-transport theory to more general forms of dynamics (e.g., discrete state spaces) is a rapidly advancing area of active research. 47,52-56

A. Exact solutions
For a one-dimensional system r = x, the entropy production (15) of the optimal-transport process can be simplified considerably as 44,45,50,51,80,114 where Q f and Q i are the final and initial quantile functions (inverse cumulative distribution functions) 77 . The entropy production is minimized if the quantiles are linearly driven between their fixed initial and final values 45,50,51,80 ; from this the probability distribution can be computed, and then the Fokker-Planck equation inverted to determine the potential V tot (x, t) to be applied to achieve the control that minimizes the entropy production: Although this calculation is often analytically intractable, it is straightforward to compute numerically for any probability distribution.
If the initial and final distributions are Gaussian, p(r, t) = N (µ t , Σ t ) for time-dependent mean µ t and covariance Σ t , the entropy-production bound (17) is 36,115,116 with subscripts 0, t, and ∆t respectively denoting the initial, time-dependent, and final values of the corresponding variable. Equality is achieved and the entropy production is minimized when following the optimal-transport map between the initial and final distributions, which for Gaussian distributions is completely specified by the mean and covariance: Here ∆µ ≡ µ ∆t − µ 0 is the total change in mean position, I is the identity matrix, and C ≡ Σ reduces in 1D to the ratio of final and initial standard deviations. If the covariance matrix Σ is diagonal, then (21b) implies with ∆Σ Thus for diagonal covariance the optimal-transport process linearly drives the standard deviation between its endpoint values. For a detailed description of minimum-dissipation protocols for general multidimensional Gaussian distributions see Ref. 36.
In both solvable cases (1D and Gaussian) the general design principle is that the minimum-dissipation protocol linearly drives the quantiles between specified initial and final values. Linearly driving the quantiles of the probability distribution can be used as a guiding principle for designing more general minimum-dissipation protocols and arises independently for parametric control 77 .

B. Strong-trap approximation
In this section we review how to exploit the known full-control solution for Gaussian distributions in order to design minimum-dissipation protocols for strong trapping potentials on any arbitrary underlying energy landscape, initially described in Ref. 37.
We assume that the total potential V tot [r, r c (t), can be separated into a time-independent component V land (r) (the underlying energy landscape) and a quadratic trapping potential The position of the system is denoted by the vector r, K is the symmetric stiffness matrix, and superscript is the vector transpose. For a strong trapping potential, the time-independent component can be expanded up to second order about the mean particle position µ: Under these assumptions, the probability distribution can be approximated as Gaussian, p(r, t) ≈ N (µ t , Σ t ). Since the distribution is Gaussian we can use the results described by (21a) and (21b). To achieve the mean and covariance of (21a) and (21b), the trap center and stiffness must respectively satisfy If Σ is diagonal, then Σ t is given by (22), the integral in (25b) can be evaluated, and the trap stiffness obeys Equations (25a) and (25b) are explicit equations that minimize dissipation at any driving speed provided the trap is sufficiently strong. We emphasize that explicit solutions for the minimum-dissipation protocol are rarities; typically solutions require numerically solving differential equations or inverting the Fokker-Planck equation 44,50,51 . For our example system of driven barrier crossing with equal initial and final covariance, the protocol designed by (25a) slows down movement of the trap center as it crosses the barrier to compensate for the force due to the hairpin potential, and (27) tightens the trap as it crosses the barrier in order to counteract the negative curvature of the hairpin potential. Slowing down and tightening the trapping potential as the system crosses a barrier appears to be a general feature and results independently from other approximation methods 37,77 .
To achieve periodic driving we constrain the final covariance matrix after one period (during which the mean completes one rotation) to equal the initial. To minimize dissipation, the standard deviation is linearly changed between the endpoints (21b), implying for a periodic system that the standard deviation (and hence covariance) remains unchanged throughout the protocol. This is achieved when the effective stiffness is constant, i.e., If in each rotation the mean travels a distance ∆µ rot in time ∆t rot , the resultant entropy production is that of an overdamped system moving at constant velocity against viscous Stokes drag; i.e., the minimum-dissipation protocol has perfect Stokes efficiency 117 .

C. Constrained final control parameters
The techniques described in section IV minimize the entropy production for constrained final distributions. Constraining the final distribution allows us to describe periodically driven systems like models inspired by ATP synthase; however, often the end state is instead constrained by the final control-parameter values. For example, in hairpin experiments the end state is typically defined by the trap separation, and the system is typically allowed to equilibrate between subsequent unfolding/folding protocols. Since the system is allowed to equilibrate between protocols, any excess energy remaining at the end of the protocol is dissipated as heat into the environment, and the entropy production as defined by (13a) does not equal the total dissipation.
In this case the relevant thermodynamic quantity is the work (14), which accounts for the excess free energy dissipated into the environment after the protocol terminates. Work can be minimized by optimizing over the final distribution subject to constrained final control-parameter values 37,77 . In general this minimization must be performed numerically, which is straightforward for one-dimensional systems, but becomes computationally intensive for higher-dimensional systems.
For the special case of Gaussian distributions, the optimization procedure simplifies considerably since optimizing over a distribution reduces to optimizing over means and covariances. The average work for the minimum-dissipation protocol in the strong-trap approximation is 37 with Tr the trace and ∆S min prod the lower bound in (20). To find the protocol that minimizes work for constrained final control parameters, we minimize (30) with respect to the final mean µ ∆t and covariance Σ ∆t , for fixed final trap center λ ∆t and stiffness K ∆t . For a flat energy landscape (for which the strong-trap approximation is exact) the optimal mean and covariance can be solved analytically; e.g., for equal initial and final covariance, the final mean is In some more general cases (e.g., energy landscapes represented by low-order polynomials), (30) can also be minimized analytically, and in general can be solved numerically with relative ease.

V. PARAMETRIC CONTROL
While full-control solutions are convenient when possible, many applications do not permit sufficient control to fully constrain the probability distribution throughout the entire protocol. In such cases the controller is constrained by a finite set of control parameters λ(t) which can be used to drive the system. Since there are insufficient control parameters to fully control the probability distribution, the endpoints of the protocol are constrained by the controlparameter values rather than the distribution. Therefore, we focus on optimizing the excess work, which is an accurate measure of dissipation provided the system equilibrates after the protocol terminates.
In this section we assume the control is the result of time-dependent control parameters λ(t), in which case the excess work (8) is where we have defined the conjugate force f ≡ −∂V tot /∂λ to express the work explicitly in terms of the control parameters, and δ denotes a difference from the equilibrium average. In this section we denote a nonequilibrium average (dependent on the entire protocol history) as · Λ and an equilibrium average (at the current control parameter value λ(t)) as · λ(t) . Since the control-parameter protocol is externally specified, the only unknown is δf (t) Λ . This quantity is particularly difficult to deal with since the nonequilibrium average depends on the entire protocol history.
For example, Ref. 35 showed that-even in one dimension-optimization requires solving the nonlocal Euler-Lagrange equation. Therefore, there are very few cases where the minimum-work protocol can be determined analytically. Promising approaches for obtaining general solutions include optimal-transport theory with limited control 54 and advanced numerical techniques 39,57,58 .
The minimum-work protocol can be determined exactly for a Brownian particle in a harmonic trap, in which case the optimal protocol, originally derived in one dimension 35 , is identical to the full-control solution for Gaussian distributions described in section IV C. This exact solution serves as a window into the properties of minimum-dissipation protocols and gives us considerable insight into what to expect from optimized protocols: e.g., the minimum-dissipation protocols have control-parameter jumps at the start and end but remain continuous in between.
Although exact solutions are nice when possible, since general solutions are intractable we turn to approximate methods to gain insight into the general properties of minimum-dissipation protocols. The first approximation has already been described in section IV B and is valid for strong trapping potentials applied to overdamped systems 37 . In the following sections, we fill in the remaining limits of fast, weak, and slow control in order to describe optimal control for any system with any given control parameters and design interpolated protocols.

A. Fast control
Until recently, although it was generally suspected, it was not definitively known whether the jumps in the minimumdissipation protocols are a general feature. By taking the fast-driving limit, it has been shown that the minimumdissipation protocol is a step function, jumping from and to specified initial and final control-parameter endpoints to spend the entire intervening protocol duration at the control-parameter values that maximize the short-time power savings. 41 We refer to the minimum-dissipation protocols in the fast limit as short-time efficient protocols, or STEP.
In the fast limit, the excess work approaches that of an instantaneous protocol, which spends no time relaxing towards equilibrium and requires excess work proportional (up to a factor of β) to the relative entropy D(π i ||π f ) between the respective initial and final equilibrium distributions, π i and π f 41 . Spending a short duration ∆t relaxing towards equilibrium throughout the protocol results in saved work W save ≡ β −1 D(π i ||π f ) − W ex compared to an instantaneous protocol, which can be approximated as 41 in terms of the initial force-relaxation rate the rate of change of the initial mean conjugate forces at the current control-parameter values.
The saved work is maximized by the short-time efficient protocol (STEP) which spends the entire duration at the intermediate control-parameter value that maximizes the short-time power savings satisfying The STEP achieves this by two instantaneous control-parameter jumps: one at the start from the initial value to the optimal value λ STEP , and one at the end from λ STEP to the final value. The optimal protocol described in this section is general, only assuming the protocol is fast compared to the system's natural relaxation time. Additionally, for quadratic time-dependent potentials (such as the driven barrier-crossing model system) or affine control (V tot (r, λ) = V 0 (r) + λV 1 (r) where λ linearly modulates the strength of an auxiliary potential V 1 (r) added to the base potential V 0 (r)) the STEP value is always halfway between the control-parameter endpoints 54,77 . Finally, since the initial force-relaxation rate (34) is an equilibrium average and the STEP value is given by a point in control-parameter space, the STEP value is simple to determine. The STEP and the strongtrap optimal protocol are the simplest protocols to determine and can be calculated analytically in many cases or numerically evaluated with relative ease.

B. Linear response
Linear-response theory can be used to determine the minimum-dissipation protocol for weak perturbations and performs relatively well at any driving speed, well beyond its strict range of validity 59,60 . For fast driving, the minimum-dissipation protocols determined from linear-response theory have jumps at the start and end of the protocol.
As mentioned in section V, the central quantity that needs to be determined to perform any optimization of the excess work is the nonequilibrium average force. In linear response, the deviation of the nonequilibrium average force from equilibrium is approximated by the integrated equilibrium force covariance This greatly simplifies the optimization procedure, since the equilibrium average depends only on the current controlparameter value and not on the entire history of control. Substituting this into the excess work (32) gives which can be optimized directly by numerical methods and can perform well at any driving speed 59,60 .

C. Slow control
The next approximation we consider is valid for slow near-equilibrium processes. This approach generalizes the paradigm of thermodynamic geometry to stochastic thermodynamics 7,118 . The generalized friction tensor endows the space of thermodynamic states with a Riemannian metric where minimum-dissipation protocols correspond to geodesics of the friction tensor. This method is widely applicable, yields a relatively simple prescription for determining minimum-dissipation protocols, and has been extended to more general settings and different forms of control. The protocols determined from this method are continuous; however, we know from exact solutions that minimumdissipation protocols can have jump discontinuities 35 , which are never optimal within the geometric framework of slow control. This arises from the slow near-equilibrium approximation used; indeed, in the limit of a slow protocol the jumps in the exact optimal protocol become negligible.
In addition to the linear-response approximation, we assume the control parameters are driven slowly compared to the system's natural relaxation timescale, so the approximation for the nonequilibrium average force (37) simplifies to Substituting into (32) yields the leading-order contribution to the excess work 118 : in terms of the generalized friction tensor with elements In analogy with fluid dynamics, this rank-two tensor is the Stokes' friction, since it produces a drag force that depends linearly on velocity. ζ j is the Hadamard product β δf j δf λ • τ j of the conjugate-force covariance (the force fluctuations) and the integral relaxation time the characteristic time for these fluctuations to die out. For overdamped dynamics, the friction can be calculated directly from the total energy as 66 where Π eq (x, λ) ≡ x −∞ dx π eq (x , λ) is the equilibrium cumulative distribution function, ∂ λj is the partial derivative with respect to λ j , and π eq (x , λ) = exp[−βV tot (x, λ)] / dx exp[−βV tot (x, λ)] is the equilibrium probability distribution.
Within the slow-protocol approximation, the excess work is minimized by a protocol with constant excess power 118 . For a single control parameter, this amounts to proceeding with velocity dλ LR /dt ∝ ζ(λ) −1/2 , which when normalized to complete the protocol in a fixed allotted time ∆t, gives dλ LR dt = ∆λ ∆t where the overline denotes the spatial average over the naive (linear) path between the control-parameter endpoints. For multidimensional control, the minimum-dissipation protocol solves the Euler-Lagrange equation where we have adopted the Einstein convention of implied summation over all repeated indices. To illustrate the steps for determining the minimum-work protocol within this approximation, Fig. 4 presents the friction matrix previously computed for the example model system of driven barrier crossing 77 . For this system the friction matrix can be directly computed from (43) and the geodesics found by numerically solving (45) with specified initial and final control parameters, as described in Refs. 71,72. Some general properties can be observed from Fig. 4: although the friction is positive semidefinite, the off-diagonal component can become negative; geodesics never have any discontinuities (following from the Riemannian geometry); geodesics tend to avoid or slow down in regions of high friction; and for driven barrier crossing the minimum-work protocol slows down and tightens the trap as it crosses the barrier. Slowing down and tightening the trap as the system crosses the barrier is consistent with the results derived within the strong-trap approximation (section IV B); however, the designed slow-control protocol lacks the jumps at the start and end of the protocol.

Extensions to the slow-control approximation
The slow-driving approximation has been generalized and extended to transitions between nonequilibrium steady states 119,120 , discrete control 68 , and stochastic control 121 .
Previously throughout section V, we assumed that the system starts in equilibrium; however, this is not always the case in applications. Periodically driven machines such as ATP synthase are often driven for a sufficiently long time as to reach a nonequilibrium steady state, breaking the initial-equilibrium and near-equilibrium assumptions. Such systems may need to transition between steady states in order to increase or decrease their output in response to variable conditions (e.g., increasing/decreasing rate of ATP production). Although determining the correct definition of dissipation is more subtle for nonequilibrium steady states (see Refs. 122-129 for detailed discussion), the slowprotocol approximation has been generalized to slow transitions between nonequilibrium steady states making use of a near-steady state approximation 119,120 . This approximation has an analogous form to the friction-tensor approximation and can be used to determine minimum-dissipation protocols for slow transitions between nonequilibrium steady states.
So far we have assumed that the protocol is continuous; however, many biological and chemical systems convert free energy stored in nonequilibrium chemical-potential differences into useful work through a series of reactions involving binding/unbinding or catalysis of small molecules. These chemical reactions typically occur on timescales much faster than the protein conformational rearrangements they couple to. Therefore, these changes are effectively instantaneous, leading to discrete control protocols. Building off of quasistatic results 130 , it has been shown that the linear-response approximation can be applied to discretely driven systems to yield an approximation analogous to the friction tensor, which can be used to determine minimum-dissipation protocols 68 . The deterministic control considered in all other sections is a reasonable assumption for single-molecule experiments; however, autonomous systems such as ATP synthase in vivo are driven not by a deterministic controller, but by other stochastic systems. Towards a description of fully autonomous machines, several advances have been made extending the friction-tensor framework to stochastic control 121 and autonomous thermodynamic cycles 131 . A detailed discussion of protocol optimization for stochastic control, as well as relations to other bounds on dissipation 132 appears in Ref. 121.

Higher-order moments and corrections
For moderate-to-fast driving, the slow-control approximation is not sufficient to accurately determine the minimumdissipation protocol. The slow-control approximation can be generalized to higher-order moments of the work distri-bution and to account for the next-order correction to the excess work. 76 This is particularly useful for nonequilibrium free-energy estimation (section VI).
In Ref. 76, the higher-order corrections are derived for the first moment second moment and n th moment We have defined the rank-three tensor and the integral n-time covariance functions The rank-three tensor [Eq. (49)] is referred to as the supra-Stokes' tensor and is indexed by superscript (2), as it corresponds to the next-order correction to dissipation beyond the Stokes' friction (46). The higher-order moments of the excess work are higher order in control-parameter velocity and are therefore smaller for slow protocols. The most notable application of this approximation is to free-energy estimation (section VI), where the supra-Stokes' tensor determines the leading-order contribution to the bias of bidirectional estimators and can be used to design protocols that minimize bias 76 .

D. Interpolated protocols
Given theory describing minimum-dissipation control in both the slow and fast limits, a simple interpolation scheme has been developed to design protocols that reduce dissipation at all driving speeds 41,77 . The interpolated protocol has an initial jump (λ STEP − λ i )/(1 + ∆t/τ ) and a final jump (λ f − λ STEP )/(1 + ∆t/τ ), and follows the slow-control path between them, with τ the crossover duration and λ slow (t) the solution to (45). This guarantees that the protocol approaches the minimum-dissipation protocol in both the fast and slow limits.

VI. FREE-ENERGY ESTIMATION
Computational and experimental measurements of free-energy differences are essential to the determination of stable equilibrium phases of matter, relative reaction rates, and binding affinities of chemical species, 133 and the identification and design of novel protein-binding ligands for drug discovery. [134][135][136][137] Current methods for determining free-energy differences rely on costly experimentation, which can be reduced through screening with efficient computational techniques. [133][134][135][136][137][138][139][140] It has been shown that the precision and accuracy of standard free-energy estimators are reduced when estimated from a protocol inducing large dissipation, 76,141,142 and that thermodynamic geometry can be applied to improve free-energy estimates. 76,[142][143][144][145][146] Free-energy differences are often estimated by measuring the work incurred during a parametric control protocol that drives the system between control-parameter endpoints corresponding to target states. Unidirectional estimators determine the free-energy difference from the work done by a forward protocol driving from initial to final controlparameter values, while bidirectional estimators additionally use reverse protocols that drive from final to initial control-parameter values. The simplest estimator is the mean-work estimator, 141 which estimates the free-energy difference by the average work and for any non-quasistatic (finite-speed) protocol yields a biased estimate. For an unbiased estimate, the Jarzynski estimator (derived from the Jarzynski equality 147 ) estimates the free-energy difference from the exponentially averaged work. The mean-work and Jarzynski estimators can be used as either uni-or bidirectional estimators; however, if bidirectional data is available the maximum log-likelihood estimator is Bennett's acceptance ratio. 148 For a large number of samples, Bennett's acceptance ratio yields the minimum variance of any unbiased estimator. 149,150 Near equilibrium, Taylor expanding all of these estimators for small excess work (small dissipation) gives the meanwork estimator. 76 For an equal number of samples in the forward and reverse directions, near equilibrium the variance of any of these estimators is 76 where the second line follows from the fluctuation-dissipation relation for the work distribution. 76,141 The bias is given by the asymmetry between forward Λ and reverse Λ † dissipation, generated from skewed work fluctuations: Equations (52a) and (53) not only hold near equilibrium but also when only a single sample is taken in the forward and reverse directions, since then the average of a function is equal to the function of the average (e.g., e x = e x ). For a slow protocol the variance is approximated by Thus the protocol designed to reduce the variance follows geodesics of ζ j , and for one-dimensional control proceeds at velocityλ ∝ (ζ) −1/2 . This, or the related force-variance metric 142 , has been used to improve the precision of calculated binding potentials of mean force [143][144][145][146] .
Unlike the variance, the protocol that maximizes the accuracy (minimizes bias) is different for unidirectional and bidirectional estimators. For unidirectional Jarzynski and mean-work estimators, near equilibrium the minimum-bias protocol is simply the minimum-dissipation protocol (protocol that minimizes (46)) and therefore to leading order is optimized by the same protocol that minimizes (54). For a slow protocol, the bias from Bennett's acceptance ratio can be approximated as where the second line follows from (46). The protocol designed to reduce the (magnitude of) bias thus follows geodesics of the cubic Finsler metric ζ (2) j m , simplifying for one-dimensional control toλ ∝ ζ (2) −1/3 .

VII. COMPARISON BETWEEN CONTROL STRATEGIES
In this section we compare naive and designed protocols based on the methods described in the previous sections: interpolated protocols combining STEP and slow-protocol approximations (section V D), strong-trap approximation (section IV B), and full control (section IV A). We again turn to the example model system of driven barrier crossing to assess similarities, differences, and performance of the designed protocols. The naive protocol serves as a baseline which the designed protocols should outperform, and full control as a bound on what parametric control could possibly achieve. and slows down and tightens the trap as it crosses the barrier. The behavior of the designed protocols can be understood in terms of the full-control solution. In one dimension the minimum-dissipation protocol linearly drives the quantiles of the probability distribution between the initial and final distributions. In the naive protocol, since it has constant stiffness, the probability distribution spreads out as it crosses the barrier due to the negative curvature of the energy landscape at the barrier. To compensate for this, the designed protocols tighten as they cross the barrier; additionally, to compensate for the changes in the gradient of the energy landscape, the designed protocols slow down as they cross the barrier. Fig. 6 shows the excess work of the designed protocols compared to the naive protocol. For long duration (slow protocol) all of the designed protocols significantly outperform the naive, with the difference between the minimum dissipation possible from full control indistinguishable from the dissipation from the interpolated protocol in this limit. While the approximations made in the interpolated protocol become exact in the long-duration limit, the same is not true for the strong-trap approximation. As a result, the strong-control protocol has slightly higher dissipation in the long-duration limit, but would achieve the minimal dissipation in the limit of high trap stiffness. Furthermore, for intermediate protocol duration (∆t ∼ τ D ), the strong control performs the best of the approximations since the approximation does not explicitly depend on the protocol duration. For short duration (fast protocols), all the designed protocols perform similarly well.
In summary, the designed protocols perform similarly well, so it seems reasonable to choose the control strategy which is simplest to determine, provided the approximations/assumptions required are satisfied. Since the strong-trap approximation has explicit solutions for the designed protocol, it will generally be the simplest to determine; however, it is not as widely applicable as the interpolated protocols.

VIII. PERSPECTIVE AND OUTLOOK
We have reviewed optimal control in stochastic thermodynamics, from full control (section IV) to parametric control (section V). General methods are known for determining minimum-dissipation protocols for parametric control ranging from strong (section IV B) to weak (section V B) and fast (section V A) to slow (section V C). These approximations  fill out the four limits of parametric control (Fig. 2) and can be combined to design protocols that reduce dissipation at any driving speed (section V D). These designed protocols reproduce key features determined by exact solutions for Gaussian distributions and quadratic trapping potentials, such as control-parameter jumps at the start and end of the protocol (section V A) and linear driving of the distribution quantiles between specified endpoints (section IV A). For the model system of driven barrier crossing (section II), we compare the interpolated, strong-control, and fullcontrol solutions (section VII). The designed protocols significantly outperform the naive (linear) protocol. Strong control has explicit solutions for the minimum-dissipation protocol, making it the simplest approximation to use; however, it is only applicable to overdamped dynamics with strong trapping potentials. Figure 2 shows the limits in which we have solutions for the minimum-dissipation protocol for parametric control. The linear-response and slow-driving approximations have been applied to several different types of systems and control (section V C). A promising area of future study would be to explore if extensions and generalizations can be made to strong and fast control. Indeed, it has recently been shown that the fast-protocol approximation can be extended to quantum systems, classical Hamiltonian dynamics, and optimization of the variance of the work distribution. 151 Another extension to consider is to allow for position-dependent diffusivity, which generically arises when a highdimensional system is represented in a lower-dimensional space. 152 For example, DNA hairpin experiments and molecular dynamics simulations take place in three spatial dimensions but are often represented as one-dimensional diffusive processes. 18,153,154 This requires averaging out the behavior in the other two dimensions, and any inhomogeneity in these dimensions will result in position-dependent effective diffusivity in the one-dimensional representation. Measured diffusivities in molecular dynamics simulations often vary with position 153 , and although accurate measurements remain a technical challenge, some hairpin experiments report a position-dependent diffusivity. 155 Position-dependent diffusivity can alter the kinetic transition state of protein folding 156 , which will impact the design of minimumdissipation protocols. Therefore, position-dependent diffusivity may be an important consideration in some applications. The minimum-dissipation protocol under full control in inhomogeneous environments (e.g., position-dependent diffusivity) has been treated in detail in Ref. 157.
An important aspect of designing protocols to minimize dissipation in both experiments and theory is the choice of control parameters and number of control parameters. For the model system of driven barrier crossing, it has been shown that designed protocols with control over both trap center and stiffness (two-parameter control) significantly reduces dissipation compared to designed protocols that can only adjust the trap center (single-parameter control). 77 However, adding even more control parameters does not significantly reduce dissipation any further, because this system is well approximated as a one-dimensional Gaussian, for which the full-control solution only requires two parameters to adjust the mean and variance (section IV A). Although this phenomenon is well explained for onedimensional overdamped systems, considerably less is known more generally. For example, Ref. 72 compared one-, two-, and four-parameter control of the Ising model and found significant qualitative and quantitative differences between the designed protocols. Since the full-control solution for this system is not yet known, we cannot explain this phenomenon the same way we did for driven barrier crossings. Beyond simply the number of control parameters, it remains an open question as to which control parameters are the most important when designing protocols to minimize dissipation. The choice of control parameters and number of control parameters will be important for experimental and computational applications, such as free-energy estimation.
Leveraging optimal-transport theory appears to be a promising approach towards a deeper understanding of optimal control in stochastic thermodynamics. Full-control solutions based on optimal-transport theory were initially only applicable to continuous classical systems with overdamped dynamics 44,46 ; however, recent studies have begun to expand their range of applicability to discrete-state and quantum systems 47,[52][53][54][55][56] .
The full-control solutions give an idealized view of optimal control and will always outperform parametric control, but can be used as the ideal solution that parametric control can strive towards and draw insight from. For example, linearly driving the quantiles between the initial and final endpoints is the minimum-dissipation protocol for a onedimensional system under full control; this is reasonably well reproduced by parametric control of driven barrier crossing (section VII). Furthermore, it has recently been shown that optimal-transport theory can be leveraged to design minimum-dissipation protocols under parametric control 54 , which is a promising new technique for determining exact minimum-dissipation protocols at any driving speed or strength of driving.
In this review we focused on protocols which reduce the average work or entropy production; however, higher-order moments (e.g., variance or skewness), or individual trajectory optimization 158 may also be of interest for strongly fluctuating systems. Ref. 159 compared for a breathing harmonic trap the protocols that minimize work variance with those that minimize average work, finding that minimum-average-work and minimum-work-variance protocols qualitatively differ for intermediate-to-fast driving speeds. In contrast, for slow near-equilibrium control, the same protocols minimize average work and work variance (section V C), so a more complete description of minimum-variance protocols far from equilibrium is desirable.
This review has focused on systems driven by external control parameters, which is adequate for describing the experimental systems discussed in section II; however, this does not accurately describe autonomous machines. For example, ATP synthase in vivo is driven by a proton gradient across the mitochondrial membrane which drives the coupled F o and F 1 components. In this system, coupling between the components means that none of the components can be treated as external, and thus it is not obvious how the present discussion of optimal control applies to such autonomous machines. Some first steps towards the description of autonomous molecular machines is discussed in section V C; however, more work is required towards a full description of efficient autonomous machines. 160 The majority of the studies on optimal control in stochastic thermodynamics have focused on theoretical understanding and simple toy models. As we have shown in this review, there is a deep understanding of minimum-dissipation protocols at both slow and fast driving speed and both weak and strong driving strength (Fig. 2). It would be interesting to see how these results apply to real physical systems and if they are able to achieve the promising results predicted by theory. The two most straightforward applications are to relatively simple experimental model systems (section II) in an analogous fashion to Ref. 75, and to free-energy estimation as discussed in section VI, in a similar fashion to Refs. 143-146.