Minimal Action Method using Optimal Control

The Minimal Action Method is a numerical technique for finding the most probable transition pathway between stable states in stochastic dynamical systems. It achieves this by minimizing an action functional that represents the path's deviation from the deterministic dynamics, effectively identifying the path of least resistance through the system's landscape. This tutorial demonstrates how to implement MAM as an optimal control problem.

Required Packages

using OptimalControl
using NLPModelsIpopt
using Plots, Printf

Problem Setup

We'll consider a 2D system with a double-well flow, called the Maier-Stein model. It is a famous benchmark problem as it exhibits non-gradient dynamics with two stable equilibrium points at (-1,0) and (1,0), connected by a non-trivial transition path. The system's deterministic dynamics are given by:

Define the vector field

f(u, v) = [u - u^3 - 10*u*v^2, -(1 - u^2)*v]
f(x) = f(x...)

Optimal Control Formulation

The minimal action path minimizes the deviation from the deterministic dynamics:

function ocp(T)
    action = @def begin
        t ∈ [0, T], time
        x ∈ R², state
        u ∈ R², control
        x(0) == [-1, 0]    # Starting point (left well)
        x(T) == [1, 0]     # End point (right well)
        ẋ(t) == u(t)       # Path dynamics
        ∫(sum((u(t) - f(x(t))) .^ 2)) → min  # Minimize deviation from deterministic flow
    end
    return action
end

Initial Guess

We provide an initial guess for the path using a simple interpolation:

T = 50 # Time horizon

x1(t) = -(1 - t/T) + t/T # Linear interpolation for x₁
x2(t) = 0.3(-x1(t)^2 + 1) # Parabolic guess for x₂

x(t) = [x1(t), x2(t)]
u(t) = f(x(t))

init = (state=x, control=u) # Initial guess

Solving the Problem

We solve the problem in two steps for better accuracy:

sol = solve(ocp(T); init=init, grid_size=50) # First solve with coarse grid
sol = solve(ocp(T); init=sol, grid_size=1000) # Refine solution
objective(sol) # Objective value

0.249426627800687

Visualizing Results

Let's plot the solution trajectory and phase space:

plot(sol)

MLP = state(sol).(time_grid(sol))
scatter(
    first.(MLP),
    last.(MLP);
    title="Minimal Action Path",
    xlabel="u",
    ylabel="v",
    label="Transition path",
) # Phase space plot

The resulting path shows the most likely transition between the two stable states given a transient time $T=50$, minimizing the action functional while respecting the system's dynamics.

Minimize with respect to T

To find the maximum likelihood path, we also need to minimize the transient time T. Hence, we perform a discrete continuation over the parameter T by solving the optimal control problem over a continuous range of final times T, using each solution to initialize the next problem.

objectives = []
Ts = range(1, 100, 30)
sol = solve(ocp(Ts[1]); display=false, init=init, grid_size=50)
println(" Time   Objective     Iterations")
for T in Ts
    global sol = solve(ocp(T); display=false, init=sol, grid_size=1000, tol=1e-8)
    @printf("%6.2f  %9.6e  %d\n", T, objective(sol), iterations(sol))
    push!(objectives, objective(sol))
end

 Time   Objective     Iterations
  1.00  4.076020e+00  2
  4.41  5.353696e-01  75
  7.83  3.026580e-01  17
 11.24  2.632190e-01  16
 14.66  2.543237e-01  13
 18.07  2.516442e-01  22
 21.48  2.505987e-01  28
 24.90  2.501063e-01  34
 28.31  2.498425e-01  37
 31.72  2.496876e-01  62
 35.14  2.495905e-01  51
 38.55  2.495268e-01  43
 41.97  2.494835e-01  45
 45.38  2.494535e-01  58
 48.79  2.494325e-01  84
 52.21  2.494177e-01  150
 55.62  2.494075e-01  61
 59.03  2.494008e-01  129
 62.45  2.493966e-01  138
 65.86  2.493964e-01  33
 69.28  2.493942e-01  225
 72.69  2.493952e-01  139
 76.10  2.493974e-01  128
 79.52  2.494005e-01  90
 82.93  2.494045e-01  100
 86.34  2.494094e-01  135
 89.76  2.494148e-01  168
 93.17  2.494209e-01  138
 96.59  2.494276e-01  122
100.00  2.494348e-01  139

T_min = Ts[argmin(objectives)]
plt1 = scatter(Ts, log10.(objectives); xlabel="Time", label="Objective (log10)")
vline!(plt1, [T_min]; label="Minimum", z_order=:back)
plt2 = scatter(
    Ts[10:30], log10.(objectives[10:30]); xlabel="Time", label="Objective (log10)"
)
vline!(plt2, [T_min]; label="Minimum", z_order=:back)
plot(plt1, plt2; layout=(2, 1), size=(800, 800))

This page was generated using Literate.jl.