diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl new file mode 100644 index 000000000..c323f554c --- /dev/null +++ b/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 diff --git a/src/BattMo.jl b/src/BattMo.jl index 0cf90b73a..997f6a3cb 100644 --- a/src/BattMo.jl +++ b/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") diff --git a/src/input/defaults/full_simulation_input/chen_2020.json b/src/input/defaults/full_simulation_input/chen_2020.json index 24e7ecea9..df8b8db36 100644 --- a/src/input/defaults/full_simulation_input/chen_2020.json +++ b/src/input/defaults/full_simulation_input/chen_2020.json @@ -203,7 +203,7 @@ "OutputPath": "", "InMemoryReports": 10, "ReportLevel": 0, - "OutputSubstrates": false + "OutputSubstates": false } } -} \ No newline at end of file +} diff --git a/src/input/defaults/full_simulation_input/chen_2020_p4d.json b/src/input/defaults/full_simulation_input/chen_2020_p4d.json index d74e1f342..dbbae3b3b 100644 --- a/src/input/defaults/full_simulation_input/chen_2020_p4d.json +++ b/src/input/defaults/full_simulation_input/chen_2020_p4d.json @@ -209,7 +209,7 @@ "OutputPath": "", "InMemoryReports": 10, "ReportLevel": 0, - "OutputSubstrates": false + "OutputSubstates": false } } -} \ No newline at end of file +} diff --git a/src/input/defaults/solver_settings/default.json b/src/input/defaults/solver_settings/default.json index 67d117798..27d81d8b8 100644 --- a/src/input/defaults/solver_settings/default.json +++ b/src/input/defaults/solver_settings/default.json @@ -38,6 +38,6 @@ "OutputPath": "", "InMemoryReports": 10, "ReportLevel": 0, - "OutputSubstrates": false + "OutputSubstates": false } -} \ No newline at end of file +} diff --git a/src/input/defaults/solver_settings/direct.json b/src/input/defaults/solver_settings/direct.json index 67d117798..27d81d8b8 100644 --- a/src/input/defaults/solver_settings/direct.json +++ b/src/input/defaults/solver_settings/direct.json @@ -38,6 +38,6 @@ "OutputPath": "", "InMemoryReports": 10, "ReportLevel": 0, - "OutputSubstrates": false + "OutputSubstates": false } -} \ No newline at end of file +} diff --git a/src/input/defaults/solver_settings/iterative.json b/src/input/defaults/solver_settings/iterative.json index f69629b16..0ce45b048 100644 --- a/src/input/defaults/solver_settings/iterative.json +++ b/src/input/defaults/solver_settings/iterative.json @@ -40,6 +40,6 @@ "OutputPath": "", "InMemoryReports": 10, "ReportLevel": 0, - "OutputSubstrates": false + "OutputSubstates": false } -} \ No newline at end of file +} diff --git a/src/input/formatter.jl b/src/input/formatter.jl index b8f163879..04930758e 100644 --- a/src/input/formatter.jl +++ b/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, diff --git a/src/input/meta_data/cycling_protocol.jl b/src/input/meta_data/cycling_protocol.jl index 288c83c7e..bad8cad33 100644 --- a/src/input/meta_data/cycling_protocol.jl +++ b/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 diff --git a/src/input/meta_data/model_settings.jl b/src/input/meta_data/model_settings.jl index baeeb8f4c..48c9f0217 100644 --- a/src/input/meta_data/model_settings.jl +++ b/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".""", ), ) diff --git a/src/input/meta_data/solver_settings.jl b/src/input/meta_data/solver_settings.jl index 3453e53d2..c1f16d66a 100644 --- a/src/input/meta_data/solver_settings.jl +++ b/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 diff --git a/src/input/schemas/get_schema.jl b/src/input/schemas/get_schema.jl index 025473bad..f0c9b51ef 100644 --- a/src/input/schemas/get_schema.jl +++ b/src/input/schemas/get_schema.jl @@ -472,11 +472,14 @@ function get_schema_cycling_protocol(model_settings::ModelSettings) "UpperVoltageLimit" => create_property(parameter_meta, "UpperVoltageLimit"), "InitialControl" => create_property(parameter_meta, "InitialControl"), "CurrentChangeLimit" => create_property(parameter_meta, "CurrentChangeLimit"), + "CVCurrentCutoff" => create_property(parameter_meta, "CVCurrentCutoff"), "VoltageChangeLimit" => create_property(parameter_meta, "VoltageChangeLimit"), "AmbientTemperature" => create_property(parameter_meta, "AmbientTemperature"), "InitialTemperature" => create_property(parameter_meta, "InitialTemperature"), "Times" => create_property(parameter_meta, "Times"), "Currents" => create_property(parameter_meta, "Currents"), + "Duration" => create_property(parameter_meta, "Duration"), + "Steps" => create_property(parameter_meta, "Steps"), ), "required" => ["Protocol"], "allOf" => [ @@ -579,6 +582,24 @@ function get_schema_cycling_protocol(model_settings::ModelSettings) ], ), ), + Dict( + "if" => Dict("properties" => Dict("Protocol" => Dict("const" => "Rest"))), + "then" => Dict( + "required" => [ + "InitialStateOfCharge", + "Duration", + ], + ), + ), + Dict( + "if" => Dict("properties" => Dict("Protocol" => Dict("const" => "Sequence"))), + "then" => Dict( + "required" => [ + "InitialStateOfCharge", + "Steps", + ], + ), + ), ], ) required = schema["required"] @@ -735,7 +756,7 @@ function get_schema_solver_settings() "OutputPath" => create_property(parameter_meta, "OutputPath"), "InMemoryReports" => create_property(parameter_meta, "InMemoryReports"), "ReportLevel" => create_property(parameter_meta, "ReportLevel"), - "OutputSubstrates" => create_property(parameter_meta, "OutputSubstrates"), + "OutputSubstates" => create_property(parameter_meta, "OutputSubstates"), ), ), ), diff --git a/src/input/schemas/lithium_ion/CyclingProtocol.json b/src/input/schemas/lithium_ion/CyclingProtocol.json index 06e7c82b1..448589957 100644 --- a/src/input/schemas/lithium_ion/CyclingProtocol.json +++ b/src/input/schemas/lithium_ion/CyclingProtocol.json @@ -7,11 +7,13 @@ "type": "string", "enum": [ "CC", - "cccv", + "CCCV", "Function", - "InputCurrentSeries" + "InputCurrentSeries", + "Rest", + "Sequence" ], - "description": "Type of cycling procedure used to cycle a cell. E.g., Constant Current ('CC'), Constant Current - Constant Voltage ('cccv'), a custom Function, or a prescribed current time series ('InputCurrentSeries')." + "description": "Type of cycling procedure used to cycle a cell. E.g., Constant Current ('CC'), Constant Current - Constant Voltage ('CCCV'), a custom Function, a prescribed current time series ('InputCurrentSeries'), a rest step ('Rest'), or an ordered sequence of supported steps ('Sequence')." }, "TotalTime": { "type": "number", @@ -105,6 +107,14 @@ "context_type_iri": "https://w3id.org/emmo/domain/electrochemistry#electrochemistry_71f10616_15eb_4dc4_bc8d_ffaac3838af2", "unit": "A·s⁻¹" }, + "CVCurrentCutoff": { + "type": "number", + "minimum": 0.0, + "description": "Absolute current cutoff used to leave a constant-voltage charge step.", + "unit": "A", + "unit_name": "emmo:Ampere", + "unit_iri": "https://w3id.org/emmo#Ampere" + }, "VoltageChangeLimit": { "type": "number", "minimum": 0.0001, @@ -136,6 +146,21 @@ "minItems": 1, "description": "Current values corresponding to each time point in the InputCurrentSeries protocol, in amperes. Positive values indicate discharge; negative values indicate charge.", "unit": "A" + }, + "Duration": { + "type": "number", + "minimum": 0.0, + "maximum": 1000000000.0, + "description": "Duration of a rest protocol or rest sequence step, in seconds.", + "unit": "s" + }, + "Steps": { + "type": "array", + "items": { + "type": "object" + }, + "minItems": 1, + "description": "Ordered list of cycling protocol steps used by the Sequence protocol." } }, "required": [ @@ -275,6 +300,34 @@ "UpperVoltageLimit" ] } + }, + { + "if": { + "properties": { + "Protocol": { + "const": "Rest" + } + } + }, + "then": { + "required": [ + "Duration" + ] + } + }, + { + "if": { + "properties": { + "Protocol": { + "const": "Sequence" + } + } + }, + "then": { + "required": [ + "Steps" + ] + } } ] } diff --git a/src/input/schemas/parameters_meta_data.json b/src/input/schemas/parameters_meta_data.json index 3b988b317..e143cea9e 100644 --- a/src/input/schemas/parameters_meta_data.json +++ b/src/input/schemas/parameters_meta_data.json @@ -125,9 +125,11 @@ "CC", "CCCV", "Function", - "InputCurrentSeries" + "InputCurrentSeries", + "Rest", + "Sequence" ], - "description": "Type of cycling procedure used to cycle a cell. E.g., Constant Current ('CC'), Constant Current - Constant Voltage ('CCCV'), or a prescribed current time series ('InputCurrentSeries')." + "description": "Type of cycling procedure used to cycle a cell. E.g., Constant Current ('CC'), Constant Current - Constant Voltage ('CCCV'), a custom Function, a prescribed current time series ('InputCurrentSeries'), a rest step ('Rest'), or an ordered sequence of supported steps ('Sequence')." }, "FunctionName": { "type": "string", @@ -579,6 +581,14 @@ "context_type_iri": "https://w3id.org/emmo/domain/electrochemistry#electrochemistry_71f10616_15eb_4dc4_bc8d_ffaac3838af2", "unit": "A·s⁻¹" }, + "CVCurrentCutoff": { + "type": "number", + "minimum": 0.0, + "description": "Absolute current cutoff used to leave a constant-voltage charge step.", + "unit": "A", + "unit_name": "emmo:Ampere", + "unit_iri": "https://w3id.org/emmo#Ampere" + }, "VoltageChangeLimit": { "type": "number", "minimum": 0.0001, @@ -614,4 +624,4 @@ }, "required": [], "additionalProperties": true -} \ No newline at end of file +} diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index f5d79b996..352283c9a 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -7,7 +7,9 @@ export SimpleCVPolicy, CyclingCVPolicy, OperationalMode, - InputCurrentPolicy + InputCurrentPolicy, + RestPolicy, + SequencePolicy ################################ # Define the operational modes # @@ -16,19 +18,17 @@ export @enum OperationalMode cc_discharge1 cc_discharge2 cc_charge1 cv_charge2 rest discharge charging discharging none function getSymbol(ctrlType::OperationalMode) - if ctrlType == cc_discharge1 - symb = :cc_discharge1 + return :cc_discharge1 elseif ctrlType == cc_discharge2 - symb = :cc_discharge2 + return :cc_discharge2 elseif ctrlType == cc_charge1 - symb = :cc_charge1 + return :cc_charge1 elseif ctrlType == cv_charge2 - symb = :cv_charge2 + return :cv_charge2 + else + error("Unsupported CCCV control type: $ctrlType") end - - return symb - end ############################################# @@ -69,6 +69,55 @@ Jutul.number_of_cells(::CurrentAndVoltageDomain) = 1 # Types for the different policies # #################################### +abstract type AbstractSequenceStep end + +struct PolicyStep{P <: AbstractPolicy} <: AbstractSequenceStep + policy::P +end + +mutable struct RestPolicy{R} <: AbstractPolicy + duration::R +end + +mutable struct SequencePolicy{R} <: AbstractPolicy + steps::Vector{AbstractSequenceStep} + ImaxDischarge::R + ImaxCharge::R + use_ramp_up::Bool + rampup_time::R +end + +function SequencePolicy( + steps::Vector{AbstractSequenceStep}, + ImaxDischarge::R, + ImaxCharge::R, + use_ramp_up::Bool, + rampup_time::R, + ) where {R <: Real} + return SequencePolicy{R}(steps, ImaxDischarge, ImaxCharge, use_ramp_up, rampup_time) +end + +function SequencePolicy( + steps::AbstractVector, + ImaxDischarge::Real, + ImaxCharge::Real, + use_ramp_up::Bool, + rampup_time::Real, + ) + normalized_steps = AbstractSequenceStep[] + for step in steps + if step isa AbstractSequenceStep + push!(normalized_steps, step) + elseif step isa AbstractPolicy + push!(normalized_steps, PolicyStep(step)) + else + error("Sequence step $(typeof(step)) must be an AbstractPolicy or AbstractSequenceStep") + end + end + T = promote_type(typeof(ImaxDischarge), typeof(ImaxCharge), typeof(rampup_time)) + return SequencePolicy{T}(normalized_steps, T(ImaxDischarge), T(ImaxCharge), use_ramp_up, T(rampup_time)) +end + ## A policy is used to compute the next control from the current control and state """ Simple constant current. Stops when lower cut-off value is reached @@ -115,10 +164,6 @@ mutable struct SimpleCVPolicy{R} <: AbstractPolicy end end -""" No policy means that the control is kept fixed throughout the simulation -""" -struct NoPolicy <: AbstractPolicy end - """ Function Policy @@ -128,7 +173,7 @@ struct FunctionPolicy <: AbstractPolicy function FunctionPolicy(function_name::String; file_path::Union{Nothing, String} = nothing) current_function = setup_function_from_function_name(function_name; file_path = file_path) - return new{}(current_function) + return new(current_function) end end @@ -176,10 +221,12 @@ mutable struct CyclingCVPolicy{R, I} <: AbstractPolicy upperCutoffVoltage::R dIdtLimit::R dEdtLimit::R + cvCurrentCutoff::Union{R, Missing} initialControl::OperationalMode numberOfCycles::I tolerances::Any use_ramp_up::Bool + rampup_time::R current_function::Any end @@ -193,8 +240,10 @@ function CyclingCVPolicy( numberOfCycles; ImaxDischarge = 0 * lowerCutoffVoltage, ImaxCharge = 0 * lowerCutoffVoltage, - use_ramp_up::Bool, + use_ramp_up::Bool = false, + rampup_time = zero(lowerCutoffVoltage), current_function = missing, + cv_current_cutoff = missing, ) if initialControl == "charging" @@ -205,6 +254,10 @@ function CyclingCVPolicy( error("InitialControl $initialControl not recognized") end + if isnothing(cv_current_cutoff) + cv_current_cutoff = missing + end + tolerances = ( cc_discharge1 = 1.0e-4, cc_discharge2 = 0.9, @@ -219,10 +272,12 @@ function CyclingCVPolicy( upperCutoffVoltage, dIdtLimit, dEdtLimit, + cv_current_cutoff, initialControl, numberOfCycles, tolerances, use_ramp_up, + rampup_time, current_function, ) end @@ -280,10 +335,19 @@ function Jutul.select_parameters!( return S[:ImaxCharge] = ImaxCharge() end +function Jutul.select_parameters!( + S, + system::CurrentAndVoltageSystem{SequencePolicy{R}}, + model::SimulationModel, + ) where {R} + S[:ImaxDischarge] = ImaxDischarge() + return S[:ImaxCharge] = ImaxCharge() +end -########################################################################################################### -# Definition of the controller and some basic utility functions. The controller will be part of the state # -########################################################################################################### + +####################################### +# Types for the different controllers # +####################################### ## A controller provides the information to exert the current control @@ -331,6 +395,12 @@ mutable struct CcCvController{R, I <: Integer} <: Controller numberOfCycles::I dEdt::Union{R, Missing} dIdt::Union{R, Missing} + ramp_start_time::R + ramp_duration::R + ramp_start_target::R # target value at beginning of ramp (actual current or actual voltage) + ramp_end_target::R # target the controller should reach at the end of the ramp (eg ImaxDischarge, lowerCutoffVoltage etc) + ramp_target_is_voltage::Bool # current or voltage + ramp_active::Bool # to keep track if the controller is in the transition end @@ -338,7 +408,7 @@ function CcCvController() maincontroller = SimpleControllerCV() - return CcCvController(maincontroller, 0, missing, missing) + return CcCvController(maincontroller, 0, missing, missing, 0.0, 0.0, 0.0, 0.0, false, false) end @@ -361,7 +431,7 @@ function Base.setproperty!(c::CcCvController, f::Symbol, v) end end -@inline function Jutul.numerical_type(x::CCController{R, I}) where {R, I} +@inline function Jutul.numerical_type(x::CCController{I, R}) where {I, R} return R end @@ -391,125 +461,121 @@ InputCurrentController() = InputCurrentController(0.0, 0.0, false) return R end +## RestController -""" -Function to create (deep) copy of simple controller -""" -function copyController!(cv_copy::CCController, cv::CCController) - - cv_copy.numberOfCycles = cv.numberOfCycles - cv_copy.target = cv.target - cv_copy.time = cv.time - cv_copy.target_is_voltage = cv.target_is_voltage - return cv_copy.ctrlType = cv.ctrlType - -end - - -""" -Function to create (deep) copy of function controller -""" -function copyController!(cv_copy::FunctionController, cv::FunctionController) - - cv_copy.target = cv.target - cv_copy.time = cv.time - return cv_copy.target_is_voltage = cv.target_is_voltage - -end - -""" -Function to create (deep) copy of simple controller -""" -function copyController!(cv_copy::SimpleControllerCV, cv::SimpleControllerCV) - - cv_copy.target = cv.target - cv_copy.time = cv.time - cv_copy.target_is_voltage = cv.target_is_voltage - return cv_copy.ctrlType = cv.ctrlType - +mutable struct RestController{R} <: Controller + target::R + time::R + target_is_voltage::Bool + ctrlType::OperationalMode end -""" -Function to create (deep) copy of CC-CV controller -""" -function copyController!(cv_copy::CcCvController, cv::CcCvController) - - copyController!(cv_copy.maincontroller, cv.maincontroller) - return cv_copy.numberOfCycles = cv.numberOfCycles +RestController() = RestController(0.0, 0.0, false, rest) +@inline function Jutul.numerical_type(x::RestController{R}) where {R} + return R end -""" -Overload function to copy simple controller -""" -function Base.copy(cv::FunctionController) - - cv_copy = FunctionController() - copyController!(cv_copy, cv) - - return cv_copy +## SequenceController +mutable struct SequenceController{R} <: Controller + # Keep data for all the controllers to get a concrete controller for Jutul + step_index::Int + step_start_time::R + target::R + # Global elapsed controller time. step_start_time tracks the start of the + # active sequence step, so local step time is time - step_start_time. + time::R + target_is_voltage::Bool + ctrlType::Any + numberOfCycles::Int + dEdt::Union{R, Missing} + dIdt::Union{R, Missing} + ramp_start_time::R + ramp_duration::R + ramp_start_target::R + ramp_end_target::R + ramp_target_is_voltage::Bool + ramp_active::Bool end +# Construct a sequence controller initialized at the first step with no active ramp. +SequenceController() = SequenceController(1, 0.0, 0.0, 0.0, false, missing, 0, missing, missing, 0.0, 0.0, 0.0, 0.0, false, false) -""" -Overload function to copy simple controller -""" -function Base.copy(cv::CCController) - - cv_copy = CCController() - copyController!(cv_copy, cv) - - return cv_copy - +# Return the numeric scalar type used by the sequence controller. +@inline function Jutul.numerical_type(x::SequenceController{R}) where {R} + return R end -""" -Overload function to copy simple controller -""" -function Base.copy(cv::SimpleControllerCV) - cv_copy = SimpleControllerCV() - copyController!(cv_copy, cv) - - return cv_copy +############################################## +# Copy helpers for the different controllers # +############################################## +function copyController!(dst::T, src::T) where {T <: Controller} + for field in fieldnames(T) + setfield!(dst, field, getfield(src, field)) + end + return dst end -""" -Overload function to copy CC-CV controller -""" -function Base.copy(cv::CcCvController) - - cv_copy = CcCvController() - copyController!(cv_copy, cv) - - return cv_copy - +function copy_controller_field(dst_value, src_value) + if dst_value isa Bool && src_value isa Bool + return src_value + elseif dst_value isa Number && src_value isa Number + return zero(dst_value) + src_value + else + return src_value + end end -""" -Function to create (deep) copy of input current controller -""" -function copyController!(cv_copy::InputCurrentController, cv::InputCurrentController) - - cv_copy.target = cv.target - cv_copy.time = cv.time - return cv_copy.target_is_voltage = cv.target_is_voltage - +function copyController!(dst::T, src::S) where {T <: Controller, S <: Controller} + # Adjoint tracing may copy CCController{..., Float64} into CCController{..., Dual}. + fieldnames(T) == fieldnames(S) || error("Cannot copy controller $(S) into $(T)") + for field in fieldnames(T) + setfield!(dst, field, copy_controller_field(getfield(dst, field), getfield(src, field))) + end + return dst +end + +function copyController!(dst::CcCvController, src::CcCvController) + # Need its own copyController since it contains a nested mutable controller + copyController!(dst.maincontroller, src.maincontroller) + dst.numberOfCycles = src.numberOfCycles + dst.dEdt = src.dEdt + dst.dIdt = src.dIdt + dst.ramp_start_time = src.ramp_start_time + dst.ramp_duration = src.ramp_duration + dst.ramp_start_target = src.ramp_start_target + dst.ramp_end_target = src.ramp_end_target + dst.ramp_target_is_voltage = src.ramp_target_is_voltage + dst.ramp_active = src.ramp_active + return dst +end + +function Base.copy(c::T) where {T <: Controller} + return T((getfield(c, field) for field in fieldnames(T))...) +end + +function Base.copy(c::CcCvController{R, I}) where {R, I} + return CcCvController( + copy(c.maincontroller), + c.numberOfCycles, + c.dEdt, + c.dIdt, + c.ramp_start_time, + c.ramp_duration, + c.ramp_start_target, + c.ramp_end_target, + c.ramp_target_is_voltage, + c.ramp_active, + ) end -""" -Overload function to copy input current controller -""" -function Base.copy(cv::InputCurrentController) - - cv_copy = InputCurrentController() - copyController!(cv_copy, cv) - - return cv_copy +################################################# +# Minimum output variables for the control model # +################################################# -end """ We add the controller in the output """ @@ -586,17 +652,197 @@ function getInitCurrent(policy::InputCurrentPolicy) return policy.current_function(policy.times[1]) end +function getInitCurrent(policy::RestPolicy) + return zero(policy.duration) +end + +# Return the initial current prescribed by the first step in a sequence. +function getInitCurrent(policy::SequencePolicy) + return getInitCurrent(sequence_step_policy(policy.steps[1])) +end + function getInitCurrent(model::CurrentAndVoltageModel) return getInitCurrent(model.system.policy) end +# Unwrap a policy step to get the policy that should control the active step. +sequence_step_policy(step::PolicyStep) = step.policy +# Allow already-unwrapped policies to be used as sequence steps. +sequence_step_policy(policy::AbstractPolicy) = policy + +##################################### +# Helpers for sequence policy steps # +##################################### + +# Return the policy object for the currently active sequence step. +function active_sequence_policy(policy::SequencePolicy, controller::SequenceController) + return sequence_step_policy(policy.steps[controller.step_index]) +end + +# Initialize the shared sequence controller fields for a rest step. +function initialize_sequence_step_controller!(controller::SequenceController, step_policy::RestPolicy) + controller.target = 0.0 + controller.target_is_voltage = false + controller.ctrlType = rest + controller.numberOfCycles = 0 + controller.dEdt = missing + return controller.dIdt = missing +end + +# Initialize the shared sequence controller fields for a CC step. +function initialize_sequence_step_controller!(controller::SequenceController, step_policy::CCPolicy) + if step_policy.initialControl == "discharging" + controller.ctrlType = "discharging" + controller.target = step_policy.ImaxDischarge + elseif step_policy.initialControl == "charging" + controller.ctrlType = "charging" + controller.target = -step_policy.ImaxCharge + else + error("Initial control $(step_policy.initialControl) is not recognized") + end + controller.target_is_voltage = false + controller.numberOfCycles = 0 + controller.dEdt = missing + return controller.dIdt = missing +end + +# Initialize the shared sequence controller fields for a CCCV step. +function initialize_sequence_step_controller!(controller::SequenceController, step_policy::CyclingCVPolicy) + if step_policy.initialControl == discharging + controller.ctrlType = cc_discharge1 + elseif step_policy.initialControl == charging + controller.ctrlType = cc_charge1 + else + error("Initial control $(step_policy.initialControl) is not recognized") + end + controller.target = 0.0 + controller.target_is_voltage = false + controller.numberOfCycles = 0 + controller.dEdt = missing + return controller.dIdt = missing +end + +# A rest step is complete once its step-local duration has elapsed. +function sequence_step_complete(policy::RestPolicy, controller, state) + return controller.time - controller.step_start_time >= policy.duration +end + +# A CC step is complete when its voltage limit or cycle count condition is reached. +function sequence_step_complete(policy::CCPolicy, controller, state) + if policy.numberOfCycles == 0 + if policy.initialControl == "charging" + return state.ElectricPotential[1] >= policy.upperCutoffVoltage + elseif policy.initialControl == "discharging" + return state.ElectricPotential[1] <= policy.lowerCutoffVoltage + else + error("Initial control $(policy.initialControl) is not recognized") + end + else + return controller.numberOfCycles >= policy.numberOfCycles + end +end + +# A CCCV step is complete when its requested cycle count has been reached. +function sequence_step_complete(policy::CyclingCVPolicy, controller, state) + return controller.numberOfCycles >= policy.numberOfCycles +end + +# A sequence is complete when the active step index has advanced past the final step. +function sequence_complete(policy::SequencePolicy, controller::SequenceController) + return controller.step_index > length(policy.steps) +end + +# The nominal target for a rest step is zero current. +function sequence_nominal_target(policy::RestPolicy, controller::SequenceController) + return (target = 0.0, target_is_voltage = false) +end + +# The nominal target for a CC step is the current setpoint for the active direction. +function sequence_nominal_target(policy::CCPolicy, controller::SequenceController) + if controller.ctrlType == "discharging" + target = policy.ImaxDischarge + elseif controller.ctrlType == "charging" + target = -policy.ImaxCharge + else + error("ctrlType $(controller.ctrlType) not recognized") + end + return (target = target, target_is_voltage = false) +end + +# The nominal target for a CCCV step is the current or voltage setpoint for its active mode. +function sequence_nominal_target(policy::CyclingCVPolicy, controller::SequenceController) + if controller.ctrlType == cc_discharge1 + target = policy.ImaxDischarge + target_is_voltage = false + elseif controller.ctrlType == cc_discharge2 + target = 0.0 + target_is_voltage = false + elseif controller.ctrlType == cc_charge1 + target = -policy.ImaxCharge + target_is_voltage = false + elseif controller.ctrlType == cv_charge2 + target = policy.upperCutoffVoltage + target_is_voltage = true + else + error("ctrlType $(controller.ctrlType) not recognized") + end + return (target = target, target_is_voltage = target_is_voltage) +end + +# Advance to the next sequence step and start a transition ramp when requested. +function advance_sequence_step!(controller::SequenceController, policy::SequencePolicy) + previous_target = controller.target + previous_target_is_voltage = controller.target_is_voltage + + controller.step_index += 1 + controller.step_start_time = controller.time + if !sequence_complete(policy, controller) + initialize_sequence_step_controller!(controller, active_sequence_policy(policy, controller)) + nominal = sequence_nominal_target(active_sequence_policy(policy, controller), controller) + if policy.use_ramp_up && policy.rampup_time > 0 && previous_target_is_voltage == nominal.target_is_voltage && previous_target != nominal.target + controller.ramp_start_time = controller.time + controller.ramp_duration = policy.rampup_time + controller.ramp_start_target = previous_target + controller.ramp_end_target = nominal.target + controller.ramp_target_is_voltage = nominal.target_is_voltage + controller.ramp_active = true + controller.target = previous_target + controller.target_is_voltage = previous_target_is_voltage + else + controller.ramp_active = false + end + else + controller.ramp_active = false + end + return controller +end + +# Apply an active sequence transition ramp to the controller target. +function apply_sequence_transition_ramp!(controller::SequenceController) + if controller.ramp_active + elapsed = controller.time - controller.ramp_start_time + if elapsed < controller.ramp_duration + delta = controller.ramp_end_target - controller.ramp_start_target + controller.target = controller.ramp_start_target + currentFun(elapsed, delta, controller.ramp_duration) + controller.target_is_voltage = controller.ramp_target_is_voltage + else + controller.target = controller.ramp_end_target + controller.target_is_voltage = controller.ramp_target_is_voltage + controller.ramp_active = false + end + end + return controller +end + ###################################################### # Setup the initial policy from the input parameters # ###################################################### +setup_initial_control_policy!(policy::AbstractPolicy, input, parameters) = nothing + function setup_initial_control_policy!(policy::CCPolicy, input, parameters) cycling_protocol = input.cycling_protocol @@ -624,7 +870,8 @@ function setup_initial_control_policy!(policy::CCPolicy, input, parameters) if haskey(cycling_protocol, "UpperVoltageLimit") policy.upperCutoffVoltage = cycling_protocol["UpperVoltageLimit"] - elseif haskey(cycling_protocol, "LowerVoltageLimit") + end + if haskey(cycling_protocol, "LowerVoltageLimit") policy.lowerCutoffVoltage = cycling_protocol["LowerVoltageLimit"] end policy.ImaxCharge = only(parameters[:Control][:ImaxCharge]) @@ -633,10 +880,6 @@ function setup_initial_control_policy!(policy::CCPolicy, input, parameters) end -function setup_initial_control_policy!(policy::FunctionPolicy, input, parameters) - -end - function setup_initial_control_policy!(policy::SimpleCVPolicy, input, parameters) cycling_protocol = input.cycling_protocol @@ -676,6 +919,7 @@ function setup_initial_control_policy!(policy::CyclingCVPolicy, input, parameter cFun(time) = currentFun(time, Imax, tup) policy.current_function = cFun + policy.rampup_time = tup end @@ -691,6 +935,15 @@ function setup_initial_control_policy!(policy::InputCurrentPolicy, input, parame # Nothing additional needs to be set up from parameters. end +function setup_initial_control_policy!(policy::RestPolicy, input, parameters) + return nothing +end + +# Sequence steps are fully configured before setup, so no additional initial policy data is needed. +function setup_initial_control_policy!(policy::SequencePolicy, input, parameters) + return nothing +end + ################################### # Special primary variable update # ################################### @@ -722,6 +975,29 @@ function Jutul.update_primary_variable!(state, p::CurrentVar, state_symbol, mode end +# Bound sequence current updates by the largest current used by any step in the sequence. +function Jutul.update_primary_variable!(state, p::CurrentVar, state_symbol, model::P, dx, w) where {P <: CurrentAndVoltageModel{<:SequencePolicy}} + + entity = associated_entity(p) + active = active_entities(model.domain, entity, for_variables = true) + v = state[state_symbol] + + ImaxDischarge = model.system.policy.ImaxDischarge + ImaxCharge = model.system.policy.ImaxCharge + Imax = max(ImaxCharge, ImaxDischarge) + + abs_max = 0.2 * Imax + rel_max = relative_increment_limit(p) + maxval = maximum_value(p) + minval = minimum_value(p) + scale = variable_scale(p) + return @inbounds for i in eachindex(active) + a_i = active[i] + v[a_i] = update_value(v[a_i], w * dx[i], abs_max, rel_max, minval, maxval, scale) + end + +end + ####################################### # Helper functions for control switch # ####################################### @@ -746,10 +1022,7 @@ function setupRegionSwitchFlags(policy::Union{CyclingCVPolicy, CCPolicy}, state, end - E = only(state.ElectricPotential) - I = only(state.Current) - if ctrlType == cc_discharge1 || ctrlType == "discharging" @@ -774,13 +1047,19 @@ function setupRegionSwitchFlags(policy::Union{CyclingCVPolicy, CCPolicy}, state, elseif ctrlType == cv_charge2 - dIdt = state.Controller.dIdt - if !ismissing(dIdt) - before = abs(dIdt) > dIdtMin * (1 + tol) - after = abs(dIdt) < dIdtMin * (1 - tol) + if !ismissing(policy.cvCurrentCutoff) + I = abs(only(state.Current)) + before = I > policy.cvCurrentCutoff + after = I < policy.cvCurrentCutoff else - before = false - after = false + dIdt = state.Controller.dIdt + if !ismissing(dIdt) + before = abs(dIdt) > dIdtMin * (1 + tol) + after = abs(dIdt) < dIdtMin * (1 - tol) + else + before = false + after = false + end end else @@ -807,7 +1086,19 @@ function check_constraints(model, storage) ctrlType = state[:Controller].ctrlType ctrlType0 = state0[:Controller].ctrlType - if policy isa CyclingCVPolicy + if policy isa SequencePolicy + if sequence_complete(policy, controller) + return true + end + policy = active_sequence_policy(policy, controller) + end + + if policy isa RestPolicy + return true + elseif policy isa CyclingCVPolicy + if controller.ramp_active + return true + end nextCtrlType = getNextCtrlTypecccv(ctrlType0) elseif policy isa CCPolicy if ctrlType == "discharging" @@ -869,6 +1160,19 @@ function Jutul.update_values!(old::InputCurrentController, new::InputCurrentCont end +function Jutul.update_values!(old::RestController, new::RestController) + + return copyController!(old, new) + +end + +# Copy the sequence controller when Jutul updates controller values. +function Jutul.update_values!(old::SequenceController, new::SequenceController) + + return copyController!(old, new) + +end + """ In addition to update the values in all primary variables, we need also to update the values in the controller. We do that by specializing the method perform_step_solve_impl! """ @@ -883,21 +1187,10 @@ end """ We need to add the specific treatment of the controller variables """ -function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{CyclingCVPolicy{T1, T2}}, T3, T4}) where {T1, T2, T3, T4} - - invoke( - reset_state_to_previous_state!, - Tuple{ - typeof(storage), - SimulationModel, - }, +function Jutul.reset_state_to_previous_state!( storage, - model, - ) - return copyController!(storage.state[:Controller], storage.state0[:Controller]) -end - -function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{CCPolicy{T1}}, T3, T4}) where {T1, T3, T4} + model::SimulationModel{CurrentAndVoltageDomain, <:CurrentAndVoltageSystem, T3, T4}, + ) where {T3, T4} invoke( reset_state_to_previous_state!, @@ -911,21 +1204,6 @@ function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{Cu return copyController!(storage.state[:Controller], storage.state0[:Controller]) end -function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{InputCurrentPolicy{T1}}, T3, T4}) where {T1, T3, T4} - - invoke( - reset_state_to_previous_state!, - Tuple{ - typeof(storage), - SimulationModel, - }, - storage, - model, - ) - return copyController!(storage.state[:Controller], storage.state0[:Controller]) - -end - function update_controller!(state, state0, policy::AbstractPolicy, dt) @@ -969,9 +1247,64 @@ function update_control_type_in_controller!(state, state0, policy::SimpleCVPolic end +##################################### +# Helpers for CCCV transition ramps # +##################################### + +function cccv_nominal_target(policy::CyclingCVPolicy, ctrlType) + if ctrlType == cc_discharge1 + target = policy.ImaxDischarge + target_is_voltage = false + elseif ctrlType == cc_discharge2 + target = 0.0 + target_is_voltage = false + elseif ctrlType == cc_charge1 + target = -policy.ImaxCharge + target_is_voltage = false + elseif ctrlType == cv_charge2 + target = policy.upperCutoffVoltage + target_is_voltage = true + else + error("ctrlType $ctrlType not recognized") + end + return (target = target, target_is_voltage = target_is_voltage) +end + +function start_cccv_transition_ramp!(controller, policy::CyclingCVPolicy, ctrlType, E, I) + nominal = cccv_nominal_target(policy, ctrlType) + controller.ramp_start_time = controller.time + controller.ramp_duration = policy.rampup_time + controller.ramp_start_target = nominal.target_is_voltage ? E : I + controller.ramp_end_target = nominal.target + controller.ramp_target_is_voltage = nominal.target_is_voltage + controller.ramp_active = policy.use_ramp_up && policy.rampup_time > 0 && controller.ramp_start_target != controller.ramp_end_target + return controller +end + +function apply_cccv_transition_ramp!(controller) + if controller.ramp_active + elapsed = controller.time - controller.ramp_start_time + if elapsed < controller.ramp_duration + delta = controller.ramp_end_target - controller.ramp_start_target + controller.target = controller.ramp_start_target + currentFun(elapsed, delta, controller.ramp_duration) + controller.target_is_voltage = controller.ramp_target_is_voltage + else + controller.target = controller.ramp_end_target + controller.target_is_voltage = controller.ramp_target_is_voltage + controller.ramp_active = false + end + end + return controller +end + """ Implementation of the cycling CC-CV policy """ + +################################################## +# Control-type updates for the different policies # +################################################## + function update_control_type_in_controller!(state, state0, policy::CyclingCVPolicy, dt) E = only(value(state[:ElectricPotential])) @@ -987,6 +1320,10 @@ function update_control_type_in_controller!(state, state0, policy::CyclingCVPoli ctrlType0 = state0.Controller.ctrlType + if state0.Controller.ramp_active + return controller.ctrlType = ctrlType0 + end + nextCtrlType = getNextCtrlTypecccv(ctrlType0) rsw00 = setupRegionSwitchFlags(policy, state0, ctrlType0) @@ -1038,6 +1375,10 @@ function update_control_type_in_controller!(state, state0, policy::CyclingCVPoli end + if ctrlType != ctrlType0 && !controller.ramp_active + start_cccv_transition_ramp!(controller, policy, ctrlType, E, I) + end + return controller.ctrlType = ctrlType end @@ -1120,10 +1461,27 @@ function update_control_type_in_controller!(state, state0, policy::CCPolicy, dt) end +function update_control_type_in_controller!(state, state0, policy::RestPolicy, dt) + controller = state.Controller + controller.time = state0.Controller.time + dt + controller.target_is_voltage = false + return controller.ctrlType = rest +end -############################################################# -# Functions to update the values in the controller in state # -############################################################# +# Delegate sequence control-type updates to the currently active step policy. +function update_control_type_in_controller!(state, state0, policy::SequencePolicy, dt) + controller = state.Controller + if sequence_complete(policy, controller) + return nothing + end + step_policy = active_sequence_policy(policy, controller) + return update_control_type_in_controller!(state, state0, step_policy, dt) +end + + +################################################## +# Control-type update for input-current policies # +################################################## """ Implementation of the InputCurrentPolicy. @@ -1182,6 +1540,10 @@ end # Once the controller has been assigned the given control, we adjust the target value which is used in the equation # assembly +################################################### +# Target-value updates for the different policies # +################################################### + function update_values_in_controller!(state, policy::CCPolicy) controller = state.Controller @@ -1348,7 +1710,8 @@ function update_values_in_controller!(state, policy::CyclingCVPolicy) controller.target_is_voltage = target_is_voltage - return controller.target = target + controller.target = target + return apply_cccv_transition_ramp!(controller) end @@ -1375,6 +1738,102 @@ function update_values_in_controller!(state, policy::InputCurrentPolicy) end +############################################# +# Target-value updates for sequence steps # +############################################# + +# Update a sequence rest step target using the standalone rest target logic. +function update_sequence_values_in_controller!(state, policy::RestPolicy) + return update_values_in_controller!(state, policy) +end + +# Update a sequence CC step target using step-local time for any ramp function. +function update_sequence_values_in_controller!(state, policy::CCPolicy) + controller = state.Controller + ctrlType = controller.ctrlType + cf = policy.current_function + + if controller.numberOfCycles == 0 && controller.ctrlType == policy.initialControl && !ismissing(cf) + local_time = controller.time - controller.step_start_time + I_t = cf isa Real ? cf : cf(local_time) + + if ctrlType == "discharging" + target = I_t + elseif ctrlType == "charging" + target = -I_t + else + error("ctrlType $ctrlType not recognized") + end + else + if ctrlType == "discharging" + target = policy.ImaxDischarge + elseif ctrlType == "charging" + target = -policy.ImaxCharge + else + error("ctrlType $ctrlType not recognized") + end + end + + controller.target_is_voltage = false + return controller.target = target +end + +# Update a sequence CCCV step target using step-local time for any ramp function. +function update_sequence_values_in_controller!(state, policy::CyclingCVPolicy) + controller = state[:Controller] + ctrlType = controller.ctrlType + cf = policy.current_function + + if ctrlType == cc_discharge1 + if controller.numberOfCycles == 0 && controller.ctrlType == policy.initialControl && !ismissing(cf) + local_time = controller.time - controller.step_start_time + target = cf isa Real ? cf : cf(local_time) + else + target = policy.ImaxDischarge + end + target_is_voltage = false + + elseif ctrlType == cc_discharge2 + target = 0.0 + target_is_voltage = false + + elseif ctrlType == cc_charge1 + if controller.numberOfCycles == 0 && controller.ctrlType == policy.initialControl && !ismissing(cf) + local_time = controller.time - controller.step_start_time + target = -(cf isa Real ? cf : cf(local_time)) + else + target = -policy.ImaxCharge + end + target_is_voltage = false + + elseif ctrlType == cv_charge2 + target = policy.upperCutoffVoltage + target_is_voltage = true + + else + error("ctrlType $ctrlType not recognized") + end + + controller.target_is_voltage = target_is_voltage + return controller.target = target +end + +function update_values_in_controller!(state, policy::RestPolicy) + controller = state.Controller + controller.target_is_voltage = false + return controller.target = 0.0 +end + +# Update the active sequence step target and apply any sequence transition ramp. +function update_values_in_controller!(state, policy::SequencePolicy) + controller = state[:Controller] + if sequence_complete(policy, controller) + return nothing + end + update_sequence_values_in_controller!(state, active_sequence_policy(policy, controller)) + return apply_sequence_transition_ramp!(controller) +end + ############################# # Assembly of the equations # ############################# @@ -1404,7 +1863,6 @@ end function Jutul.update_equation_in_entity!(v, i, state, state0, eq::CurrentEquation, model, dt, ldisc = local_discretization(eq, i)) - # Sign is strange here due to cross term? I = only(state.Current) phi = only(state.ElectricPotential) @@ -1417,92 +1875,90 @@ end # Function to update the controller part in state after convergence # ##################################################################### -""" Update after convergence. Here, we copy the controller to state0 and count the total number of cycles in case of CyclingCVPolicy -""" -function Jutul.update_after_step!(storage, domain::CurrentAndVoltageDomain, model::CurrentAndVoltageModel, dt, forces; time = NaN) - - ctrl = storage.state[:Controller] - - policy = model.system.policy - - return if policy isa CyclingCVPolicy - - initctrl = policy.initialControl - - ctrlType = ctrl.ctrlType - - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles +function update_cycle_count!(ctrl, ctrl0, initialControl::OperationalMode) + ctrlType = ctrl.ctrlType + ctrlType0 = ctrl0.ctrlType + ncycles = ctrl0.numberOfCycles - copyController!(storage.state0[:Controller], ctrl) - - if initctrl == charging - - if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) - ncycles = ncycles + 1 - end + if initialControl == charging + if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) + ncycles += 1 + end + elseif initialControl == discharging + if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) + ncycles += 1 + end + else + error("Initial control $initialControl is not recognized") + end - elseif initctrl == discharging + ctrl.numberOfCycles = ncycles + return ctrl +end - if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) - ncycles = ncycles + 1 - end +function update_cycle_count!(ctrl, ctrl0, initialControl::String) + ctrlType = ctrl.ctrlType + ctrlType0 = ctrl0.ctrlType + ncycles = ctrl0.numberOfCycles + if initialControl == "charging" + if ctrlType0 == "discharging" && ctrlType == "charging" + ncycles += 1 end + elseif initialControl == "discharging" + if ctrlType0 == "charging" && ctrlType == "discharging" + ncycles += 1 + end + else + error("Initial control $initialControl is not recognized") + end - ctrl.numberOfCycles = ncycles + ctrl.numberOfCycles = ncycles + return ctrl +end - elseif policy isa SimpleCVPolicy +""" Update after convergence. Copy the controller to state0 and update cycle counts. """ +function Jutul.update_after_step!(storage, domain::CurrentAndVoltageDomain, model::CurrentAndVoltageModel, dt, forces; time = NaN) - copyController!(storage.state0[:Controller], ctrl) + ctrl = storage.state[:Controller] + ctrl0 = storage.state0[:Controller] + policy = model.system.policy - elseif policy isa FunctionPolicy + if policy isa CyclingCVPolicy + update_cycle_count!(ctrl, ctrl0, policy.initialControl) + return copyController!(ctrl0, ctrl) - copyController!(storage.state0[:Controller], ctrl) + elseif policy isa SimpleCVPolicy || policy isa FunctionPolicy || policy isa InputCurrentPolicy || policy isa RestPolicy + return copyController!(ctrl0, ctrl) elseif policy isa CCPolicy + if policy.numberOfCycles > 0 + update_cycle_count!(ctrl, ctrl0, policy.initialControl) + end + return copyController!(ctrl0, ctrl) - if policy.numberOfCycles == 0 - copyController!(storage.state0[:Controller], ctrl) - - else - - initctrl = policy.initialControl - - ctrlType = ctrl.ctrlType - - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles - - copyController!(storage.state0[:Controller], ctrl) - - if initctrl == "charging" - - if ctrlType0 == "discharging" && ctrlType == "charging" - ncycles = ncycles + 1 - end - - elseif initctrl == "discharging" - - if ctrlType0 == "charging" && ctrlType == "discharging" - ncycles = ncycles + 1 - end + elseif policy isa SequencePolicy + if sequence_complete(policy, ctrl) + return copyController!(ctrl0, ctrl) + end - end + step_policy = active_sequence_policy(policy, ctrl) - ctrl.numberOfCycles = ncycles + if step_policy isa CyclingCVPolicy + update_cycle_count!(ctrl, ctrl0, step_policy.initialControl) + elseif step_policy isa CCPolicy && step_policy.numberOfCycles > 0 + update_cycle_count!(ctrl, ctrl0, step_policy.initialControl) end - elseif policy isa InputCurrentPolicy + if sequence_step_complete(step_policy, ctrl, storage.state) + advance_sequence_step!(ctrl, policy) + end - copyController!(storage.state0[:Controller], ctrl) + return copyController!(ctrl0, ctrl) else - error("Policy $(typeof(policy)) not recognized") - end - end ######################################################################## @@ -1596,6 +2052,25 @@ function Jutul.initialize_extra_state_fields!(state, ::Any, model::CurrentAndVol target, time = promote(target, time) state[:Controller] = InputCurrentController(target, time, target_is_voltage) + elseif policy isa RestPolicy + + time = 0.0 + target = 0.0 + target_is_voltage = false + + target, time = promote(target, time) + state[:Controller] = RestController(target, time, target_is_voltage, rest) + + elseif policy isa SequencePolicy + + time = 0.0 + target = 0.0 + target_is_voltage = false + + target, time = promote(target, time) + state[:Controller] = SequenceController(1, time, target, time, target_is_voltage, missing, 0, missing, missing, time, 0.0, target, target, target_is_voltage, false) + initialize_sequence_step_controller!(state[:Controller], active_sequence_policy(policy, state[:Controller])) + end end @@ -1636,9 +2111,20 @@ end # Helper function to compute control value # ############################################ -""" -sineup rampup function -""" +function currentFun(t::Real, inputI::Real, tup::Real = 0.1) + + t, inputI, tup, val = promote(t, inputI, tup, 0.0) + + if t <= tup + val = sineup(0.0, inputI, 0.0, tup, t) + else + val = inputI + end + + return val + +end + function sineup(y1, y2, x1, x2, x) #SINEUP Creates a sine ramp function # diff --git a/src/models/formatter.jl b/src/models/formatter.jl index abdbc0246..fd4b2e7ea 100644 --- a/src/models/formatter.jl +++ b/src/models/formatter.jl @@ -451,6 +451,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 return ( @@ -718,6 +721,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, diff --git a/src/models/full_battery_models/intercalation_battery.jl b/src/models/full_battery_models/intercalation_battery.jl index 7e6437263..6c28f3b1d 100644 --- a/src/models/full_battery_models/intercalation_battery.jl +++ b/src/models/full_battery_models/intercalation_battery.jl @@ -212,6 +212,9 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) protocol = cycling_protocol["Protocol"] + cap = min(computeElectrodeCapacity(model_neam, :NegativeElectrodeActiveMaterial), computeElectrodeCapacity(model_peam, :PositiveElectrodeActiveMaterial)) + con = Constants() + if protocol == "CC" initial_control = cycling_protocol["InitialControl"] @@ -222,7 +225,6 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) CRate = get(cycling_protocol, "CRate", 0.0) # Capacity goes into actual realized rates, so check the type for AD/promotion - cap = min(computeElectrodeCapacity(model_neam, :NegativeElectrodeActiveMaterial), computeElectrodeCapacity(model_peam, :PositiveElectrodeActiveMaterial)) T_i = promote_type(typeof(DRate), typeof(CRate), typeof(cap), T) policy = CCPolicy( @@ -237,6 +239,11 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) elseif protocol == "CCCV" + cv_current_cutoff = get(cycling_protocol, "CVCurrentCutoff", missing) + if isnothing(cv_current_cutoff) + cv_current_cutoff = missing + end + policy = CyclingCVPolicy( cycling_protocol["LowerVoltageLimit"], cycling_protocol["UpperVoltageLimit"], @@ -245,6 +252,7 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) cycling_protocol["InitialControl"], cycling_protocol["TotalNumberOfCycles"]; use_ramp_up = use_ramp_up, + cv_current_cutoff = cv_current_cutoff, ) elseif protocol == "Function" @@ -263,6 +271,88 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) policy = InputCurrentPolicy(times, currents, lower_v, upper_v) + elseif protocol == "Rest" + + policy = RestPolicy(cycling_protocol["Duration"]) + + elseif protocol == "Sequence" + + steps = AbstractSequenceStep[] + ImaxDischarge = zero(cap / con.hour) + ImaxCharge = zero(cap / con.hour) + + for step in cycling_protocol["Steps"] + step_protocol = step["Protocol"] + + if step_protocol == "CC" + DRate = get(step, "DRate", 0.0) + CRate = get(step, "CRate", 0.0) + Idis = (cap / con.hour) * DRate + Ichg = (cap / con.hour) * CRate + ImaxDischarge = max(ImaxDischarge, Idis) + ImaxCharge = max(ImaxCharge, Ichg) + + step_policy = CCPolicy( + get(step, "TotalNumberOfCycles", 0), + step["InitialControl"], + step["LowerVoltageLimit"], + step["UpperVoltageLimit"], + use_ramp_up; + ImaxDischarge = Idis, + ImaxCharge = Ichg, + T = promote_type(typeof(Idis), typeof(Ichg), typeof(step["LowerVoltageLimit"]), typeof(step["UpperVoltageLimit"])), + ) + if use_ramp_up + Imax = step["InitialControl"] == "charging" ? Ichg : Idis + tup = Float64(input.simulation_settings["RampUpTime"]) + step_policy.current_function = time -> currentFun(time, Imax, tup) + end + push!(steps, PolicyStep(step_policy)) + + elseif step_protocol == "CCCV" + DRate = step["DRate"] + CRate = step["CRate"] + Idis = (cap / con.hour) * DRate + Ichg = (cap / con.hour) * CRate + ImaxDischarge = max(ImaxDischarge, Idis) + ImaxCharge = max(ImaxCharge, Ichg) + + cv_current_cutoff = get(step, "CVCurrentCutoff", missing) + if isnothing(cv_current_cutoff) + cv_current_cutoff = missing + end + + step_policy = CyclingCVPolicy( + step["LowerVoltageLimit"], + step["UpperVoltageLimit"], + step["CurrentChangeLimit"], + step["VoltageChangeLimit"], + step["InitialControl"], + get(step, "TotalNumberOfCycles", 1); + ImaxDischarge = Idis, + ImaxCharge = Ichg, + use_ramp_up = use_ramp_up, + rampup_time = use_ramp_up ? Float64(input.simulation_settings["RampUpTime"]) : zero(Idis), + cv_current_cutoff = cv_current_cutoff, + ) + if use_ramp_up + Imax = step["InitialControl"] == "charging" ? Ichg : Idis + tup = Float64(input.simulation_settings["RampUpTime"]) + step_policy.current_function = time -> currentFun(time, Imax, tup) + end + push!(steps, PolicyStep(step_policy)) + + elseif step_protocol == "Rest" + push!(steps, PolicyStep(RestPolicy(step["Duration"]))) + + else + error("Sequence step protocol $step_protocol not recognized.") + end + end + + rampup_time = use_ramp_up ? Float64(input.simulation_settings["RampUpTime"]) : zero(ImaxDischarge) + policy = SequencePolicy(steps, ImaxDischarge, ImaxCharge, use_ramp_up, rampup_time) + else error("controlPolicy not recognized.") @@ -822,6 +912,30 @@ function set_parameters( parameters[:Control] = setup_parameters(multimodel[:Control]) + elseif protocol == "Rest" + + parameters[:Control] = setup_parameters(multimodel[:Control]) + + elseif protocol == "Sequence" + + cap = computeCellCapacity(multimodel) + con = Constants() + + ImaxDischarge = 0.0 + ImaxCharge = 0.0 + + for step in cycling_protocol["Steps"] + if step["Protocol"] == "CC" || step["Protocol"] == "CCCV" + ImaxDischarge = max(ImaxDischarge, (cap / con.hour) * get(step, "DRate", 0.0)) + ImaxCharge = max(ImaxCharge, (cap / con.hour) * get(step, "CRate", 0.0)) + end + end + + prm_control[:ImaxDischarge] = ImaxDischarge + prm_control[:ImaxCharge] = ImaxCharge + + parameters[:Control] = setup_parameters(multimodel[:Control], prm_control) + else error("control policy $controlPolicy not recognized") end @@ -1055,25 +1169,18 @@ end function setup_initial_state(input, model::IntercalationBattery) if hasproperty(input, :initial_state) && !isnothing(input.initial_state) + # Copy to keep the physical state state = deepcopy(input.initial_state) - # Reset controller - if haskey(state, :Control) - ctrl = state[:Control] - if haskey(ctrl, :Controller) - controller = ctrl[:Controller] - for fd in fieldnames(typeof(controller)) - value = getfield(controller, fd) - if !isa(value, String) - value = zero(value) - end - setfield!(controller, fd, value) - end - end - ctrl[:Current] = [getInitCurrent(model.multimodel[:Control])] - end + # Reset control + fresh_input = merge(input, (; initial_state = nothing)) + fresh_state = setup_initial_state(fresh_input, model) + + state[:Control][:Controller] = deepcopy(fresh_state[:Control][:Controller]) + state[:Control][:Current] = deepcopy(fresh_state[:Control][:Current]) return state + end multimodel = model.multimodel diff --git a/src/output/output_format.jl b/src/output/output_format.jl index 913b47888..052867fa5 100644 --- a/src/output/output_format.jl +++ b/src/output/output_format.jl @@ -564,7 +564,7 @@ function get_output_metrics( controller = states[1][:Control][:Controller] - if isa(controller, FunctionController) || isa(controller, InputCurrentController) + if isa(controller, FunctionController) || isa(controller, InputCurrentController) || isa(controller, RestController) available_quantities = Dict() else cycle_array = hasproperty(states[1][:Control][:Controller], :numberOfCycles) ? [state[:Control][:Controller].numberOfCycles for state in states] : nothing diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 648c027cf..fdaaca5c4 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -2,9 +2,36 @@ export Simulation export solve export setup_config +struct ControlRampTimestepSelector <: Jutul.AbstractTimestepSelector + timestep::Float64 +end + +Jutul.pick_first_timestep(::ControlRampTimestepSelector, sim, config, dT, forces) = Inf +Jutul.pick_cut_timestep(::ControlRampTimestepSelector, sim, config, dt, dT, forces, reports, cut_count) = dt +Jutul.valid_timestep(::ControlRampTimestepSelector, dt) = dt + +function Jutul.pick_next_timestep(sel::ControlRampTimestepSelector, sim, config, dt_prev, dT, forces, reports, current_reports, step_index, new_step) + controller = Jutul.get_simulator_storage(sim).state.Control.Controller + if hasfield(typeof(controller), :ramp_active) && controller.ramp_active + remaining = controller.ramp_start_time + controller.ramp_duration - controller.time + if remaining > 0 + return min(sel.timestep, remaining) + end + end + return dT +end + +function control_rampup_time(policy) + if policy isa Union{CyclingCVPolicy, SequencePolicy} + return policy.rampup_time + else + return 0.0 + end +end + """ - abstract type AbstractSimulation + abstract type AbstractSimulation Abstract type for Simulation structs. Subtypes of `AbstractSimulation` represent specific simulation configurations. """ @@ -284,6 +311,7 @@ end info_level = 0, end_report = info_level > -1, include_initial_state = false, + output_substates = false, kwargs...) Solves a battery `Simulation` problem by executing the simulation workflow defined in `solve_simulation`. @@ -298,6 +326,9 @@ Solves a battery `Simulation` problem by executing the simulation workflow defin prepended to the output time series (`Time`, `Voltage`, `Current`, capacity, etc.). This eliminates the gap between the initial condition and the first solver output, improving plots and RMSE comparisons when restarting from a saved state. Default is `false`. +- `output_substates::Bool` (optional): If `true`, the returned BattMo output is sampled at + every accepted solver substate instead of only at report timesteps. The output format is + unchanged; time-dependent arrays simply contain more time points. Default is `false`. - `kwargs...`: Additional keyword arguments forwarded to `solve_simulation`. # Behavior @@ -326,6 +357,7 @@ function solve( solver_settings = get_default_solver_settings(problem.model), logger = nothing, include_initial_state = false, + output_substates::Bool = false, kwargs..., ) @@ -358,6 +390,7 @@ function solve( solver_settings, logger, include_initial_state, + output_substates, validate, accept_invalid, kwargs..., @@ -405,6 +438,7 @@ function solve_simulation( solver_settings, logger = nothing, include_initial_state = false, + output_substates::Bool = false, validate::Bool = true, accept_invalid::Bool = false, kwargs..., @@ -436,12 +470,20 @@ function solve_simulation( ), parameters; solver_settings, + rampup_steps = get(simulation_settings, "RampUpSteps", nothing), + rampup_reference_timestep = get(simulation_settings, "TimeStepDuration", nothing), validate, accept_invalid, logger, kwargs..., ) + # Jutul stores accepted solver substates when output_substates is enabled. + output_substates = output_substates || cfg[:output_substates] + if output_substates + cfg[:output_substates] = true + end + # Setup hook if given hook = get(kwargs, :hook, nothing) @@ -457,7 +499,12 @@ function solve_simulation( end - jutul_states, jutul_reports = simulate(state0, simulator, timesteps; forces = forces, config = cfg, kwargs...) + tmp_states, jutul_reports = simulate(state0, simulator, timesteps; forces = forces, config = cfg, kwargs...) + + output_states = tmp_states + if output_substates + output_states, _, _ = Jutul.expand_to_ministeps(tmp_states, jutul_reports[eachindex(tmp_states)]) + end # perfom decoupled thermal simulation if needed @@ -468,13 +515,20 @@ function solve_simulation( thermal_parameters = thermal_sim.parameters thermal_model = thermal_sim.model - decoupled_thermal_states = solve_decoupled_thermal_model(model, cycling_protocol, thermal_parameters, parameters, thermal_model, jutul_states, global_maps, timesteps) + thermal_timesteps = if output_substates + diff(vcat(0.0, [state[:Control][:Controller].time for state in output_states])) + else + timesteps + end + decoupled_thermal_states = solve_decoupled_thermal_model(model, cycling_protocol, thermal_parameters, parameters, thermal_model, output_states, global_maps, thermal_timesteps) - # Add decoupled thermal states to the jutul_states - for i in eachindex(jutul_states) - jutul_states[i][:ThermalModel] = decoupled_thermal_states[i][:state] - jutul_reports[i][:ThermalModel] = decoupled_thermal_states[i][:report] + # Add decoupled thermal states to the output states + for i in eachindex(output_states) + output_states[i][:ThermalModel] = decoupled_thermal_states[i][:state] + if !output_substates + jutul_reports[i][:ThermalModel] = decoupled_thermal_states[i][:report] + end end @@ -482,7 +536,7 @@ function solve_simulation( jutul_output = ( - states = jutul_states, + states = output_states, reports = jutul_reports, solver_configuration = cfg, multimodel = model.multimodel, @@ -503,7 +557,7 @@ function solve_simulation( # Optionally include the initial state as the first data point if include_initial_state jutul_output_for_ts = ( - states = vcat([deepcopy(state0)], jutul_states), + states = vcat([deepcopy(state0)], output_states), reports = jutul_reports, solver_configuration = cfg, multimodel = model.multimodel, @@ -631,7 +685,7 @@ function kwarg_dict(; kwargs...) "OutputReports" => get(kwargs, :output_reports, nothing), "InMemoryReports" => get(kwargs, :in_memory_reports, nothing), "ReportLevel" => get(kwargs, :report_level, nothing), - "OutputSubstrates" => get(kwargs, :output_substates, nothing), + "OutputSubstates" => get(kwargs, :output_substates, nothing), ), ) @@ -670,6 +724,8 @@ function solver_configuration( parameters; use_model_scaling::Bool = true, solver_settings, + rampup_steps = nothing, + rampup_reference_timestep = nothing, logger = nothing, validate::Bool = true, accept_invalid::Bool = false, @@ -708,6 +764,7 @@ function solver_configuration( linear_solver_dict = solver_settings["LinearSolver"] output = solver_settings["Output"] verbose = solver_settings["Verbose"] + output_substates_setting = get(output, "OutputSubstates", false) relaxation = non_linear_solver["Relaxation"] if relaxation == "NoRelaxation" @@ -720,11 +777,17 @@ function solver_configuration( timestep_selector = non_linear_solver["TimeStepSelectors"] if timestep_selector == "TimestepSelector" - timesel = [TimestepSelector()] + timesel = Jutul.AbstractTimestepSelector[TimestepSelector()] else error("Timestep selector $(timestep_selector) not recognized. Only 'TimestepSelector' is currently implemented.") end + control_ramp_time = control_rampup_time(multimodel[:Control].system.policy) + if !isnothing(rampup_steps) && rampup_steps > 0 && control_ramp_time > 0 + ramp_dt = Float64(control_ramp_time) / rampup_steps + push!(timesel, ControlRampTimestepSelector(ramp_dt)) + end + if linear_solver_dict["Method"] == "UserDefined" linear_solver = overwritten_settings.linear_solver else @@ -763,9 +826,15 @@ function solver_configuration( output_path = output["OutputPath"] == "" ? nothing : output["OutputPath"], in_memory_reports = output["InMemoryReports"], report_level = output["ReportLevel"], - output_substates = output["OutputSubstrates"], + output_substates = output_substates_setting, ) + if !isnothing(rampup_reference_timestep) && !isnothing(rampup_steps) && rampup_steps > 0 && control_ramp_time > 0 + ramp_dt = Float64(control_ramp_time) / rampup_steps + ramp_decrease = ramp_dt / Float64(rampup_reference_timestep) + cfg[:timestep_max_decrease] = min(cfg[:timestep_max_decrease], ramp_decrease) + end + if !isempty(non_linear_solver["Tolerances"]) cfg[:tolerances] = non_linear_solver["Tolerances"] end @@ -789,13 +858,13 @@ function solver_configuration( end end - if multimodel[:Control].system.policy isa CyclingCVPolicy || multimodel[:Control].system.policy isa CCPolicy + if multimodel[:Control].system.policy isa CyclingCVPolicy || multimodel[:Control].system.policy isa CCPolicy || multimodel[:Control].system.policy isa RestPolicy || multimodel[:Control].system.policy isa SequencePolicy if multimodel[:Control].system.policy isa CyclingCVPolicy - cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) - elseif multimodel[:Control].system.policy isa CCPolicy && multimodel[:Control].system.policy.numberOfCycles > 0 cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) + elseif multimodel[:Control].system.policy isa SequencePolicy + cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) end function post_hook(done, report, sim, dt, forces, max_iter, cfg) @@ -841,6 +910,11 @@ function solver_configuration( report[:stopnow] = false end end + elseif multimodel[:Control].system.policy isa RestPolicy + report[:stopnow] = s.state.Control.Controller.time >= m[:Control].system.policy.duration + + elseif multimodel[:Control].system.policy isa SequencePolicy + report[:stopnow] = sequence_complete(m[:Control].system.policy, s.state.Control.Controller) end return (done, report) @@ -1014,24 +1088,32 @@ function setup_timesteps( if cycling_protocol["TotalNumberOfCycles"] == 0 if cycling_protocol["InitialControl"] == "discharging" - CRate = cycling_protocol["DRate"] + rate = cycling_protocol["DRate"] else - CRate = cycling_protocol["CRate"] + rate = cycling_protocol["CRate"] end con = Constants() - totalTime = 1.1 * con.hour / CRate + totalTime = 1.1 * con.hour / rate dt = simulation_settings["TimeStepDuration"] if haskey(input.model_settings, "RampUp") nr = simulation_settings["RampUpSteps"] + timesteps = compute_rampup_timesteps(totalTime, dt, nr) else - nr = 1 + # Set to false to use only regular dt steps for no-ramp CC. + use_initial_half_step = true + timesteps = use_initial_half_step ? [dt / 2] : Float64[] + remaining_time = totalTime - sum(timesteps) + n = Int64(floor(remaining_time / dt)) + append!(timesteps, repeat([dt], n)) + dt_final = totalTime - sum(timesteps) + if dt_final > 0 + push!(timesteps, dt_final) + end end - timesteps = compute_rampup_timesteps(totalTime, dt, nr) - else ncycles = cycling_protocol["TotalNumberOfCycles"] @@ -1059,9 +1141,14 @@ function setup_timesteps( totalTime = ncycles * 2.5 * (1 * con.hour / CRate + 1 * con.hour / DRate) dt = simulation_settings["TimeStepDuration"] - n = Int64(floor(totalTime / dt)) - timesteps = repeat([dt], n) + if haskey(input.model_settings, "RampUp") + nr = simulation_settings["RampUpSteps"] + timesteps = compute_rampup_timesteps(totalTime, dt, nr) + else + n = Int64(floor(totalTime / dt)) + timesteps = repeat([dt], n) + end elseif protocol == "Function" @@ -1077,6 +1164,38 @@ function setup_timesteps( timesteps = diff(series_times) timesteps = timesteps[timesteps .> 0.0] + elseif protocol == "Rest" + + totalTime = cycling_protocol["Duration"] + dt = simulation_settings["TimeStepDuration"] + if haskey(input.model_settings, "RampUp") + nr = simulation_settings["RampUpSteps"] + timesteps = compute_rampup_timesteps(totalTime, dt, nr) + else + n = Int64(floor(totalTime / dt)) + timesteps = repeat([dt], n) + dt_final = totalTime - sum(timesteps) + if dt_final > 0 + push!(timesteps, dt_final) + end + end + + elseif protocol == "Sequence" + + timesteps = Float64[] + for step in cycling_protocol["Steps"] + step_protocol = CyclingProtocol( + merge( + Dict( + "InitialStateOfCharge" => get(cycling_protocol, "InitialStateOfCharge", 0.5), + "TotalNumberOfCycles" => get(step, "TotalNumberOfCycles", step["Protocol"] == "CCCV" ? 1 : 0), + ), + step, + ) + ) + append!(timesteps, setup_timesteps((cycling_protocol = step_protocol, simulation_settings = simulation_settings, model_settings = input.model_settings))) + end + else error("Protocol $protocol not recognized") @@ -1116,21 +1235,3 @@ function compute_rampup_timesteps(time::Real, dt::Real, n::Integer = 8) return dT end - -#################### -# Current function # -#################### - -function currentFun(t::Real, inputI::Real, tup::Real = 0.1) - - t, inputI, tup, val = promote(t, inputI, tup, 0.0) - - if t <= tup - val = sineup(0.0, inputI, 0.0, tup, t) - else - val = inputI - end - - return val - -end diff --git a/src/utils/curve_analysis.jl b/src/utils/curve_analysis.jl new file mode 100644 index 000000000..fff80ecc5 --- /dev/null +++ b/src/utils/curve_analysis.jl @@ -0,0 +1,158 @@ +export trapz, rmse, compute_dqdv + +""" + trapz(x, y) + +Integrate `y` with respect to `x` using the trapezoidal rule. +""" +function trapz(x, y) + length(x) == length(y) || throw(DimensionMismatch("x and y must have the same number of elements")) + + length(x) >= 2 || throw(ArgumentError("x must contain at least two points")) + + s = zero((x[2] - x[1]) * (y[1] + y[2])) + + @inbounds for i in 1:(length(x) - 1) + s += (x[i + 1] - x[i]) * (y[i] + y[i + 1]) + end + + return s / 2 +end + + +""" + rmse(x0, y0, x1, y1) + +Compute the root mean squared error using trapezoidal integration`. +""" + +function rmse( + x1, y1, + x2, y2; + extrap::Bool = false, + truncate::Bool = false, + ) + + length(x1) == length(y1) || throw(DimensionMismatch("x1 and y1 must have the same number of elements")) + length(x2) == length(y2) || throw(DimensionMismatch("x2 and y2 must have the same number of elements")) + + function strictly_increasing(x) + @inbounds for i in 2:length(x) + x[i] > x[i - 1] || return false + end + return true + end + + strictly_increasing(x1) || throw(ArgumentError("x1 must be strictly increasing")) + strictly_increasing(x2) || throw(ArgumentError("x2 must be strictly increasing")) + + # Optionally truncate both datasets to their common x-range. + # This preserves ordering and does not sort. + if truncate + xmin = max(first(x1), first(x2)) + xmax = min(last(x1), last(x2)) + + xmin < xmax || throw(ArgumentError("datasets have no overlapping x-range")) + + idx1 = findall(x -> xmin <= x <= xmax, x1) + idx2 = findall(x -> xmin <= x <= xmax, x2) + + x1 = x1[idx1] + y1 = y1[idx1] + x2 = x2[idx2] + y2 = y2[idx2] + end + + length(x1) >= 2 || throw(ArgumentError("x1 must contain at least two points after truncation")) + length(x2) >= 2 || throw(ArgumentError("x2 must contain at least two points after truncation")) + + if length(x1) > length(x2) # Could also use median spacing + # x1 is finer: interpolate y2 onto x1 + x = x1 + ya = y1 + # yb = Jutul.linear_interp(x2, y2, x) + interp = Jutul.get_1d_interpolator(x2, y2) + yb = interp.(x) + else + # x2 is finer: interpolate y1 onto x2 + x = x2 + ya = y2 + # yb = Jutul.linear_interp(x1, y1, x) + interp = Jutul.get_1d_interpolator(x1, y1) + yb = interp.(x) + end + + all(isfinite, yb) || throw(ArgumentError("Interpolated values are not finite. Check that the finer grid is within the range of the coarser grid, or use extrap=true or truncate=true.")) + + # Compute L2 error on the finer grid. + l2sq = trapz(x, @. (ya - yb)^2) + l2 = sqrt(l2sq) + wl2 = sqrt(l2sq / (last(x) - first(x))) + + return wl2 +end + + +function gaussian_smooth(values; sigma = 1.0, truncate = 4.0) + length(values) == 1 && return copy(values) + radius = round(Int, truncate * sigma) + offsets = collect(-radius:radius) + weights = exp.(-0.5 .* (offsets ./ sigma) .^ 2) + weights ./= sum(weights) + n = length(values) + + function reflected_index(i) + while i < 1 || i > n + i = i < 1 ? 1 - i : 2 * n + 1 - i + end + return i + end + + return [ + sum(weights[j] * values[reflected_index(i + offsets[j])] for j in eachindex(weights)) + for i in eachindex(values) + ] +end + +""" + compute_dqdv(q, v; bin_size = nothing, nbins = nothing, smooth = true) + +Compute a differential-capacity curve using an equal-width voltage histogram. +Either `bin_size` or `nbins` must be provided. If both are provided, `nbins` +takes precedence. Gaussian smoothing with `sigma = 1` is enabled by default. + +Returns `(voltage, dqdv)`. For decreasing-voltage curves, both arrays are +reversed and the differential capacity is negated. +""" +function compute_dqdv(q, v; bin_size = nothing, nbins = nothing, smooth = true) + length(q) == length(v) || throw(ArgumentError("q and v must have equal length.")) + length(q) >= 2 || throw(ArgumentError("At least two points are required.")) + all(isfinite, q) && all(isfinite, v) || throw(ArgumentError("q and v must contain only finite values.")) + isnothing(bin_size) && isnothing(nbins) && throw(ArgumentError("Either bin_size or nbins must be provided.")) + + vmin, vmax = extrema(v) + vmax > vmin || throw(ArgumentError("Voltage values must span a non-zero interval.")) + if isnothing(nbins) + bin_size > 0 || throw(ArgumentError("bin_size must be positive.")) + nbins = floor(Int, (vmax - vmin) / bin_size) + end + nbins isa Integer && nbins > 0 || throw(ArgumentError("nbins must be a positive integer.")) + + bin_width = (vmax - vmin) / nbins + counts = zeros(Float64, nbins) + for voltage in v + index = voltage == vmax ? nbins : floor(Int, (voltage - vmin) / bin_width) + 1 + counts[index] += 1 + end + + counts ./= length(v) * bin_width + counts .*= maximum(q) - minimum(q) + smooth && (counts = gaussian_smooth(counts)) + + bin_centers = collect(range(vmin + bin_width / 2, step = bin_width, length = nbins)) + if v[end] > v[1] + return (bin_centers, counts) + else + return (reverse(bin_centers), -reverse(counts)) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 5661c289d..bdaaa8e25 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using Test include("test_voltage_breakdown.jl") include("test_input_current_series.jl") include("test_initial_state.jl") +include("test_sequence_protocols.jl") include("test_sei.jl") include("test_loader.jl") # include("test_matlab_input.jl") diff --git a/test/test_capacity_calculation.jl b/test/test_capacity_calculation.jl index 6fa2d9c8c..65a598be1 100644 --- a/test/test_capacity_calculation.jl +++ b/test/test_capacity_calculation.jl @@ -9,7 +9,9 @@ using Test cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") cycling_protocol = load_cycling_protocol(; from_default_set = "cccv") - model_setup = LithiumIonBattery() + ms = load_model_settings(; from_default_set = "p2d") + delete!(ms, "RampUp") + model_setup = LithiumIonBattery(; model_settings = ms) sim = Simulation(model_setup, cell_parameters, cycling_protocol) diff --git a/test/test_cycling_protocols.jl b/test/test_cycling_protocols.jl index 515b3c53e..6ed305f28 100644 --- a/test/test_cycling_protocols.jl +++ b/test/test_cycling_protocols.jl @@ -1,6 +1,117 @@ using BattMo using Test +@testset "CCCV switch ramp-up" begin + policy = BattMo.CyclingCVPolicy( + 3.0, + 4.0, + 1.0e-4, + 1.0e-4, + "charging", + 1; + ImaxDischarge = 1.0, + ImaxCharge = 1.0, + use_ramp_up = true, + rampup_time = 100.0, + ) + + controller0 = BattMo.CcCvController() + controller0.ctrlType = BattMo.cc_charge1 + controller0.time = 0.0 + controller0.target = -1.0 + + controller = copy(controller0) + + state0 = ( + Controller = controller0, + ElectricPotential = [4.0], + Current = [-1.0], + ) + state = ( + Controller = controller, + ElectricPotential = [4.01], + Current = [-0.8], + ) + + BattMo.update_control_type_in_controller!(state, state0, policy, 10.0) + + @test controller.ctrlType == BattMo.cv_charge2 + @test controller.ramp_active + @test controller.ramp_target_is_voltage + + BattMo.update_values_in_controller!(state, policy) + + @test controller.target_is_voltage + @test controller.target ≈ 4.01 + + controller.time = 60.0 + BattMo.update_values_in_controller!(state, policy) + + @test 4.0 < controller.target < 4.01 +end + +@testset "CCCV CVCurrentCutoff" begin + policy = BattMo.CyclingCVPolicy( + 3.0, + 4.0, + 1.0e-4, + 1.0e-4, + "charging", + 1; + ImaxDischarge = 1.0, + ImaxCharge = 1.0, + cv_current_cutoff = 0.05, + ) + + controller = BattMo.CcCvController() + controller.ctrlType = BattMo.cv_charge2 + controller.dIdt = 1.0 + + before_state = ( + Controller = controller, + ElectricPotential = [4.0], + Current = [-0.06], + ) + before_flags = BattMo.setupRegionSwitchFlags(policy, before_state, BattMo.cv_charge2) + @test before_flags.beforeSwitchRegion + @test !before_flags.afterSwitchRegion + + after_state = ( + Controller = controller, + ElectricPotential = [4.0], + Current = [-0.04], + ) + after_flags = BattMo.setupRegionSwitchFlags(policy, after_state, BattMo.cv_charge2) + @test !after_flags.beforeSwitchRegion + @test after_flags.afterSwitchRegion +end + +@testset "CC without ramp-up keeps one small initial timestep" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + delete!(model_settings, "RampUp") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 50.0 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "CC", + "InitialStateOfCharge" => 0.5, + "InitialControl" => "discharging", + "DRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 3.9, + "UpperVoltageLimit" => 4.1, + ) + ) + + sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) + + @test sim.time_steps[1] == simulation_settings["TimeStepDuration"] / 2 + @test sim.time_steps[2] == simulation_settings["TimeStepDuration"] +end + @testset "Crate" begin @test begin diff --git a/test/test_initial_state.jl b/test/test_initial_state.jl index 3d3df94bb..af0debce5 100644 --- a/test/test_initial_state.jl +++ b/test/test_initial_state.jl @@ -1,9 +1,6 @@ using BattMo using Test -trapz(x, y) = sum((y[1:(end - 1)] .+ y[2:end]) .* diff(x)) / 2 -rmse(x, y0, y1) = sqrt(trapz(x, (y1 .- y0) .^ 2) / (x[end] - x[1])) - @testset "Restart from saved state reproduces full discharge-charge simulation" begin # Parameters rate = 1.0 @@ -130,6 +127,6 @@ rmse(x, y0, y1) = sqrt(trapz(x, (y1 .- y0) .^ 2) / (x[end] - x[1])) @test t_merge ≈ t_ref atol = 1.0e-14 rtol = 0.0 @test I_merge ≈ output_ref.time_series["Current"] atol = 1.0e-7 rtol = 0.0 @test maximum(abs.(E_merge .- E_ref)) <= 1.0e-6 - @test rmse(t_ref, E_ref, E_merge) <= 1.0e-4 + @test BattMo.rmse(t_ref, E_ref, t_merge, E_merge) <= 1.0e-4 end end diff --git a/test/test_sei.jl b/test/test_sei.jl index 70f55f4a7..547eb8f55 100644 --- a/test/test_sei.jl +++ b/test/test_sei.jl @@ -14,6 +14,7 @@ using Test cell_parameters = load_cell_parameters(; from_file_path = file_path_cell) cycling_protocol = load_cycling_protocol(; from_file_path = file_path_cycling) model_settings = load_model_settings(; from_file_path = file_path_model) + delete!(model_settings, "RampUp") simulation_settings = load_simulation_settings(; from_file_path = file_path_simulation) @@ -35,7 +36,7 @@ using Test voltage_drop = states["NegativeElectrode"]["Interphase"]["VoltageDrop"] - @test length(sei_thickness[:, 2]) ≈ 2629 atol = 0 + @test length(sei_thickness[:, 2]) == 2629 @test sei_thickness[2, 2] ≈ 1.0000000339846753e-8 atol = 1.0e-1 @test voltage_drop[2, 2] ≈ -0.001401323958482127 atol = 1.0e-1 diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl new file mode 100644 index 000000000..72427f685 --- /dev/null +++ b/test/test_sequence_protocols.jl @@ -0,0 +1,394 @@ +using BattMo +using Jutul: report_timesteps +using Test + +function control_voltage(output) + return [only(state[:Control][:ElectricPotential]) for state in output.jutul_output.states] +end + +function control_current(output) + return [only(state[:Control][:Current]) for state in output.jutul_output.states] +end + +function sequence_step_index(output) + return [state[:Control][:Controller].step_index for state in output.jutul_output.states] +end + +function control_controller(output) + return [state[:Control][:Controller] for state in output.jutul_output.states] +end + +@testset "Rest protocol" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 50 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Rest", + "InitialStateOfCharge" => 0.5, + "Duration" => 100.0, + ) + ) + + output = solve(Simulation(model, cell_parameters, cycling_protocol; simulation_settings); info_level = -1) + + @test output.time_series["Time"][end] > 0.0 + @test all(abs.(output.time_series["Current"]) .< 1.0e-8) +end + +@testset "Sequence Rest steps use global controller time and local step start time" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + delete!(model_settings, "RampUp") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 10.0 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 0.5, + "Steps" => [ + Dict("Protocol" => "Rest", "Duration" => 20.0), + Dict("Protocol" => "Rest", "Duration" => 30.0), + Dict("Protocol" => "Rest", "Duration" => 40.0), + ], + ) + ) + + output = solve( + Simulation(model, cell_parameters, cycling_protocol; simulation_settings); + info_level = -1, + output_substates = true, + ) + + controllers = control_controller(output) + step_index = sequence_step_index(output) + report_time = cumsum(report_timesteps(output.jutul_output.reports; ministeps = true)) + + @test length(controllers) == length(report_time) + # When the final sequence step completes, the controller advances once more: + # step_index == length(Steps) + 1 means "sequence complete". + @test sort(unique(step_index)) == [1, 2, 3, 4] + + for k in 1:3 + indices = findall(==(k), step_index) + @test !isempty(indices) + if !isempty(indices) + step_start = sum(cycling_protocol["Steps"][i]["Duration"] for i in 1:(k - 1); init = 0.0) + @test all(controller -> isapprox(controller.step_start_time, step_start; atol = 1.0e-12, rtol = 0.0), controllers[indices]) + @test isapprox(controllers[indices[end]].time, report_time[indices[end]]; atol = 1.0e-12, rtol = 0.0) + @test controllers[indices[end]].time - controllers[indices[end]].step_start_time <= cycling_protocol["Steps"][k]["Duration"] + 1.0e-12 + end + end + + terminal = findall(==(4), step_index) + @test !isempty(terminal) + if !isempty(terminal) + sequence_duration = sum(step["Duration"] for step in cycling_protocol["Steps"]) + @test all(controller -> isapprox(controller.time, sequence_duration; atol = 1.0e-12, rtol = 0.0), controllers[terminal]) + end +end + +@testset "Sequence Rest advances when convergence is checked before solve" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 50.0 + simulation_settings["RampUpSteps"] = 5 + + solver_settings = load_solver_settings(; from_default_set = "direct") + solver_settings["NonLinearSolver"]["CheckBeforeSolve"] = true + solver_settings["NonLinearSolver"]["AlwaysUpdateSecondary"] = true + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 0.99, + "Steps" => [ + Dict( + "Protocol" => "CC", + "InitialControl" => "discharging", + "DRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 4.0, + "UpperVoltageLimit" => 4.2, + ), + Dict("Protocol" => "Rest", "Duration" => 600.0), + Dict( + "Protocol" => "CC", + "InitialControl" => "charging", + "CRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 3.0, + "UpperVoltageLimit" => 4.1, + ), + ], + ) + ) + + output = solve( + Simulation(model, cell_parameters, cycling_protocol; simulation_settings); + solver_settings, + info_level = -1, + output_substates = true, + ) + + controllers = control_controller(output) + step_index = sequence_step_index(output) + report_time = cumsum(report_timesteps(output.jutul_output.reports; ministeps = true)) + + # Motivation: during a rest step, CheckBeforeSolve can accept a timestep + # without solving because the system is already converged. The sequence + # controller must still update global time on that accepted timestep; + # otherwise the rest duration is never reached and later steps are skipped. + @test 3 in step_index + @test issorted([controller.time for controller in controllers]) + @test controllers[end].time > cycling_protocol["Steps"][2]["Duration"] + + rest = findall(==(2), step_index) + charge = findall(==(3), step_index) + @test !isempty(rest) + @test !isempty(charge) + if !isempty(rest) && !isempty(charge) + rest_start = controllers[rest[1]].step_start_time + rest_duration = cycling_protocol["Steps"][2]["Duration"] + # The final stored rest state can be the state immediately before the + # transition. The first charge state records the global time at which + # the rest step completed in step_start_time. + @test isapprox(controllers[charge[1]].step_start_time, rest_start + rest_duration; atol = 1.0e-12, rtol = 0.0) + end +end + +@testset "Sequence protocol with CC and Rest steps" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 10 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 0.99, + "Steps" => [ + Dict( + "Protocol" => "CC", + "InitialControl" => "discharging", + "DRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 4.0, + "UpperVoltageLimit" => 4.1, + ), + Dict( + "Protocol" => "Rest", + "Duration" => 20.0, + ), + ], + ) + ) + + output = solve(Simulation(model, cell_parameters, cycling_protocol; simulation_settings); info_level = -1) + current = output.time_series["Current"] + + @test any(current .> 0.1) + @test first(current) < maximum(current) + @test any(abs.(current) .< 1.0e-8) +end + +@testset "Sequence protocol advances through CC discharge Rest and CC charge" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + delete!(model_settings, "RampUp") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 20.0 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 1.0, + "Steps" => [ + Dict( + "Protocol" => "CC", + "InitialControl" => "discharging", + "DRate" => 1.0, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 3.0, + "UpperVoltageLimit" => 4.2, + ), + Dict( + "Protocol" => "Rest", + "Duration" => 120.0, + ), + Dict( + "Protocol" => "CC", + "InitialControl" => "charging", + "CRate" => 1.0, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 3.0, + "UpperVoltageLimit" => 3.8, + ), + ], + ) + ) + + output = solve( + Simulation(model, cell_parameters, cycling_protocol; simulation_settings); + info_level = -1, + output_substates = true, + ) + + voltage = control_voltage(output) + current = control_current(output) + step_index = sequence_step_index(output) + + discharge = findall(==(1), step_index) + rest = findall(==(2), step_index) + charge = findall(==(3), step_index) + + @test !isempty(discharge) + @test !isempty(rest) + @test !isempty(charge) + + if !isempty(discharge) + @test voltage[last(discharge)] < voltage[first(discharge)] + @test all(current[discharge] .> 0.0) + end + if length(rest) > 2 + rest_tail = rest[cld(length(rest), 2):end] + @test all(abs.(current[rest_tail]) .< 1.0e-8) + @test maximum(voltage[rest_tail]) - minimum(voltage[rest_tail]) < 0.05 + end + if !isempty(charge) + @test voltage[last(charge)] > voltage[first(charge)] + println("current during charge: ", current[charge]) + @test all(current[charge[2:end]] .< 0.0) + end +end + +@testset "Sequence protocol applies ramp-up at each current-controlled step" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 50.0 + simulation_settings["RampUpSteps"] = 3 + + model = LithiumIonBattery(; model_settings) + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 0.5, + "Steps" => [ + Dict( + "Protocol" => "CC", + "InitialControl" => "discharging", + "DRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 3.9, + "UpperVoltageLimit" => 4.1, + ), + Dict( + "Protocol" => "Rest", + "Duration" => 20.0, + ), + Dict( + "Protocol" => "CCCV", + "InitialControl" => "charging", + "CRate" => 1.0, + "DRate" => 1.0, + "TotalNumberOfCycles" => 1, + "LowerVoltageLimit" => 3.0, + "UpperVoltageLimit" => 4.0, + "CurrentChangeLimit" => 1.0e-4, + "VoltageChangeLimit" => 1.0e-4, + ), + ], + ) + ) + + sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) + sequence_policy = sim.model.multimodel[:Control].system.policy + cc_policy = BattMo.sequence_step_policy(sequence_policy.steps[1]) + rest_policy = BattMo.sequence_step_policy(sequence_policy.steps[2]) + cccv_policy = BattMo.sequence_step_policy(sequence_policy.steps[3]) + + @test !ismissing(cc_policy.current_function) + @test !ismissing(cccv_policy.current_function) + @test cc_policy.current_function(0.0) == 0.0 + @test cccv_policy.current_function(0.0) == 0.0 + @test cc_policy.current_function(simulation_settings["RampUpTime"]) ≈ cc_policy.ImaxDischarge + @test cccv_policy.current_function(simulation_settings["RampUpTime"]) ≈ cccv_policy.ImaxCharge + + ramp_timesteps = BattMo.compute_rampup_timesteps( + 1.1 * BattMo.Constants().hour / cycling_protocol["Steps"][1]["DRate"], + simulation_settings["TimeStepDuration"], + simulation_settings["RampUpSteps"], + ) + rest_ramp_timesteps = BattMo.compute_rampup_timesteps( + cycling_protocol["Steps"][2]["Duration"], + simulation_settings["TimeStepDuration"], + simulation_settings["RampUpSteps"], + ) + rest_index = length(ramp_timesteps) + 1 + next_step_index = rest_index + length(rest_ramp_timesteps) + + @test sim.time_steps[rest_index:(rest_index + simulation_settings["RampUpSteps"] - 1)] == rest_ramp_timesteps[1:simulation_settings["RampUpSteps"]] + @test sim.time_steps[next_step_index:(next_step_index + simulation_settings["RampUpSteps"] - 1)] == ramp_timesteps[1:simulation_settings["RampUpSteps"]] +end + +@testset "Sequence transition ramp is controlled by RampUp model setting" begin + cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") + model_settings = load_model_settings(; from_default_set = "p2d") + RampUp_setting = model_settings["RampUp"] + delete!(model_settings, "RampUp") + model_settings["Rampup"] = RampUp_setting + simulation_settings = load_simulation_settings(; from_default_set = "p2d") + simulation_settings["TimeStepDuration"] = 10.0 + simulation_settings["RampUpTime"] = 10.0 + simulation_settings["RampUpSteps"] = 3 + + cycling_protocol = CyclingProtocol( + Dict( + "Protocol" => "Sequence", + "InitialStateOfCharge" => 0.99, + "Steps" => [ + Dict( + "Protocol" => "CC", + "InitialControl" => "discharging", + "DRate" => 0.1, + "TotalNumberOfCycles" => 0, + "LowerVoltageLimit" => 4.0, + "UpperVoltageLimit" => 4.1, + ), + Dict( + "Protocol" => "Rest", + "Duration" => 20.0, + ), + ], + ) + ) + + model = LithiumIonBattery(; model_settings) + output = solve(Simulation(model, cell_parameters, cycling_protocol; simulation_settings); info_level = -1) + raw_states = output.jutul_output.states + step_index = [state[:Control][:Controller].step_index for state in raw_states] + transition_index = findfirst(i -> step_index[i - 1] == 1 && step_index[i] == 2, 2:length(step_index)) + current = output.time_series["Current"] + + @test transition_index !== nothing + ramp_window = current[transition_index:min(end, transition_index + simulation_settings["RampUpSteps"] + 1)] + @test any(0.0 .< ramp_window .< current[transition_index - 1]) + @test any(abs.(current[transition_index:end]) .< 1.0e-8) + + model_settings_without_ramp = load_model_settings(; from_default_set = "p2d") + delete!(model_settings_without_ramp, "RampUp") + model_without_ramp = LithiumIonBattery(; model_settings = model_settings_without_ramp) + sim_without_ramp = Simulation(model_without_ramp, cell_parameters, cycling_protocol; simulation_settings) + + @test !sim_without_ramp.model.multimodel[:Control].system.policy.use_ramp_up +end