Optimal Control of Engineered Swift Equilibration of Nanomechanical Oscillators

TL;DR

The paper tackles how to optimally steer stochastic nanomechanical oscillators between states by distinguishing transitions that minimize dissipation from those that minimize work, within a Bolza-form optimal-control framework. It models the underdamped, Gaussian dynamics of a one-dimensional oscillator with time-dependent stiffness and force parameters, and treats the control inputs with convex or barrier penalties to ensure well-posedness and numerical tractability. Central results show that, for Gaussian states, optimal protocols exhibit turnpike behaviour: in the bulk of the horizon the dynamics align to a universal centre manifold determined solely by the running cost, with boundary layers near the horizon ends; this explains why boundary jumps are neither required nor physical. Numerically, engineered swift equilibration protocols outperform minimum-entropy-production strategies in many scenarios, and the centre-manifold analysis provides a robust, scalable way to predict the qualitative behaviour of optimal controls across penalty choices and regimes.

Abstract

We propose a reformulation of the problem of optimally controlled transitions in stochastic thermodynamics. We impose that any terminal cost specified by a thermodynamic functional should depend only on state variables and not on control protocols, according to the canonical Bolza form. In this way, we can unambiguously discriminate between transitions at minimum dissipation between genuine equilibrium states, and transitions at minimum work driving a system from a genuine equilibrium to a non-equilibrium state. For underdamped dynamics subject to a mechanical force, genuine equilibrium means a Maxwell-Boltzmann probability distribution defining a vanishing current velocity. Transitions at minimum dissipation between equilibria are a model of optimal swift engineered equilibration. Remarkably, we show that transitions at minimum work do not directly imply explicit boundary conditions on terminal values of parameters of the mechanical force and on control protocols. Thus, the problem often discussed in the literature, that optimal protocols need terminal jumps to satisfy boundary conditions, completely disappears. The quantitative properties of optimal controls are entirely determined by the form of the penalty modelling an experimental setup. More generally, we use centre manifold theory to analytically account for the tendency of optimal controls to exhibit a turnpike property: optimal protocols in the bulk of the control horizon tend to converge to a universal centre manifold determined only by the running cost. Exponential deviations from the centre manifold occur at the ends of the control horizon in order to satisfy the boundary conditions. Our findings are supported numerically.

Figures (10)

• Figure 1: Engineered swift equilibration minimising the entropy production \ref{['ep:cost']}, subject to the harmonic penalty \ref{['penalty:harmonic']}. The solution is computed using a direct optimisation on the cost functional \ref{['ep:cost']} (blue) with InfiniteOpt.jl PuZhHoZa2022 and solved using Interior Point Optimisation IPOpt WacA2009 and using $g=0.01$. The central manifold solution is shown for decreasing values of $g$: $g=0.01$ orange, $g=0.001$ green, $g=0.0001$ maroon; and is found by numerically integrating the set of differential equations \ref{['sf:nonuniv']}, \ref{['sf:univ']} with a collocation method from the package DifferentialEquations.jl rackauckas2017differentialequationsLobatto_Jay2015. We fix $\Lambda=\sqrt{2},\,\varepsilon=1,\,t_{\mathfrak{f}}=3$ and boundary conditions are described in \ref{['ep_num:bc']}, with $\sigma^2_0 = 1$ and $\sigma^2_{t_{\mathfrak{f}}}=2$. Results produced using the code available at github_link
• Figure 2: Engineered swift equilibration minimising the entropy production \ref{['ep:cost']}. We find the solution by a direct optimisation of the cost functional with a hard penalty \ref{['penalty:hard']} (blue) and the solution of system of differential equations specifying the first order optimality conditions \ref{['ext:eqs']} with a logarithmic penalty \ref{['penalty:log']} (orange-dashed). We fix $\Lambda=9,\,\varepsilon=1,\,t_{\mathfrak{f}}=3$ and use $g=0.001$ and $g=0$ for the logarithmic and hard penalties respectively. The boundary conditions are given by \ref{['noneq_bc']}, with $\sigma^2_0 = 1$ and $\sigma^2_{t_{\mathfrak{f}}}=2$. Numerical methods are as described in Fig. \ref{['fig:harmonic']}.
• Figure 3: Engineered swift equilibration minimising the entropy production \ref{['ep:cost']} for a contraction. We use a harmonic penalty \ref{['penalty:harmonic']} and find the solution with a direct optimisation of the cost functional (blue lines) and as a solution of the first order conditions (orange, dashed). We fix $\Lambda=\sqrt{2},\,\varepsilon=1,\,t_{\mathfrak{f}}=3$ and $g=0.01$. Boundary conditions are given by \ref{['noneq_bc']}, with $\sigma^2_0 = 1$ and $\sigma^2_{t_{\mathfrak{f}}}=1/2$. Numerical methods are as described in Fig. \ref{['fig:harmonic']}.
• Figure 4: Entropy production as a function of the time horizon $t_f$ where the stiffness is a state \ref{['S1']} (blue, triangle) and where the stiffness is a control \ref{['S2']} (orange, circle). When the stiffness is a state, we use a hard penalty $-\Lambda\leq\lambda_t\leq\Lambda$. Panels (a) and (b) constrain the stiffness $\mathscr{k}_{t}$ in the interval $0.2\leq\mathscr{k}_{t}\leq 1.2$. Panels (c) and (d) have the constraint $-0.1\leq\mathscr{k}_{t}\leq 1.2$, allowing for negative values of $\mathscr{k}_{t}$. To model \ref{['S2']}, only $\mathscr{k}_{t}$ is constrained and there is no constraint on $\lambda_t$ (its time derivative). Results are computed by a direct optimisation of the cost functional. For all plots we use $g=0$ and $\varepsilon=1$; in (a) and (c) we set the size of the hard penalty $\Lambda =1$; for (b) and (d) we set $\Lambda =10$. The squared Wasserstein-2 distance $\mathcal{W}_2$ between the initial and final position marginals (blue, dashed) is estimated numerically using $n=20\,000$ independent samples using simulated evolutions of the dynamics controlled by $\mathscr{k}_t$ The blue shaded region indicates the standard error of the estimator. The numerical methods used are as those in Fig. \ref{['fig:harmonic']}.
• Figure 5: Comparison of engineered swift equilibration at minimum entropy production when the stiffness is the control \ref{['S2']} (orange) and when the stiffness is a state \ref{['S1']} (blue, green). The stiffness $\mathscr{k}_{t}$ is constrained in the interval $0.2\leq\mathscr{k}_{t}\leq 1.2$ in both cases. We use a hard penalty \ref{['penalty:hard']} to model the case when the stiffness is a state, with $\Lambda=1$ (blue) and $\Lambda=10$ (green). Results are computed by a direct optimisation (see Fig. \ref{['fig:harmonic']}) with parameters $\varepsilon=1$ and $g=0$. We use boundary conditions \ref{['ep_num:bc']} when the stiffness is a state (SI) and remove the boundary conditions on $\mathscr{k}_{t}$ for stiffness as the control (case \ref{['S2']}). We use $\sigma^2_0 = 1$ and $\sigma^2_{t_{\mathfrak{f}}}=2$. Final quantities are smoothed by a convolution with box filter over an interval of approximately $0.067$