Allow for "sequential control" such as [cc, cccv, rest, ...]

examples/example_sequence_control.jl

@@ -0,0 +1,247 @@
+using BattMo
+using Jutul: si_unit
+
+###
+# Example: Sequence control with CC, Rest, and CCCV steps
+#
+# This example demonstrates the Sequence cycling protocol. The sequence below
+# starts with a short constant-current discharge, rests, and then charges using
+# a CCCV step. Ramp-up is configured globally in the simulation settings and is
+# applied to the current-controlled parts of the CC and CCCV steps.
+###
+
+cell_parameters = load_cell_parameters(; from_default_set = "chen_2020")
+base_simulation_settings = load_simulation_settings(; from_default_set = "p2d")
+base_simulation_settings["TimeStepDuration"] = 360.0
+
+cc = Dict(
+    "Protocol" => "CC",
+    "InitialControl" => "discharging",
+    "DRate" => 1.0,
+    "TotalNumberOfCycles" => 0,
+    "LowerVoltageLimit" => 3.0,
+    "UpperVoltageLimit" => 4.2,
+)
+rest = Dict(
+    "Protocol" => "Rest",
+    "Duration" => 1.0 * si_unit("hour")
+)
+cccv = Dict(
+    "Protocol" => "CCCV",
+    "InitialControl" => "charging",
+    "CRate" => 1.5,
+    "DRate" => 1.0,
+    "TotalNumberOfCycles" => 1,
+    "LowerVoltageLimit" => 3.0,
+    "UpperVoltageLimit" => 3.8,
+    "CurrentChangeLimit" => 1.0e-4,
+    "VoltageChangeLimit" => 1.0e-4,
+)
+
+cycling_protocol = CyclingProtocol(
+    Dict(
+        "Protocol" => "Sequence",
+        "InitialStateOfCharge" => 1.0,
+        "Steps" => [cc, rest, cccv]
+    )
+)
+
+function run_sequence_case(use_rampup::Bool)
+    model_settings = load_model_settings(; from_default_set = "p2d")
+    simulation_settings = deepcopy(base_simulation_settings)
+
+    if use_rampup
+        label = "W rampup"
+    else
+        delete!(model_settings, "RampUp")
+        label = "WO rampup"
+    end
+
+    model = LithiumIonBattery(; model_settings = model_settings)
+    sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings)
+    output = solve(sim; info_level = -1, output_substates = false)
+    substate_output = solve(sim; info_level = -1, output_substates = true)
+
+    return (
+        label = label,
+        output = output,
+        substate_output = substate_output,
+    )
+end
+
+results = [run_sequence_case(use_rampup) for use_rampup in (true, false)]
+
+doplot = true
+if doplot
+
+    using GLMakie
+
+    fig = Figure(size = (900, 760))
+    hour = si_unit("hour")
+    min = si_unit("minute")
+
+    function t_pw_const(times)
+        timestep_time = Float64[]
+        timestep_vals = Float64[]
+        dt = diff(times)
+        for i in 2:length(times)
+            push!(timestep_time, times[i - 1])
+            push!(timestep_time, times[i])
+            push!(timestep_vals, dt[i - 1])
+            push!(timestep_vals, dt[i - 1])
+        end
+        return timestep_time, timestep_vals
+    end
+
+    report_colors = [:dodgerblue3, :darkorange2]
+    substate_colors = [:seagreen4, :purple3]
+
+    ax_voltage = Axis(
+        fig[1, 1],
+        title = "Voltage: report output and substate output",
+        xlabel = "Time / h",
+        ylabel = "Voltage / V",
+    )
+    ax_current = Axis(
+        fig[2, 1],
+        title = "Current: report output and substate output",
+        xlabel = "Time / h",
+        ylabel = "Current / A",
+    )
+    ax_timestep = Axis(
+        fig[3, 1],
+        title = "Accepted solver timestep",
+        xlabel = "Time / h",
+        ylabel = "Accepted solver timestep / min",
+    )
+    linkxaxes!(ax_voltage, ax_current, ax_timestep)
+
+    for (case_index, result) in enumerate(results)
+        label = result.label
+        output = result.output
+        substate_output = result.substate_output
+        report_color = report_colors[case_index]
+        substate_color = substate_colors[case_index]
+
+        scatter!(
+            ax_voltage,
+            output.time_series["Time"] ./ hour,
+            output.time_series["Voltage"],
+            label = "$label: output",
+            color = report_color,
+            markersize = 8,
+        )
+        scatter!(
+            ax_current,
+            output.time_series["Time"] ./ hour,
+            output.time_series["Current"],
+            label = "$label: output",
+            color = report_color,
+            markersize = 8,
+        )
+        output_timestep_time, output_timestep_vals = t_pw_const(output.time_series["Time"])
+        lines!(
+            ax_timestep,
+            output_timestep_time ./ hour,
+            output_timestep_vals ./ min,
+            label = "$label: output",
+            color = report_color,
+            linewidth = 3.0,
+        )
+
+        # output_substates=true shows accepted solver substates, including adaptive
+        # timestep cuts from convergence and control transitions, not only ramp-up steps.
+        time = substate_output.time_series["Time"]
+        scatterlines!(
+            ax_voltage,
+            time ./ hour,
+            substate_output.time_series["Voltage"],
+            label = "$label: substates",
+            color = substate_color,
+            linewidth = 3.0,
+            markersize = 8,
+            marker = :circle,
+        )
+        scatterlines!(
+            ax_current,
+            time ./ hour,
+            substate_output.time_series["Current"],
+            label = "$label: substates",
+            color = substate_color,
+            linewidth = 3.0,
+            markersize = 8,
+            marker = :circle,
+        )
+        substate_timestep_time, substate_timestep_vals = t_pw_const(time)
+        lines!(
+            ax_timestep,
+            substate_timestep_time ./ hour,
+            substate_timestep_vals ./ min,
+            label = "$label: substates",
+            color = substate_color,
+            linewidth = 3.0,
+        )
+    end
+
+    axislegend(ax_voltage, position = :rb)
+    axislegend(ax_current, position = :rb)
+    axislegend(ax_timestep, position = :rt)
+
+    fig_cases = Figure(size = (1100, 760))
+
+    function plot_voltage_current!(position, output, title)
+        time = output.time_series["Time"] ./ hour
+        voltage = output.time_series["Voltage"]
+        current = output.time_series["Current"]
+
+        ax_voltage = Axis(
+            position,
+            title = title,
+            xlabel = "Time / h",
+            ylabel = "Voltage / V",
+        )
+        ax_current = Axis(
+            position,
+            yaxisposition = :right,
+            ylabel = "Current / A",
+        )
+        hidespines!(ax_current, :l, :t, :b)
+        hidexdecorations!(ax_current, grid = false)
+        linkxaxes!(ax_voltage, ax_current)
+
+        voltage_line = lines!(
+            ax_voltage,
+            time,
+            voltage,
+            color = :dodgerblue3,
+            linewidth = 3.0,
+        )
+        current_line = lines!(
+            ax_current,
+            time,
+            current,
+            color = :firebrick3,
+            linewidth = 3.0,
+        )
+        return Legend(
+            position,
+            [voltage_line, current_line],
+            ["Voltage", "Current"],
+            halign = :right,
+            valign = :top,
+            tellwidth = false,
+            tellheight = false,
+        )
+    end
+
+    plot_voltage_current!(fig_cases[1, 1], results[1].output, "W rampup: output")
+    plot_voltage_current!(fig_cases[1, 2], results[1].substate_output, "W rampup: substates")
+    plot_voltage_current!(fig_cases[2, 1], results[2].output, "WO rampup: output")
+    plot_voltage_current!(fig_cases[2, 2], results[2].substate_output, "WO rampup: substates")
+
+    screen1 = GLMakie.Screen()
+    screen2 = GLMakie.Screen()
+    GLMakie.display(screen1, fig)
+    GLMakie.display(screen2, fig_cases)
+
+end

