Mocca.jl

Appearance

Optimization of a cyclic process ​

We optimise the four-stage vacuum swing adsorption (VSA) cycle of Haghpanah et al. (2013) at cyclic steady state, twice over:

  1. maximise the CO₂ recovery, keeping the purity of the CO₂ product at 90% or above;

  2. from that design, minimise the energy per tonne of CO₂ captured, keeping the purity at 90% and the recovery at 60% or above.

Each evaluation finds the cyclic steady state, starting from the last one, and its gradient with respect to the design variables from adjoints (see the gradients example). A gradient-based optimiser (BFGS within the bounds) then needs a few tens of evaluations, where the genetic algorithm of the paper evaluates thousands. For more details about the VSA modelling, see the process from its parts example.

julia
import Mocca
using Mocca.Units
using CairoMakie

The process and its design variables ​

Four of the decision variables of section 7.1.1 of the paper: the feed velocity, the evacuation pressure P_L, the blowdown pressure P_I and the adsorption time. The bounds of the pressures and the velocity are those of Table 5 of the paper. The paper also varies the blowdown and evacuation times, which we leave out to keep the example short.

Since a duration is among the variables, a stage of a different duration keeps the time steps it has now, stretched in proportion, so that the performance changes smoothly with the duration.

julia
process = Mocca.four_stage_vsa_process()
Mocca.add_variable!(process, :v_feed, ("adsorption", :velocity); bounds = (0.1, 2.0), description = "interstitial feed velocity [m/s]")
Mocca.add_variable!(process, :P_L, ("evacuation", :to); bounds = (0.05bar, 0.5bar), description = "evacuation pressure")
Mocca.add_variable!(process, :P_I, ("blowdown", :to); bounds = (0.05bar, 0.5bar), description = "blowdown pressure")
Mocca.add_variable!(process, :t_ads, ("adsorption", :duration); bounds = (10.0, 100.0), description = "adsorption time [s]")
Process four-stage VSA: CO2/N2 feed at 1.0 bar, 298.15 K, capturing CO2
  Bed: fixed bed of zeolite 13X, 30 cells
    bottom V_feed (feed), from boundary at 1.0 bar; V_vacuum (extract), to boundary at 0.1 bar
    top    V_product (product), to boundary at 1.0 bar
  cycle of 4 stages, 100.0 s:
    pressurisation                        15.0 s   V_feed
    adsorption                            15.0 s   V_feed, V_product
    blowdown                              30.0 s   V_product
    evacuation                            40.0 s   V_vacuum
  variables: v_feed, P_L, P_I, t_ads

The blowdown must end above the evacuation pressure, which the optimiser takes as a linear constraint on the variables, P_L − P_I ≤ −0.01 bar:

julia
blowdown_above_evacuation = [(Dict(:P_L => 1.0, :P_I => -1.0), -0.01bar)]
1-element Vector{Tuple{Dict{Symbol, Float64}, Float64}}:
 (Dict(:P_I => -1.0, :P_L => 1.0), -1000.0)

The design we start from, and a line for each design the optimiser tries:

julia
show_design(values) = join(["$k = " * (k in (:P_L, :P_I) ? "$(round(v / bar, digits = 3)) bar" : "$(round(v, sigdigits = 3))")
    for (k, v) in pairs(values)], ", ")
show_performance(m) = "purity $(round(100m.purity, digits = 1))%, recovery $(round(100m.recovery, digits = 1))%, " *
    "$(round(m.specific_energy, digits = 1)) kWh/t"
report(values, r) = println("  ", show_design(values), ": ", show_performance(r.metrics))

start = Mocca.performance(process)
println("start: ", show_performance(start))
start: purity 79.5%, recovery 36.7%, 141.3 kWh/t

1. Maximise the recovery ​

The purity target is a penalty, so the optimum may miss it slightly.

julia
max_recovery = Mocca.optimise!(process; maximise = :recovery, at_least = (purity = 0.9,),
    linear_constraints = blowdown_above_evacuation, max_it = 12, log = report)
