Minimal Action Method using Optimal Control

The Minimal Action Method (MAM) 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, using the classical Maier-Stein model as a benchmark example.

Required Packages

using OptimalControl
using NLPModelsIpopt
using Plots, Printf

Problem Statement

We aim to find the most probable transition path between two stable states of a stochastic dynamical system. For a system with deterministic dynamics $f(x)$ and small noise, the transition path minimizes the action functional:

\[S[x(\cdot), u(\cdot)] = \int_0^T \|u(t) - f(x(t))\|^2 \, dt\]

subject to the path dynamics:

\[\dot{x}(t) = u(t), \quad x(0) = x_0, \quad x(T) = x_f\]

where $x_0$ and $x_f$ are the initial and final states, and $T$ is the transition time.

Problem Setup

We 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 with the @init macro:

# Time horizon
T = 50

# Helper functions for initial state guess
L(t) = -(1 - t/T) + t/T      # Linear interpolation from -1 to 1
P(t) = 0.3*(-L(t)^2 + 1)     # Parabolic arc (approximates saddle crossing)

init = @init ocp(T) begin
    # Linear interpolation for x₁
    x₁(t) := L(t)
    # Parabolic guess for x₂
    x₂(t) := P(t)
    # Control from vector field
    u(t) := f(L(t), P(t))
end

Solving the Problem

We solve the problem in two steps for better accuracy:

# First solve with coarse grid
sol = solve(ocp(T); init=init, grid_size=50)

# Refine solution with finer grid
sol = solve(ocp(T); init=sol, grid_size=1000)

# Objective value
objective(sol)

0.24942662080398548

Visualizing Results

Let's plot the solution trajectory and phase space:

plot(sol)

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

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.

# Continuation function to avoid global variables
function continuation_mam(Ts; init_guess=init)
    objectives = Float64[]
    iterations_list = Int[]
    current_sol = init_guess

    opts = (display=false, grid_size=1000, tol=1e-8)
    for T in Ts
        current_sol = solve(ocp(T); init=current_sol, opts...)
        push!(objectives, objective(current_sol))
        push!(iterations_list, iterations(current_sol))
    end

    return objectives, iterations_list, current_sol
end

# Perform continuation
Ts = range(50, 100, 101)
objectives, iters, final_sol = continuation_mam(Ts)

# Display results
println(" Time   Objective     Iterations")
for i in eachindex(Ts)
    @printf("%6.2f  %9.6e  %d\n", Ts[i], objectives[i], iters[i])
end

 Time   Objective     Iterations
 50.00  2.494266e-01  434
 50.50  2.494244e-01  46
 51.00  2.494223e-01  55
 51.50  2.494203e-01  53
 52.00  2.494184e-01  49
 52.50  2.494169e-01  31
 53.00  2.494150e-01  74
 53.50  2.494134e-01  61
 54.00  2.494119e-01  62
 54.50  2.494104e-01  61
 55.00  2.494091e-01  61
 55.50  2.494078e-01  58
 56.00  2.494066e-01  64
 56.50  2.494055e-01  67
 57.00  2.494044e-01  34
 57.50  2.494034e-01  58
 58.00  2.494025e-01  52
 58.50  2.494016e-01  66
 59.00  2.494008e-01  52
 59.50  2.494000e-01  63
 60.00  2.493993e-01  54
 60.50  2.493992e-01  15
 61.00  2.493981e-01  92
 61.50  2.493975e-01  53
 62.00  2.493970e-01  36
 62.50  2.493966e-01  65
 63.00  2.493962e-01  71
 63.50  2.493958e-01  74
 64.00  2.493955e-01  70
 64.50  2.493952e-01  44
 65.00  2.493949e-01  62
 65.50  2.493947e-01  52
 66.00  2.493945e-01  57
 66.50  2.493944e-01  64
 67.00  2.493943e-01  64
 67.50  2.493942e-01  75
 68.00  2.493942e-01  63
 68.50  2.493942e-01  50
 69.00  2.493942e-01  67
 69.50  2.493942e-01  71
 70.00  2.493943e-01  86
 70.50  2.493944e-01  75
 71.00  2.493945e-01  90
 71.50  2.493947e-01  86
 72.00  2.493949e-01  44
 72.50  2.493951e-01  71
 73.00  2.493954e-01  92
 73.50  2.493956e-01  53
 74.00  2.493959e-01  65
 74.50  2.493962e-01  58
 75.00  2.493966e-01  43
 75.50  2.493969e-01  64
 76.00  2.493973e-01  70
 76.50  2.493977e-01  90
 77.00  2.493981e-01  63
 77.50  2.493985e-01  90
 78.00  2.493990e-01  57
 78.50  2.493995e-01  81
 79.00  2.494000e-01  53
 79.50  2.494005e-01  67
 80.00  2.494010e-01  77
 80.50  2.494016e-01  65
 81.00  2.494022e-01  93
 81.50  2.494031e-01  41
 82.00  2.494034e-01  111
 82.50  2.494040e-01  47
 83.00  2.494046e-01  14
 83.50  2.494053e-01  76
 84.00  2.494060e-01  48
 84.50  2.494067e-01  83
 85.00  2.494074e-01  79
 85.50  2.494081e-01  44
 86.00  2.494088e-01  80
 86.50  2.494096e-01  102
 87.00  2.494104e-01  80
 87.50  2.494111e-01  100
 88.00  2.494119e-01  91
 88.50  2.494127e-01  60
 89.00  2.494136e-01  53
 89.50  2.494144e-01  102
 90.00  2.494152e-01  51
 90.50  2.494161e-01  105
 91.00  2.494170e-01  105
 91.50  2.494179e-01  117
 92.00  2.494188e-01  121
 92.50  2.494197e-01  110
 93.00  2.494206e-01  108
 93.50  2.494216e-01  117
 94.00  2.494225e-01  48
 94.50  2.494235e-01  118
 95.00  2.494245e-01  107
 95.50  2.494254e-01  22
 96.00  2.494264e-01  111
 96.50  2.494275e-01  122
 97.00  2.494286e-01  33
 97.50  2.494296e-01  51
 98.00  2.494306e-01  183
 98.50  2.494316e-01  98
 99.00  2.494327e-01  114
 99.50  2.494338e-01  91
100.00  2.494348e-01  127

We can now analyze the results and identify the optimal transition time:

# Find optimal time
idx_min = argmin(objectives)
T_min = Ts[idx_min]
obj_min = objectives[idx_min]

@printf("Optimal transition time: T* = %.2f\n", T_min)
@printf("Minimal action: S* = %.6e\n", obj_min)

Optimal transition time: T* = 68.50
Minimal action: S* = 2.493942e-01

Let us visualize the evolution of the objective function with respect to the transition time:

plt1 = scatter(Ts, log10.(objectives), xlabel="Time", label="Objective (log10)")
vline!(plt1, [T_min], label="Minimum at T=$(round(T_min, digits=1))", z_order=:back)
plot(plt1, size=(800,400))