src/BattMo.jl

@@ -214,6 +214,7 @@ include("models/full_battery_models/lithium_ion.jl")
 include("models/full_battery_models/sodium_ion.jl")
 
 include("utils/handy_functions.jl")
+include("utils/curve_analysis.jl")
 
 include("input/loader.jl")
 include("input/defaults.jl")

src/input/defaults/full_simulation_input/chen_2020.json

@@ -203,7 +203,7 @@
             "OutputPath": "",
             "InMemoryReports": 10,
             "ReportLevel": 0,
-            "OutputSubstrates": false
+            "OutputSubstates": false
         }
     }
-}
+}

src/input/defaults/full_simulation_input/chen_2020_p4d.json

@@ -209,7 +209,7 @@
             "OutputPath": "",
             "InMemoryReports": 10,
             "ReportLevel": 0,
-            "OutputSubstrates": false
+            "OutputSubstates": false
         }
     }
-}
+}

src/input/defaults/solver_settings/default.json

@@ -38,6 +38,6 @@
         "OutputPath": "",
         "InMemoryReports": 10,
         "ReportLevel": 0,
-        "OutputSubstrates": false
+        "OutputSubstates": false
     }
-}
+}

src/input/defaults/solver_settings/direct.json

@@ -38,6 +38,6 @@
         "OutputPath": "",
         "InMemoryReports": 10,
         "ReportLevel": 0,
-        "OutputSubstrates": false
+        "OutputSubstates": false
     }
-}
+}

src/input/defaults/solver_settings/iterative.json

@@ -40,6 +40,6 @@
         "OutputPath": "",
         "InMemoryReports": 10,
         "ReportLevel": 0,
-        "OutputSubstrates": false
+        "OutputSubstates": false
     }
-}
+}

src/input/formatter.jl