println("$(max_recovery.evaluations) evaluations: ", show_design(max_recovery.values))
println(show_performance(max_recovery.metrics))
  v_feed = 1.0, P_L = 0.1 bar, P_I = 0.2 bar, t_ads = 15.0: purity 79.5%, recovery 36.7%, 141.3 kWh/t
  v_feed = 1.03, P_L = 0.09 bar, P_I = 0.183 bar, t_ads = 19.5: purity 85.6%, recovery 39.2%, 130.2 kWh/t
  v_feed = 1.04, P_L = 0.087 bar, P_I = 0.172 bar, t_ads = 21.3: purity 88.0%, recovery 40.0%, 128.1 kWh/t
  v_feed = 1.05, P_L = 0.085 bar, P_I = 0.164 bar, t_ads = 22.4: purity 89.5%, recovery 40.5%, 127.5 kWh/t
  v_feed = 1.05, P_L = 0.083 bar, P_I = 0.162 bar, t_ads = 22.6: purity 89.9%, recovery 41.5%, 127.2 kWh/t
  v_feed = 1.06, P_L = 0.08 bar, P_I = 0.161 bar, t_ads = 22.7: purity 90.1%, recovery 43.6%, 126.6 kWh/t
  v_feed = 1.06, P_L = 0.074 bar, P_I = 0.159 bar, t_ads = 22.9: purity 90.3%, recovery 47.8%, 126.5 kWh/t
  v_feed = 1.06, P_L = 0.055 bar, P_I = 0.153 bar, t_ads = 23.5: purity 90.6%, recovery 60.2%, 132.7 kWh/t
  v_feed = 1.05, P_L = 0.05 bar, P_I = 0.151 bar, t_ads = 23.2: purity 90.2%, recovery 64.0%, 137.4 kWh/t
LBFGS: Line search at max step size, Wolfe conditions not satisfied for this step
  v_feed = 0.998, P_L = 0.05 bar, P_I = 0.149 bar, t_ads = 22.1: purity 89.3%, recovery 64.9%, 140.9 kWh/t
  v_feed = 1.04, P_L = 0.05 bar, P_I = 0.151 bar, t_ads = 23.0: purity 90.0%, recovery 64.2%, 138.0 kWh/t
  v_feed = 1.02, P_L = 0.05 bar, P_I = 0.15 bar, t_ads = 22.6: purity 89.7%, recovery 64.5%, 139.4 kWh/t
  v_feed = 1.04, P_L = 0.05 bar, P_I = 0.15 bar, t_ads = 22.9: purity 90.0%, recovery 64.2%, 138.3 kWh/t
  v_feed = 1.03, P_L = 0.05 bar, P_I = 0.15 bar, t_ads = 22.8: purity 89.9%, recovery 64.3%, 138.8 kWh/t
  v_feed = 1.0, P_L = 0.05 bar, P_I = 0.147 bar, t_ads = 22.7: purity 89.9%, recovery 64.5%, 140.1 kWh/t
  v_feed = 0.957, P_L = 0.05 bar, P_I = 0.14 bar, t_ads = 22.4: purity 89.9%, recovery 64.7%, 143.0 kWh/t
  v_feed = 0.509, P_L = 0.05 bar, P_I = 0.082 bar, t_ads = 19.4: purity 92.2%, recovery 61.8%, 203.7 kWh/t
  v_feed = 0.778, P_L = 0.05 bar, P_I = 0.117 bar, t_ads = 21.2: purity 90.2%, recovery 65.2%, 157.7 kWh/t
  v_feed = 0.799, P_L = 0.05 bar, P_I = 0.121 bar, t_ads = 21.5: purity 90.1%, recovery 65.2%, 155.1 kWh/t
  v_feed = 0.769, P_L = 0.05 bar, P_I = 0.125 bar, t_ads = 22.8: purity 89.8%, recovery 65.6%, 153.2 kWh/t
  v_feed = 0.782, P_L = 0.05 bar, P_I = 0.123 bar, t_ads = 22.2: purity 89.9%, recovery 65.5%, 154.0 kWh/t
21 evaluations: v_feed = 0.782, P_L = 0.05 bar, P_I = 0.123 bar, t_ads = 22.2
purity 89.9%, recovery 65.5%, 154.0 kWh/t

The evacuation pressure reaches its lower bound within ten evaluations: a deeper vacuum takes more CO₂ off the adsorbent, so the recovery rises as P_L falls. Haghpanah et al. find the same (their Figure 11), and reach 90% purity and recovery together only with P_L ≤ 0.02 bar. The blowdown pressure then sets the purity: blowing down further removes more of the nitrogen, but loses more CO₂ with it.

2. Minimise the energy ​

The optimiser starts from the design above, where the process was left. Each target is a penalty of weight (1000 by default) times the square of the shortfall, so the objective must be of order one: with the energy in units of 100 kWh/t, a shortfall of 1 percentage point of recovery costs as much as 10 kWh/t.

julia
min_energy = Mocca.optimise!(process; minimise = :specific_energy, at_least = (purity = 0.9, recovery = 0.6),
    scale = (specific_energy = 100.0,), linear_constraints = blowdown_above_evacuation, max_it = 15, log = report)
println("$(min_energy.evaluations) evaluations: ", show_design(min_energy.values))
println(show_performance(min_energy.metrics))
  v_feed = 0.782, P_L = 0.05 bar, P_I = 0.123 bar, t_ads = 22.2: purity 89.9%, recovery 65.5%, 154.0 kWh/t
  v_feed = 0.813, P_L = 0.072 bar, P_I = 0.124 bar, t_ads = 24.8: purity 93.0%, recovery 47.6%, 140.9 kWh/t
  v_feed = 0.784, P_L = 0.052 bar, P_I = 0.123 bar, t_ads = 22.4: purity 90.2%, recovery 64.3%, 152.3 kWh/t
  v_feed = 0.796, P_L = 0.06 bar, P_I = 0.128 bar, t_ads = 23.4: purity 90.9%, recovery 58.4%, 144.4 kWh/t
  v_feed = 0.791, P_L = 0.056 bar, P_I = 0.126 bar, t_ads = 22.9: purity 90.6%, recovery 60.9%, 147.4 kWh/t
  v_feed = 0.85, P_L = 0.077 bar, P_I = 0.146 bar, t_ads = 28.1: purity 91.8%, recovery 45.4%, 128.1 kWh/t
  v_feed = 0.791, P_L = 0.056 bar, P_I = 0.126 bar, t_ads = 23.0: purity 90.6%, recovery 60.8%, 147.2 kWh/t
  v_feed = 0.792, P_L = 0.057 bar, P_I = 0.126 bar, t_ads = 23.1: purity 90.7%, recovery 60.6%, 146.8 kWh/t
  v_feed = 0.795, P_L = 0.058 bar, P_I = 0.127 bar, t_ads = 23.3: purity 90.8%, recovery 59.9%, 145.6 kWh/t
  v_feed = 0.808, P_L = 0.057 bar, P_I = 0.131 bar, t_ads = 24.4: purity 90.8%, recovery 60.2%, 142.3 kWh/t
  v_feed = 0.947, P_L = 0.058 bar, P_I = 0.174 bar, t_ads = 36.6: purity 92.0%, recovery 57.3%, 119.6 kWh/t
  v_feed = 0.856, P_L = 0.057 bar, P_I = 0.146 bar, t_ads = 28.6: purity 91.2%, recovery 59.5%, 131.9 kWh/t
  v_feed = 0.899, P_L = 0.056 bar, P_I = 0.159 bar, t_ads = 32.4: purity 91.5%, recovery 59.7%, 126.2 kWh/t
  v_feed = 0.925, P_L = 0.054 bar, P_I = 0.167 bar, t_ads = 34.7: purity 91.7%, recovery 59.9%, 123.9 kWh/t
  v_feed = 1.01, P_L = 0.051 bar, P_I = 0.194 bar, t_ads = 42.6: purity 92.4%, recovery 59.0%, 118.5 kWh/t
  v_feed = 0.96, P_L = 0.053 bar, P_I = 0.178 bar, t_ads = 37.8: purity 92.0%, recovery 59.7%, 121.3 kWh/t
  v_feed = 0.979, P_L = 0.052 bar, P_I = 0.184 bar, t_ads = 39.5: purity 92.1%, recovery 59.8%, 120.5 kWh/t
  v_feed = 0.999, P_L = 0.051 bar, P_I = 0.19 bar, t_ads = 41.3: purity 92.3%, recovery 59.9%, 120.0 kWh/t
  v_feed = 1.01, P_L = 0.05 bar, P_I = 0.194 bar, t_ads = 42.5: purity 92.4%, recovery 59.9%, 119.8 kWh/t
  v_feed = 1.01, P_L = 0.05 bar, P_I = 0.194 bar, t_ads = 42.4: purity 92.3%, recovery 59.9%, 119.8 kWh/t
  v_feed = 1.01, P_L = 0.05 bar, P_I = 0.195 bar, t_ads = 42.4: purity 92.3%, recovery 60.0%, 119.8 kWh/t
  v_feed = 1.0, P_L = 0.05 bar, P_I = 0.2 bar, t_ads = 42.8: purity 92.0%, recovery 60.1%, 119.4 kWh/t
  v_feed = 0.984, P_L = 0.05 bar, P_I = 0.222 bar, t_ads = 44.5: purity 91.0%, recovery 60.3%, 117.9 kWh/t