@@ -479,6 +479,9 @@ function convert_to_parameter_sets(params::AdvancedDictInput)
         cycling_protocol["CRate"] = params["Control"]["CRate"]
         cycling_protocol["CurrentChangeLimit"] = params["Control"]["dIdtLimit"]
         cycling_protocol["VoltageChangeLimit"] = params["Control"]["dEdtLimit"]
+        if haskey(params["Control"], "CVCurrentCutoff")
+            cycling_protocol["CVCurrentCutoff"] = params["Control"]["CVCurrentCutoff"]
+        end
     end
 
     if haskey(model_settings, "ThermalModel")
@@ -751,6 +754,7 @@ function convert_parameter_sets_to_old_input_format(model_settings::ModelSetting
             "upperCutoffVoltage" => get_key_value(cycling_protocol, "UpperVoltageLimit"),
             "dIdtLimit" => get_key_value(cycling_protocol, "CurrentChangeLimit"),
             "dEdtLimit" => get_key_value(cycling_protocol, "VoltageChangeLimit"),
+            "CVCurrentCutoff" => get_key_value(cycling_protocol, "CVCurrentCutoff"),
         ),
         "NegativeElectrode" => Dict(
             "use_normed_current_collector" => false,

src/input/meta_data/cycling_protocol.jl

@@ -5,10 +5,10 @@ function get_cycling_protocol_meta_data()
     meta_data = Dict(
         "Protocol" => Dict(
             "type" => String,
-            "options" => ["CC", "CCCV", "Function", "InputCurrentSeries"],
+            "options" => ["CC", "CCCV", "Function", "InputCurrentSeries", "Rest", "Sequence"],
             "context_type" => "Protocol",
             "context_type_iri" => "https://w3id.org/emmo/domain/electrochemistry#electrochemistry_d3e2d213_d078_4b9a_8beb_62f063e57d69",
-            "description" => """Type of cycling procedure used to cycle a cell. For instance: Constant Current ("CC"), Constant Current - Constant Voltage ("CCCV"), a custom function ("Function"), or a prescribed current time series ("InputCurrentSeries").""",
+            "description" => """Type of cycling procedure used to cycle a cell. For instance: Constant Current ("CC"), Constant Current - Constant Voltage ("CCCV"), a custom function ("Function"), a prescribed current time series ("InputCurrentSeries"), a rest step ("Rest"), or a sequence of supported steps ("Sequence").""",
         ),
         "FunctionName" => Dict(
             "type" => String,
@@ -152,6 +152,14 @@ function get_cycling_protocol_meta_data()
             "type" => Real,
             "unit" => "A·s⁻¹",
         ),
+        "CVCurrentCutoff" => Dict(
+            "min_value" => 0.0,
+            "description" => "Absolute current cutoff used to leave a constant-voltage charge step.",
+            "type" => Real,
+            "unit" => "A",
+            "unit_name" => "emmo:Ampere",
+            "unit_iri" => "https://w3id.org/emmo#Ampere",
+        ),
         "VoltageChangeLimit" => Dict(
             "max_value" => 0.1,
             "min_value" => 0.0001,
@@ -172,6 +180,17 @@ function get_cycling_protocol_meta_data()
             "description" => "Current values for the InputCurrentSeries protocol",
             "unit" => "A",
         ),
+        "Duration" => Dict(
+            "type" => Real,
+            "min_value" => 0.0,
+            "max_value" => 1.0e9,
+            "description" => "Duration of a rest protocol or rest sequence step.",
+            "unit" => "s",
+        ),
+        "Steps" => Dict(
+            "type" => Vector,
+            "description" => "Ordered list of cycling protocol steps used by the Sequence protocol.",
+        ),
     )
 
     return meta_data

src/input/meta_data/model_settings.jl

@@ -76,7 +76,7 @@ function get_model_settings_meta_data()
             "context_type_iri" => "https://w3id.org/emmo/domain/electrochemistry#electrochemistry_rampup",
             "documentation" => "https://battmoteam.github.io/BattMo.jl/dev/manuals/user_guide/ramp_up",
             "category" => "ModelSettings",
-            "description" => """Type of signal of electric current used to initialize the cell simulation. Example: "Sinusoidal".""",
+            "description" => """Type of signal used to ramp control targets. When this setting is present, RampUpTime and RampUpSteps also control extra ministeps used to resolve those ramps. Remove "RampUp" from ModelSettings to disable ramping. The extra time steps when using RampUp also appear in the output. Example: "Sinusoidal".""",
         ),
     )
 

src/input/meta_data/solver_settings.jl

@@ -266,11 +266,11 @@ function get_solver_settings_meta_data()
             "category" => "SolverSettings",
             "description" => "Level of information stored in reports when written to disk.",
         ),
-        "OutputSubstrates" => Dict(
+        "OutputSubstates" => Dict(
             "type" => Bool,
             "variable_name" => "output_substates",
             "category" => "SolverSettings",
-            "description" => "Store substates (between report steps) as field on each state.",
+            "description" => "Return output at every accepted solver substate instead of only at report timesteps.",
         ),
     )
     return meta_data