23 evaluations: v_feed = 0.984, P_L = 0.05 bar, P_I = 0.222 bar, t_ads = 44.5
purity 91.0%, recovery 60.3%, 117.9 kWh/t

Comparison of the designs ​

The process is left at the last optimum. The two optimisations side by side:

julia
designs = ["start" => first(max_recovery.history), "max recovery" => max_recovery, "min energy" => min_energy]
rows = [
    "feed velocity [m/s]" => d -> d.values.v_feed,
    "evacuation P_L [bar]" => d -> d.values.P_L / bar,
    "blowdown P_I [bar]" => d -> d.values.P_I / bar,
    "adsorption time [s]" => d -> d.values.t_ads,
    "purity [%]" => d -> 100d.metrics.purity,
    "recovery [%]" => d -> 100d.metrics.recovery,
    "energy [kWh/t]" => d -> d.metrics.specific_energy,
    "productivity [mol/m³/s]" => d -> d.metrics.productivity,
]
println(rpad("", 24), join(lpad.(first.(designs), 14)))
for (label, f) in rows
    println(rpad(label, 24), join([lpad(round(f(d), sigdigits = 3), 14) for d in last.(designs)]))
end
                                 start  max recovery    min energy
feed velocity [m/s]                1.0         0.782         0.984
evacuation P_L [bar]               0.1          0.05          0.05
blowdown P_I [bar]                 0.2         0.123         0.222
adsorption time [s]               15.0          22.2          44.5
purity [%]                        79.5          89.9          91.0
recovery [%]                      36.7          65.5          60.3
energy [kWh/t]                   141.0         154.0         118.0
productivity [mol/m³/s]          0.237         0.474         0.797

Every design tried, on the purity and the energy against the recovery, with the targets dashed and each optimum starred:

julia
fig = Figure(size = (1000, 420))
ax_purity = Axis(fig[1, 1], xlabel = "CO₂ recovery [%]", ylabel = "CO₂ purity [%]")
ax_energy = Axis(fig[1, 2], xlabel = "CO₂ recovery [%]", ylabel = "Specific energy [kWh/t]")
hlines!(ax_purity, [90.0], color = :gray, linestyle = :dash)
vlines!(ax_energy, [60.0], color = :gray, linestyle = :dash)
for (r, color, label) in ((max_recovery, :steelblue, "1. max recovery"), (min_energy, :darkorange, "2. min energy"))
    recovery = [100h.metrics.recovery for h in r.history]
    for (ax, y) in ((ax_purity, [100h.metrics.purity for h in r.history]), (ax_energy, [h.metrics.specific_energy for h in r.history]))
        scatterlines!(ax, recovery, y, color = color, label = label, markersize = 6, linewidth = 0.5)
    end
    scatter!(ax_purity, [100r.metrics.recovery], [100r.metrics.purity], color = color, marker = :star5, markersize = 20)
    scatter!(ax_energy, [100r.metrics.recovery], [r.metrics.specific_energy], color = color, marker = :star5, markersize = 20)
end
axislegend(ax_purity, position = :rb)
fig

Minimising the energy gives up five points of recovery for about a quarter less energy per tonne. It does not raise the evacuation pressure, which stays at its bound: it blows down less deeply and feeds for longer and faster, so that the pump work is spread over more CO₂ captured per cycle. That also raises the productivity, so the same plant would need less adsorbent.

Example on GitHub ​

If you would like to run this example yourself, it can be downloaded from the Mocca.jl GitHub repository.

This page was generated using Literate.jl.