From cc3d8a3235dc23b879e82ff49354ed9fbd65e71d Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 09:11:14 +0200 Subject: [PATCH 01/27] Starting rest and sequence controllers --- src/input/meta_data/cycling_protocol.jl | 15 +- src/input/schemas/get_schema.jl | 20 + .../schemas/lithium_ion/CyclingProtocol.json | 50 +- src/input/schemas/parameters_meta_data.json | 8 +- src/models/current_and_voltage_boundary.jl | 463 +++++++++++++++++- .../intercalation_battery.jl | 102 +++- src/output/output_format.jl | 2 +- src/simulation/simulation.jl | 34 +- test/runtests.jl | 1 + 9 files changed, 683 insertions(+), 12 deletions(-) diff --git a/src/input/meta_data/cycling_protocol.jl b/src/input/meta_data/cycling_protocol.jl index 288c83c7e..4aca99097 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, @@ -172,6 +172,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/schemas/get_schema.jl b/src/input/schemas/get_schema.jl index 025473bad..3e3866245 100644 --- a/src/input/schemas/get_schema.jl +++ b/src/input/schemas/get_schema.jl @@ -477,6 +477,8 @@ function get_schema_cycling_protocol(model_settings::ModelSettings) "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 +581,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"] diff --git a/src/input/schemas/lithium_ion/CyclingProtocol.json b/src/input/schemas/lithium_ion/CyclingProtocol.json index 06e7c82b1..33bbae64a 100644 --- a/src/input/schemas/lithium_ion/CyclingProtocol.json +++ b/src/input/schemas/lithium_ion/CyclingProtocol.json @@ -8,10 +8,13 @@ "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", @@ -136,6 +139,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 +293,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..523bd884d 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", @@ -614,4 +616,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..e21de8f70 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 # @@ -69,6 +71,22 @@ 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 +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 @@ -280,6 +298,15 @@ 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 # @@ -391,6 +418,41 @@ InputCurrentController() = InputCurrentController(0.0, 0.0, false) return R end +## RestController + +mutable struct RestController{R} <: Controller + target::R + time::R + target_is_voltage::Bool + ctrlType::OperationalMode +end + +RestController() = RestController(0.0, 0.0, false, rest) + +@inline function Jutul.numerical_type(x::RestController{R}) where {R} + return R +end + +## SequenceController + +mutable struct SequenceController{R} <: Controller + step_index::Int + step_start_time::R + target::R + time::R + target_is_voltage::Bool + ctrlType::Any + numberOfCycles::Int + dEdt::Union{R, Missing} + dIdt::Union{R, Missing} +end + +SequenceController() = SequenceController(1, 0.0, 0.0, 0.0, false, missing, 0, missing, missing) + +@inline function Jutul.numerical_type(x::SequenceController{R}) where {R} + return R +end + """ Function to create (deep) copy of simple controller @@ -509,6 +571,53 @@ function Base.copy(cv::InputCurrentController) return cv_copy +end + +""" +Function to create (deep) copy of rest controller +""" +function copyController!(cv_copy::RestController, cv::RestController) + + 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 Base.copy(cv::RestController) + + cv_copy = RestController() + copyController!(cv_copy, cv) + + return cv_copy + +end + +""" +Function to create (deep) copy of sequence controller +""" +function copyController!(cv_copy::SequenceController, cv::SequenceController) + + cv_copy.step_index = cv.step_index + cv_copy.step_start_time = cv.step_start_time + cv_copy.target = cv.target + cv_copy.time = cv.time + cv_copy.target_is_voltage = cv.target_is_voltage + cv_copy.ctrlType = cv.ctrlType + cv_copy.numberOfCycles = cv.numberOfCycles + cv_copy.dEdt = cv.dEdt + return cv_copy.dIdt = cv.dIdt + +end + +function Base.copy(cv::SequenceController) + + cv_copy = SequenceController() + copyController!(cv_copy, cv) + + return cv_copy + end """ We add the controller in the output @@ -586,12 +695,102 @@ function getInitCurrent(policy::InputCurrentPolicy) return policy.current_function(policy.times[1]) end +function getInitCurrent(policy::RestPolicy) + return zero(policy.duration) +end + +function getInitCurrent(policy::SequencePolicy) + return getInitCurrent(sequence_step_policy(policy.steps[1])) +end + function getInitCurrent(model::CurrentAndVoltageModel) return getInitCurrent(model.system.policy) end +sequence_step_policy(step::PolicyStep) = step.policy +sequence_step_policy(policy::AbstractPolicy) = policy + +function active_sequence_policy(policy::SequencePolicy, controller::SequenceController) + return sequence_step_policy(policy.steps[controller.step_index]) +end + +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 + +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 + +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 + +function sequence_step_complete(policy::RestPolicy, controller, state) + return controller.time - controller.step_start_time >= policy.duration +end + +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 + +function sequence_step_complete(policy::CyclingCVPolicy, controller, state) + return controller.numberOfCycles >= policy.numberOfCycles +end + +function sequence_complete(policy::SequencePolicy, controller::SequenceController) + return controller.step_index > length(policy.steps) +end + +function advance_sequence_step!(controller::SequenceController, policy::SequencePolicy) + 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)) + end + return controller +end + ###################################################### # Setup the initial policy from the input parameters # @@ -691,6 +890,14 @@ 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 + +function setup_initial_control_policy!(policy::SequencePolicy, input, parameters) + return nothing +end + ################################### # Special primary variable update # ################################### @@ -722,6 +929,28 @@ function Jutul.update_primary_variable!(state, p::CurrentVar, state_symbol, mode end +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 # ####################################### @@ -807,7 +1036,16 @@ 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 nextCtrlType = getNextCtrlTypecccv(ctrlType0) elseif policy isa CCPolicy if ctrlType == "discharging" @@ -869,6 +1107,18 @@ function Jutul.update_values!(old::InputCurrentController, new::InputCurrentCont end +function Jutul.update_values!(old::RestController, new::RestController) + + return copyController!(old, new) + +end + +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! """ @@ -926,6 +1176,36 @@ function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{Cu end +function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{RestPolicy{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 Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{SequencePolicy{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) @@ -1120,6 +1400,22 @@ 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 + +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 + ############################################################# # Functions to update the values in the controller in state # @@ -1375,6 +1671,93 @@ function update_values_in_controller!(state, policy::InputCurrentPolicy) end +function update_sequence_values_in_controller!(state, policy::RestPolicy) + return update_values_in_controller!(state, policy) +end + +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 + +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 + +function update_values_in_controller!(state, policy::SequencePolicy) + controller = state[:Controller] + if sequence_complete(policy, controller) + return nothing + end + return update_sequence_values_in_controller!(state, active_sequence_policy(policy, controller)) +end + ############################# # Assembly of the equations # ############################# @@ -1497,6 +1880,63 @@ function Jutul.update_after_step!(storage, domain::CurrentAndVoltageDomain, mode copyController!(storage.state0[:Controller], ctrl) + elseif policy isa RestPolicy + + copyController!(storage.state0[:Controller], ctrl) + + elseif policy isa SequencePolicy + + if sequence_complete(policy, ctrl) + return copyController!(storage.state0[:Controller], ctrl) + end + + step_policy = active_sequence_policy(policy, ctrl) + + if step_policy isa CyclingCVPolicy + initctrl = step_policy.initialControl + + ctrlType = ctrl.ctrlType + ctrlType0 = storage.state0[:Controller].ctrlType + ncycles = storage.state0[:Controller].numberOfCycles + + if initctrl == charging + if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) + ncycles = ncycles + 1 + end + elseif initctrl == discharging + if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) + ncycles = ncycles + 1 + end + end + + ctrl.numberOfCycles = ncycles + + elseif step_policy isa CCPolicy && step_policy.numberOfCycles > 0 + initctrl = step_policy.initialControl + + ctrlType = ctrl.ctrlType + ctrlType0 = storage.state0[:Controller].ctrlType + ncycles = storage.state0[:Controller].numberOfCycles + + if initctrl == "charging" + if ctrlType0 == "discharging" && ctrlType == "charging" + ncycles = ncycles + 1 + end + elseif initctrl == "discharging" + if ctrlType0 == "charging" && ctrlType == "discharging" + ncycles = ncycles + 1 + end + end + + ctrl.numberOfCycles = ncycles + end + + if sequence_step_complete(step_policy, ctrl, storage.state) + advance_sequence_step!(ctrl, policy) + end + + copyController!(storage.state0[:Controller], ctrl) + else error("Policy $(typeof(policy)) not recognized") @@ -1596,6 +2036,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) + initialize_sequence_step_controller!(state[:Controller], active_sequence_policy(policy, state[:Controller])) + end end diff --git a/src/models/full_battery_models/intercalation_battery.jl b/src/models/full_battery_models/intercalation_battery.jl index 7e6437263..31f7eb105 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( @@ -263,6 +265,80 @@ 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) + + 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, + ) + 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 + + policy = SequencePolicy(steps, ImaxDischarge, ImaxCharge) + else error("controlPolicy not recognized.") @@ -822,6 +898,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 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 28af39704..457e45315 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -797,13 +797,15 @@ function solver_configuration( end end - if model[:Control].system.policy isa CyclingCVPolicy || model[:Control].system.policy isa CCPolicy + if model[:Control].system.policy isa CyclingCVPolicy || model[:Control].system.policy isa CCPolicy || model[:Control].system.policy isa RestPolicy || model[:Control].system.policy isa SequencePolicy if model[:Control].system.policy isa CyclingCVPolicy cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) elseif model[:Control].system.policy isa CCPolicy && model[:Control].system.policy.numberOfCycles > 0 cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) + elseif model[: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) @@ -848,6 +850,11 @@ function solver_configuration( report[:stopnow] = false end end + elseif model[:Control].system.policy isa RestPolicy + report[:stopnow] = s.state.Control.Controller.time >= m[:Control].system.policy.duration + + elseif model[:Control].system.policy isa SequencePolicy + report[:stopnow] = sequence_complete(m[:Control].system.policy, s.state.Control.Controller) end return (done, report) @@ -1077,6 +1084,31 @@ function setup_timesteps( timesteps = diff(series_times) timesteps = timesteps[timesteps .> 0.0] + elseif protocol == "Rest" + + totalTime = cycling_protocol["Duration"] + dt = simulation_settings["TimeStepDuration"] + n = Int64(floor(totalTime / dt)) + timesteps = repeat([dt], n) + dt_final = totalTime - sum(timesteps) + if dt_final > 0 + push!(timesteps, dt_final) + 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") 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") From 9758aafa9aab0802469a2798de39df59a1502a9b Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 09:15:56 +0200 Subject: [PATCH 02/27] Test sequence protocol --- test/test_sequence_protocols.jl | 104 ++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 test/test_sequence_protocols.jl diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl new file mode 100644 index 000000000..4f928bb81 --- /dev/null +++ b/test/test_sequence_protocols.jl @@ -0,0 +1,104 @@ +using BattMo +using Test + +@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 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 accepts combined CC Rest and CCCV 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"] = 50 + + 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) + + @test !isempty(sim.time_steps) + @test length(cycling_protocol["Steps"]) == 3 +end From 9fb5251998f0a62f10b758c29acc9322407a3405 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 12:37:33 +0200 Subject: [PATCH 03/27] Add rampup --- src/simulation/simulation.jl | 9 ++++- test/test_sequence_protocols.jl | 65 +++++++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+), 2 deletions(-) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 457e45315..0042c0643 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -1066,9 +1066,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" diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl index 4f928bb81..d4aa04a7b 100644 --- a/test/test_sequence_protocols.jl +++ b/test/test_sequence_protocols.jl @@ -102,3 +102,68 @@ end @test !isempty(sim.time_steps) @test length(cycling_protocol["Steps"]) == 3 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_index = length(ramp_timesteps) + 1 + next_step_index = rest_index + 1 + + @test rest_policy.duration == sim.time_steps[rest_index] + @test sim.time_steps[next_step_index:(next_step_index + simulation_settings["RampUpSteps"] - 1)] == ramp_timesteps[1:simulation_settings["RampUpSteps"]] +end From ddab960dda5764f78061f62f185c9e567f90c490 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 12:37:46 +0200 Subject: [PATCH 04/27] example demonstrating sequence control --- examples/example_sequence_control.jl | 75 ++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 examples/example_sequence_control.jl diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl new file mode 100644 index 000000000..051d7d575 --- /dev/null +++ b/examples/example_sequence_control.jl @@ -0,0 +1,75 @@ +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. +### + +model_settings = load_model_settings(; from_default_set = "p2d") +#delete!(model_settings, "Rampup") +model = LithiumIonBattery(; model_settings = model_settings) +cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") +simulation_settings = load_simulation_settings(; from_default_set = "p2d") + +simulation_settings["TimeStepDuration"] = 500.0 +# simulation_settings["RampUpTime"] = 10.0 +# simulation_settings["RampUpSteps"] = 5 +# delete!(simulation_settings, "RampUpTime") +# delete!(simulation_settings, "RampUpSteps") + +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] + ) +) + +sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) +output = solve(sim; info_level = -1) + +time = output.time_series["Time"] +current = output.time_series["Current"] +voltage = output.time_series["Voltage"] + +println("Completed sequence-control example with $(length(cycling_protocol["Steps"])) steps.") +println("Ramp-up time: $(simulation_settings["RampUpTime"]) s") +println("Ramp-up steps: $(simulation_settings["RampUpSteps"])") +println("Final time: $(round(time[end]; digits = 2)) s") +println("Current range: $(round(minimum(current); digits = 4)) A to $(round(maximum(current); digits = 4)) A") +println("Voltage range: $(round(minimum(voltage); digits = 4)) V to $(round(maximum(voltage); digits = 4)) V") + +doplot = true +if doplot + using GLMakie + plot_dashboard(output; plot_type = "simple") +end From 8263615742dce29973fc595cec5edac55fc79e92 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 12:37:58 +0200 Subject: [PATCH 05/27] rename vars --- src/simulation/simulation.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 0042c0643..de52316e8 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -1021,13 +1021,13 @@ 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"] From 6a502b5f62ba5a9e6864dcc5317c21900c286583 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 20:49:17 +0200 Subject: [PATCH 06/27] Allow for rampup at each discontinuity and new sequence step --- src/models/current_and_voltage_boundary.jl | 91 ++++++++++++++++++- .../intercalation_battery.jl | 3 +- src/simulation/simulation.jl | 31 ++++--- test/test_sequence_protocols.jl | 61 ++++++++++++- 4 files changed, 167 insertions(+), 19 deletions(-) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index e21de8f70..7e31028d5 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -85,6 +85,8 @@ mutable struct SequencePolicy{R} <: AbstractPolicy steps::Vector{AbstractSequenceStep} ImaxDischarge::R ImaxCharge::R + use_ramp_up::Bool + rampup_time::R end ## A policy is used to compute the next control from the current control and state @@ -445,9 +447,15 @@ mutable struct SequenceController{R} <: Controller 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 -SequenceController() = SequenceController(1, 0.0, 0.0, 0.0, false, missing, 0, missing, missing) +SequenceController() = SequenceController(1, 0.0, 0.0, 0.0, false, missing, 0, missing, missing, 0.0, 0.0, 0.0, 0.0, false, false) @inline function Jutul.numerical_type(x::SequenceController{R}) where {R} return R @@ -607,7 +615,13 @@ function copyController!(cv_copy::SequenceController, cv::SequenceController) cv_copy.ctrlType = cv.ctrlType cv_copy.numberOfCycles = cv.numberOfCycles cv_copy.dEdt = cv.dEdt - return cv_copy.dIdt = cv.dIdt + cv_copy.dIdt = cv.dIdt + cv_copy.ramp_start_time = cv.ramp_start_time + cv_copy.ramp_duration = cv.ramp_duration + cv_copy.ramp_start_target = cv.ramp_start_target + cv_copy.ramp_end_target = cv.ramp_end_target + cv_copy.ramp_target_is_voltage = cv.ramp_target_is_voltage + return cv_copy.ramp_active = cv.ramp_active end @@ -782,11 +796,79 @@ function sequence_complete(policy::SequencePolicy, controller::SequenceControlle return controller.step_index > length(policy.steps) end +function sequence_nominal_target(policy::RestPolicy, controller::SequenceController) + return (target = 0.0, target_is_voltage = false) +end + +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 + +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 + 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 + +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 @@ -1755,7 +1837,8 @@ function update_values_in_controller!(state, policy::SequencePolicy) if sequence_complete(policy, controller) return nothing end - return update_sequence_values_in_controller!(state, active_sequence_policy(policy, controller)) + update_sequence_values_in_controller!(state, active_sequence_policy(policy, controller)) + return apply_sequence_transition_ramp!(controller) end ############################# @@ -2052,7 +2135,7 @@ function Jutul.initialize_extra_state_fields!(state, ::Any, model::CurrentAndVol target_is_voltage = false target, time = promote(target, time) - state[:Controller] = SequenceController(1, time, target, time, target_is_voltage, missing, 0, missing, missing) + 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 diff --git a/src/models/full_battery_models/intercalation_battery.jl b/src/models/full_battery_models/intercalation_battery.jl index 31f7eb105..276867b3f 100644 --- a/src/models/full_battery_models/intercalation_battery.jl +++ b/src/models/full_battery_models/intercalation_battery.jl @@ -337,7 +337,8 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) end end - policy = SequencePolicy(steps, ImaxDischarge, ImaxCharge) + rampup_time = use_ramp_up ? Float64(input.simulation_settings["RampUpTime"]) : zero(ImaxDischarge) + policy = SequencePolicy(steps, ImaxDischarge, ImaxCharge, use_ramp_up, rampup_time) else diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index de52316e8..5cec45f9d 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -1093,24 +1093,31 @@ function setup_timesteps( totalTime = cycling_protocol["Duration"] dt = simulation_settings["TimeStepDuration"] - n = Int64(floor(totalTime / dt)) - timesteps = repeat([dt], n) - dt_final = totalTime - sum(timesteps) - if dt_final > 0 - push!(timesteps, dt_final) + 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, - )) + 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 diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl index d4aa04a7b..f462ac1ae 100644 --- a/test/test_sequence_protocols.jl +++ b/test/test_sequence_protocols.jl @@ -161,9 +161,66 @@ end 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 + 1 + next_step_index = rest_index + length(rest_ramp_timesteps) - @test rest_policy.duration == sim.time_steps[rest_index] + @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 From 3f9e61e0ddad07342d606b4b313625b083846371 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Tue, 16 Jun 2026 22:40:43 +0200 Subject: [PATCH 07/27] Work on cccv rampup --- examples/example_sequence_control.jl | 76 ++++++++--- src/models/current_and_voltage_boundary.jl | 124 ++++++++++++++++-- .../intercalation_battery.jl | 1 + src/simulation/simulation.jl | 47 ++++++- test/test_cycling_protocols.jl | 49 +++++++ 5 files changed, 265 insertions(+), 32 deletions(-) diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl index 051d7d575..9de7473ca 100644 --- a/examples/example_sequence_control.jl +++ b/examples/example_sequence_control.jl @@ -1,5 +1,5 @@ using BattMo -using Jutul: si_unit +using Jutul: expand_to_ministeps, si_unit ### # Example: Sequence control with CC, Rest, and CCCV steps @@ -11,16 +11,14 @@ using Jutul: si_unit ### model_settings = load_model_settings(; from_default_set = "p2d") -#delete!(model_settings, "Rampup") + +# Optionally disable ramp-up to see the effect of the initial ramp-up steps on the current and voltage profiles. +# delete!(model_settings, "RampUp") + model = LithiumIonBattery(; model_settings = model_settings) cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") simulation_settings = load_simulation_settings(; from_default_set = "p2d") - -simulation_settings["TimeStepDuration"] = 500.0 -# simulation_settings["RampUpTime"] = 10.0 -# simulation_settings["RampUpSteps"] = 5 -# delete!(simulation_settings, "RampUpTime") -# delete!(simulation_settings, "RampUpSteps") +simulation_settings["TimeStepDuration"] = 360.0 cc = Dict( "Protocol" => "CC", @@ -32,7 +30,7 @@ cc = Dict( ) rest = Dict( "Protocol" => "Rest", - "Duration" => 1.0*si_unit("hour") + "Duration" => 1.0 * si_unit("hour") ) cccv = Dict( "Protocol" => "CCCV", @@ -55,21 +53,57 @@ cycling_protocol = CyclingProtocol( ) sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) -output = solve(sim; info_level = -1) +output = solve(sim; info_level = -1, output_substates = true) -time = output.time_series["Time"] -current = output.time_series["Current"] -voltage = output.time_series["Voltage"] - -println("Completed sequence-control example with $(length(cycling_protocol["Steps"])) steps.") -println("Ramp-up time: $(simulation_settings["RampUpTime"]) s") -println("Ramp-up steps: $(simulation_settings["RampUpSteps"])") -println("Final time: $(round(time[end]; digits = 2)) s") -println("Current range: $(round(minimum(current); digits = 4)) A to $(round(maximum(current); digits = 4)) A") -println("Voltage range: $(round(minimum(voltage); digits = 4)) V to $(round(maximum(voltage); digits = 4)) V") +stored_states = output.jutul_output.states +raw_states, dt, = expand_to_ministeps(stored_states, output.jutul_output.reports[eachindex(stored_states)]) +time = cumsum(dt) +current = [only(state[:Control][:Current]) for state in raw_states] +voltage = [only(state[:Control][:ElectricPotential]) for state in raw_states] doplot = true if doplot using GLMakie - plot_dashboard(output; plot_type = "simple") + + sequence_step_index = [state[:Control][:Controller].step_index for state in raw_states] + sequence_time = [state[:Control][:Controller].time for state in raw_states] + step_boundary_times = Float64[] + + for i in 2:length(sequence_step_index) + if sequence_step_index[i] != sequence_step_index[i - 1] && sequence_step_index[i] <= length(cycling_protocol["Steps"]) + push!(step_boundary_times, sequence_time[i]) + end + end + + hour = si_unit("hour") + time_hours = time ./ hour + dt_minutes = dt ./ si_unit("minute") + timestep_time_hours = Float64[] + timestep_values = Float64[] + for i in eachindex(dt) + push!(timestep_time_hours, (time[i] - dt[i]) / hour) + push!(timestep_time_hours, time_hours[i]) + push!(timestep_values, dt_minutes[i]) + push!(timestep_values, dt_minutes[i]) + end + boundary_hours = step_boundary_times ./ hour + + fig = Figure(size = (900, 760)) + ax_voltage = Axis(fig[1, 1], xlabel = "Time / h", ylabel = "Voltage / V") + ax_current = Axis(fig[2, 1], xlabel = "Time / h", ylabel = "Current / A") + ax_timestep = Axis(fig[3, 1], xlabel = "Time / h", ylabel = "Timestep / min") + linkxaxes!(ax_voltage, ax_current, ax_timestep) + + scatterlines!(ax_voltage, time_hours, voltage, label = "Voltage", linewidth = 4.0, markersize = 10, markercolor = :black) + scatterlines!(ax_current, time_hours, current, label = "Current", linewidth = 4.0, markersize = 10, markercolor = :black) + lines!(ax_timestep, timestep_time_hours, timestep_values, label = "Timestep", linewidth = 3.0) + + if !isempty(boundary_hours) + vlines!(ax_voltage, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) + vlines!(ax_current, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) + vlines!(ax_timestep, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) + end + + screen = GLMakie.Screen() + GLMakie.display(screen, fig) end diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 7e31028d5..07329ea05 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -200,6 +200,7 @@ mutable struct CyclingCVPolicy{R, I} <: AbstractPolicy numberOfCycles::I tolerances::Any use_ramp_up::Bool + rampup_time::R current_function::Any end @@ -214,6 +215,7 @@ function CyclingCVPolicy( ImaxDischarge = 0 * lowerCutoffVoltage, ImaxCharge = 0 * lowerCutoffVoltage, use_ramp_up::Bool, + rampup_time = zero(lowerCutoffVoltage), current_function = missing, ) @@ -243,6 +245,7 @@ function CyclingCVPolicy( numberOfCycles, tolerances, use_ramp_up, + rampup_time, current_function, ) end @@ -310,9 +313,9 @@ function Jutul.select_parameters!( 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 @@ -360,6 +363,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 + ramp_end_target::R + ramp_target_is_voltage::Bool + ramp_active::Bool end @@ -367,7 +376,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 @@ -462,6 +471,10 @@ SequenceController() = SequenceController(1, 0.0, 0.0, 0.0, false, missing, 0, m end +############################################## +# Copy helpers for the different controllers # +############################################## + """ Function to create (deep) copy of simple controller """ @@ -505,7 +518,15 @@ 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 + cv_copy.numberOfCycles = cv.numberOfCycles + cv_copy.dEdt = cv.dEdt + cv_copy.dIdt = cv.dIdt + cv_copy.ramp_start_time = cv.ramp_start_time + cv_copy.ramp_duration = cv.ramp_duration + cv_copy.ramp_start_target = cv.ramp_start_target + cv_copy.ramp_end_target = cv.ramp_end_target + cv_copy.ramp_target_is_voltage = cv.ramp_target_is_voltage + return cv_copy.ramp_active = cv.ramp_active end @@ -633,6 +654,11 @@ function Base.copy(cv::SequenceController) return cv_copy end + +################################################# +# Minimum output variables for the control model # +################################################# + """ We add the controller in the output """ @@ -726,6 +752,10 @@ end sequence_step_policy(step::PolicyStep) = step.policy sequence_step_policy(policy::AbstractPolicy) = policy +##################################### +# Helpers for sequence policy steps # +##################################### + function active_sequence_policy(policy::SequencePolicy, controller::SequenceController) return sequence_step_policy(policy.steps[controller.step_index]) end @@ -957,6 +987,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 @@ -1128,6 +1159,9 @@ function check_constraints(model, storage) 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" @@ -1331,9 +1365,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])) @@ -1349,6 +1438,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) @@ -1400,6 +1493,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 @@ -1499,9 +1596,9 @@ function update_control_type_in_controller!(state, state0, policy::SequencePolic end -############################################################# -# Functions to update the values in the controller in state # -############################################################# +################################################## +# Control-type update for input-current policies # +################################################## """ Implementation of the InputCurrentPolicy. @@ -1560,6 +1657,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 @@ -1726,7 +1827,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 @@ -1753,6 +1855,10 @@ function update_values_in_controller!(state, policy::InputCurrentPolicy) end +############################################# +# Target-value updates for sequence steps # +############################################# + function update_sequence_values_in_controller!(state, policy::RestPolicy) return update_values_in_controller!(state, policy) end diff --git a/src/models/full_battery_models/intercalation_battery.jl b/src/models/full_battery_models/intercalation_battery.jl index 276867b3f..7a323ad37 100644 --- a/src/models/full_battery_models/intercalation_battery.jl +++ b/src/models/full_battery_models/intercalation_battery.jl @@ -321,6 +321,7 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) ImaxDischarge = Idis, ImaxCharge = Ichg, use_ramp_up = use_ramp_up, + rampup_time = use_ramp_up ? Float64(input.simulation_settings["RampUpTime"]) : zero(Idis), ) if use_ramp_up Imax = step["InitialControl"] == "charging" ? Ichg : Idis diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 5cec45f9d..b7699f89e 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. """ @@ -430,6 +457,8 @@ function solve_simulation( model.multimodel, parameters; solver_settings, + rampup_steps = get(simulation_settings, "RampUpSteps", nothing), + rampup_reference_timestep = get(simulation_settings, "TimeStepDuration", nothing), validate, accept_invalid, logger, @@ -679,6 +708,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, @@ -728,11 +759,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(model[: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 @@ -774,6 +811,12 @@ function solver_configuration( output_substates = output["OutputSubstrates"], ) + 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 diff --git a/test/test_cycling_protocols.jl b/test/test_cycling_protocols.jl index 515b3c53e..399b30de7 100644 --- a/test/test_cycling_protocols.jl +++ b/test/test_cycling_protocols.jl @@ -1,6 +1,55 @@ 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 "Crate" begin @test begin From 3ed3c24fcd0433511b1a3b86ed8cfb952e4bcb8d Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 09:02:40 +0200 Subject: [PATCH 08/27] move current function --- src/models/current_and_voltage_boundary.jl | 17 ++++++++++++++--- src/simulation/simulation.jl | 18 ------------------ 2 files changed, 14 insertions(+), 21 deletions(-) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 07329ea05..752d59893 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -2284,9 +2284,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/simulation/simulation.jl b/src/simulation/simulation.jl index b7699f89e..3a75f87d6 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -1203,21 +1203,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 From 502715723001051a39d36bd60b7cc0452a783427 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 09:03:49 +0200 Subject: [PATCH 09/27] clean up and add comments --- src/models/current_and_voltage_boundary.jl | 523 +++++++-------------- 1 file changed, 167 insertions(+), 356 deletions(-) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 752d59893..8e594ee82 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -18,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 ############################################# @@ -89,6 +87,37 @@ mutable struct SequencePolicy{R} <: AbstractPolicy 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 @@ -135,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 @@ -148,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 @@ -399,7 +424,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 @@ -464,8 +489,10 @@ mutable struct SequenceController{R} <: Controller 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) +# Return the numeric scalar type used by the sequence controller. @inline function Jutul.numerical_type(x::SequenceController{R}) where {R} return R end @@ -475,184 +502,64 @@ end # Copy helpers for the different controllers # ############################################## -""" -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 - -end - -""" -Function to create (deep) copy of CC-CV controller -""" -function copyController!(cv_copy::CcCvController, cv::CcCvController) - - copyController!(cv_copy.maincontroller, cv.maincontroller) - cv_copy.numberOfCycles = cv.numberOfCycles - cv_copy.dEdt = cv.dEdt - cv_copy.dIdt = cv.dIdt - cv_copy.ramp_start_time = cv.ramp_start_time - cv_copy.ramp_duration = cv.ramp_duration - cv_copy.ramp_start_target = cv.ramp_start_target - cv_copy.ramp_end_target = cv.ramp_end_target - cv_copy.ramp_target_is_voltage = cv.ramp_target_is_voltage - return cv_copy.ramp_active = cv.ramp_active - -end - -""" -Overload function to copy simple controller -""" -function Base.copy(cv::FunctionController) - - cv_copy = FunctionController() - copyController!(cv_copy, cv) - - return cv_copy - -end - - -""" -Overload function to copy simple controller -""" -function Base.copy(cv::CCController) - - cv_copy = CCController() - copyController!(cv_copy, cv) - - return cv_copy - -end - -""" -Overload function to copy simple controller -""" -function Base.copy(cv::SimpleControllerCV) - - cv_copy = SimpleControllerCV() - copyController!(cv_copy, cv) - - return cv_copy - -end - -""" -Overload function to copy CC-CV controller -""" -function Base.copy(cv::CcCvController) - - cv_copy = CcCvController() - copyController!(cv_copy, cv) - - return cv_copy - -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 - -end - -""" -Overload function to copy input current controller -""" -function Base.copy(cv::InputCurrentController) - - cv_copy = InputCurrentController() - copyController!(cv_copy, cv) - - return cv_copy - -end - -""" -Function to create (deep) copy of rest controller -""" -function copyController!(cv_copy::RestController, cv::RestController) - - 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 Base.copy(cv::RestController) - - cv_copy = RestController() - copyController!(cv_copy, cv) - - return cv_copy - +function copyController!(dst::T, src::T) where {T <: Controller} + for field in fieldnames(T) + setfield!(dst, field, getfield(src, field)) + end + return dst end -""" -Function to create (deep) copy of sequence controller -""" -function copyController!(cv_copy::SequenceController, cv::SequenceController) - - cv_copy.step_index = cv.step_index - cv_copy.step_start_time = cv.step_start_time - cv_copy.target = cv.target - cv_copy.time = cv.time - cv_copy.target_is_voltage = cv.target_is_voltage - cv_copy.ctrlType = cv.ctrlType - cv_copy.numberOfCycles = cv.numberOfCycles - cv_copy.dEdt = cv.dEdt - cv_copy.dIdt = cv.dIdt - cv_copy.ramp_start_time = cv.ramp_start_time - cv_copy.ramp_duration = cv.ramp_duration - cv_copy.ramp_start_target = cv.ramp_start_target - cv_copy.ramp_end_target = cv.ramp_end_target - cv_copy.ramp_target_is_voltage = cv.ramp_target_is_voltage - return cv_copy.ramp_active = cv.ramp_active - +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 Base.copy(cv::SequenceController) - - cv_copy = SequenceController() - copyController!(cv_copy, cv) - - return cv_copy - +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 ################################################# @@ -739,6 +646,7 @@ 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 @@ -749,17 +657,21 @@ function getInitCurrent(model::CurrentAndVoltageModel) 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 @@ -769,6 +681,7 @@ function initialize_sequence_step_controller!(controller::SequenceController, st 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" @@ -785,6 +698,7 @@ function initialize_sequence_step_controller!(controller::SequenceController, st 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 @@ -800,10 +714,12 @@ function initialize_sequence_step_controller!(controller::SequenceController, st 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" @@ -818,18 +734,22 @@ function sequence_step_complete(policy::CCPolicy, controller, state) 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 @@ -841,6 +761,7 @@ function sequence_nominal_target(policy::CCPolicy, controller::SequenceControlle 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 @@ -860,6 +781,7 @@ function sequence_nominal_target(policy::CyclingCVPolicy, controller::SequenceCo 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 @@ -887,6 +809,7 @@ function advance_sequence_step!(controller::SequenceController, policy::Sequence 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 @@ -908,6 +831,8 @@ 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 @@ -935,7 +860,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]) @@ -944,10 +870,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 @@ -1007,6 +929,7 @@ 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 @@ -1042,6 +965,7 @@ 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) @@ -1088,10 +1012,7 @@ function setupRegionSwitchFlags(policy::Union{CyclingCVPolicy, CCPolicy}, state, end - E = only(state.ElectricPotential) - I = only(state.Current) - if ctrlType == cc_discharge1 || ctrlType == "discharging" @@ -1229,6 +1150,7 @@ function Jutul.update_values!(old::RestController, new::RestController) end +# Copy the sequence controller when Jutul updates controller values. function Jutul.update_values!(old::SequenceController, new::SequenceController) return copyController!(old, new) @@ -1249,35 +1171,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} - - invoke( - reset_state_to_previous_state!, - Tuple{ - typeof(storage), - SimulationModel, - }, - storage, - model, - ) - 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} + model::SimulationModel{CurrentAndVoltageDomain, <:CurrentAndVoltageSystem, T3, T4}, + ) where {T3, T4} invoke( reset_state_to_previous_state!, @@ -1289,37 +1186,6 @@ function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{Cu model, ) return copyController!(storage.state[:Controller], storage.state0[:Controller]) - -end - -function Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{RestPolicy{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 Jutul.reset_state_to_previous_state!(storage, model::SimulationModel{CurrentAndVoltageDomain, CurrentAndVoltageSystem{SequencePolicy{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 @@ -1586,6 +1452,7 @@ function update_control_type_in_controller!(state, state0, policy::RestPolicy, d return controller.ctrlType = rest end +# 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) @@ -1859,10 +1726,12 @@ 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 @@ -1893,6 +1762,7 @@ function update_sequence_values_in_controller!(state, policy::CCPolicy) 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 @@ -1938,6 +1808,7 @@ function update_values_in_controller!(state, policy::RestPolicy) 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) @@ -1976,7 +1847,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) @@ -1989,149 +1859,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 - +function update_cycle_count!(ctrl, ctrl0, initialControl::OperationalMode) ctrlType = ctrl.ctrlType + ctrlType0 = ctrl0.ctrlType + ncycles = ctrl0.numberOfCycles - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles - - copyController!(storage.state0[:Controller], ctrl) - - if initctrl == charging - + if initialControl == charging if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) - ncycles = ncycles + 1 + ncycles += 1 end - - elseif initctrl == discharging - + elseif initialControl == discharging if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) - ncycles = ncycles + 1 + ncycles += 1 end - + else + error("Initial control $initialControl is not recognized") end ctrl.numberOfCycles = ncycles + return ctrl +end - elseif policy isa SimpleCVPolicy - - copyController!(storage.state0[:Controller], ctrl) - - elseif policy isa FunctionPolicy - - copyController!(storage.state0[:Controller], ctrl) - - elseif policy isa CCPolicy - - if policy.numberOfCycles == 0 - copyController!(storage.state0[:Controller], ctrl) - - else - - initctrl = policy.initialControl - +function update_cycle_count!(ctrl, ctrl0, initialControl::String) ctrlType = ctrl.ctrlType + ctrlType0 = ctrl0.ctrlType + ncycles = ctrl0.numberOfCycles - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles - - copyController!(storage.state0[:Controller], ctrl) - - if initctrl == "charging" - + if initialControl == "charging" if ctrlType0 == "discharging" && ctrlType == "charging" - ncycles = ncycles + 1 + ncycles += 1 end - - elseif initctrl == "discharging" - + elseif initialControl == "discharging" if ctrlType0 == "charging" && ctrlType == "discharging" - ncycles = ncycles + 1 + ncycles += 1 end - + else + error("Initial control $initialControl is not recognized") end ctrl.numberOfCycles = ncycles + return ctrl end - elseif policy isa InputCurrentPolicy +""" 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) + + ctrl = storage.state[:Controller] + ctrl0 = storage.state0[:Controller] + policy = model.system.policy - copyController!(storage.state0[:Controller], ctrl) + if policy isa CyclingCVPolicy + update_cycle_count!(ctrl, ctrl0, policy.initialControl) + return copyController!(ctrl0, ctrl) - elseif policy isa RestPolicy + elseif policy isa SimpleCVPolicy || policy isa FunctionPolicy || policy isa InputCurrentPolicy || policy isa RestPolicy + return copyController!(ctrl0, ctrl) - copyController!(storage.state0[:Controller], ctrl) + elseif policy isa CCPolicy + if policy.numberOfCycles > 0 + update_cycle_count!(ctrl, ctrl0, policy.initialControl) + end + return copyController!(ctrl0, ctrl) elseif policy isa SequencePolicy - if sequence_complete(policy, ctrl) - return copyController!(storage.state0[:Controller], ctrl) + return copyController!(ctrl0, ctrl) end step_policy = active_sequence_policy(policy, ctrl) if step_policy isa CyclingCVPolicy - initctrl = step_policy.initialControl - - ctrlType = ctrl.ctrlType - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles - - if initctrl == charging - if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) - ncycles = ncycles + 1 - end - elseif initctrl == discharging - if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) - ncycles = ncycles + 1 - end - end - - ctrl.numberOfCycles = ncycles - + update_cycle_count!(ctrl, ctrl0, step_policy.initialControl) elseif step_policy isa CCPolicy && step_policy.numberOfCycles > 0 - initctrl = step_policy.initialControl - - ctrlType = ctrl.ctrlType - ctrlType0 = storage.state0[:Controller].ctrlType - ncycles = storage.state0[:Controller].numberOfCycles - - if initctrl == "charging" - if ctrlType0 == "discharging" && ctrlType == "charging" - ncycles = ncycles + 1 - end - elseif initctrl == "discharging" - if ctrlType0 == "charging" && ctrlType == "discharging" - ncycles = ncycles + 1 - end - end - - ctrl.numberOfCycles = ncycles + update_cycle_count!(ctrl, ctrl0, step_policy.initialControl) end 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 ######################################################################## From c2fa3ae6f050675a933c4eaf66fc28fad620ae8c Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 10:42:31 +0200 Subject: [PATCH 10/27] RampUp documentation --- src/input/meta_data/model_settings.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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".""", ), ) From 87190a4b8023847038c36772713235ffe92a8785 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 10:42:51 +0200 Subject: [PATCH 11/27] Disable rampup for the tests to pass since we now return extra timesteps --- test/test_capacity_calculation.jl | 4 +++- test/test_sei.jl | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) 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_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 From 212620a4ae0d921e49b6af4b9104518ca9ab1f9d Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 13:18:30 +0200 Subject: [PATCH 12/27] Add comments --- src/models/current_and_voltage_boundary.jl | 41 +++++++++++----------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 8e594ee82..5dea07c0e 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -390,10 +390,10 @@ mutable struct CcCvController{R, I <: Integer} <: Controller 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 + 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 @@ -472,6 +472,7 @@ end ## 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 @@ -529,7 +530,7 @@ function copyController!(dst::T, src::S) where {T <: Controller, S <: Controller end function copyController!(dst::CcCvController, src::CcCvController) -# Need its own copyController since it contains a nested mutable controller + # Need its own copyController since it contains a nested mutable controller copyController!(dst.maincontroller, src.maincontroller) dst.numberOfCycles = src.numberOfCycles dst.dEdt = src.dEdt @@ -1860,46 +1861,46 @@ end ##################################################################### function update_cycle_count!(ctrl, ctrl0, initialControl::OperationalMode) - ctrlType = ctrl.ctrlType + ctrlType = ctrl.ctrlType ctrlType0 = ctrl0.ctrlType ncycles = ctrl0.numberOfCycles if initialControl == charging - if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) + if (ctrlType0 == cc_discharge1 || ctrlType0 == cc_discharge2) && (ctrlType == cc_charge1 || ctrlType == cv_charge2) ncycles += 1 - end + end elseif initialControl == discharging - if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) + if (ctrlType0 == cc_charge1 || ctrlType0 == cv_charge2) && (ctrlType == cc_discharge1 || ctrlType == cc_discharge2) ncycles += 1 - end + end else error("Initial control $initialControl is not recognized") - end + end - ctrl.numberOfCycles = ncycles + ctrl.numberOfCycles = ncycles return ctrl end function update_cycle_count!(ctrl, ctrl0, initialControl::String) - ctrlType = ctrl.ctrlType + ctrlType = ctrl.ctrlType ctrlType0 = ctrl0.ctrlType ncycles = ctrl0.numberOfCycles if initialControl == "charging" - if ctrlType0 == "discharging" && ctrlType == "charging" + if ctrlType0 == "discharging" && ctrlType == "charging" ncycles += 1 - end + end elseif initialControl == "discharging" - if ctrlType0 == "charging" && ctrlType == "discharging" + if ctrlType0 == "charging" && ctrlType == "discharging" ncycles += 1 - end + end else error("Initial control $initialControl is not recognized") - end + end - ctrl.numberOfCycles = ncycles + ctrl.numberOfCycles = ncycles return ctrl - end +end """ 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) From 689322464e030381bdad2e950a34ab87671f5164 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 13:18:37 +0200 Subject: [PATCH 13/27] Remove cccv. Only use CCCV --- src/input/schemas/lithium_ion/CyclingProtocol.json | 1 - 1 file changed, 1 deletion(-) diff --git a/src/input/schemas/lithium_ion/CyclingProtocol.json b/src/input/schemas/lithium_ion/CyclingProtocol.json index 33bbae64a..65dd2b7e1 100644 --- a/src/input/schemas/lithium_ion/CyclingProtocol.json +++ b/src/input/schemas/lithium_ion/CyclingProtocol.json @@ -7,7 +7,6 @@ "type": "string", "enum": [ "CC", - "cccv", "CCCV", "Function", "InputCurrentSeries", From 537e45bd561e1b2c1b2a13584795a19cde9dc3a7 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Wed, 17 Jun 2026 13:26:17 +0200 Subject: [PATCH 14/27] fix a few things from the merge --- src/simulation/simulation.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 336c7ab2d..5caf8bbd2 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -756,7 +756,7 @@ function solver_configuration( error("Timestep selector $(timestep_selector) not recognized. Only 'TimestepSelector' is currently implemented.") end - control_ramp_time = control_rampup_time(model[:Control].system.policy) + 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)) @@ -832,12 +832,12 @@ function solver_configuration( end end - if model[:Control].system.policy isa CyclingCVPolicy || model[:Control].system.policy isa CCPolicy || model[:Control].system.policy isa RestPolicy || model[:Control].system.policy isa SequencePolicy - if model[:Control].system.policy isa CyclingCVPolicy + 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 model[:Control].system.policy isa SequencePolicy + elseif multimodel[:Control].system.policy isa SequencePolicy cfg[:tolerances][:global_convergence_check_function] = (model, storage) -> check_constraints(model, storage) end @@ -884,10 +884,10 @@ function solver_configuration( report[:stopnow] = false end end - elseif model[:Control].system.policy isa RestPolicy + elseif multimodel[:Control].system.policy isa RestPolicy report[:stopnow] = s.state.Control.Controller.time >= m[:Control].system.policy.duration - elseif model[:Control].system.policy isa SequencePolicy + elseif multimodel[:Control].system.policy isa SequencePolicy report[:stopnow] = sequence_complete(m[:Control].system.policy, s.state.Control.Controller) end From 96753b6c40b4d17fc05a4be76c12ffcdeb36a1d5 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 12:40:39 +0200 Subject: [PATCH 15/27] add curve_analysis for rmse --- src/BattMo.jl | 1 + src/utils/curve_analysis.jl | 158 ++++++++++++++++++++++++++++++++++++ test/test_initial_state.jl | 5 +- 3 files changed, 160 insertions(+), 4 deletions(-) create mode 100644 src/utils/curve_analysis.jl 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/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/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 From 50568717b2737828e178a5df15512b2a1d3d5efd Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 12:40:49 +0200 Subject: [PATCH 16/27] Hopefully correct init state setup --- .../intercalation_battery.jl | 23 +++++++------------ 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/models/full_battery_models/intercalation_battery.jl b/src/models/full_battery_models/intercalation_battery.jl index 7a323ad37..4a4dec6db 100644 --- a/src/models/full_battery_models/intercalation_battery.jl +++ b/src/models/full_battery_models/intercalation_battery.jl @@ -1157,25 +1157,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 From c1772402f1a44806147aad7e427f0c64f52aafc5 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 15:02:36 +0200 Subject: [PATCH 17/27] work on example on seq control --- examples/example_sequence_control.jl | 87 +++++++++++++++++----------- 1 file changed, 52 insertions(+), 35 deletions(-) diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl index 9de7473ca..c0295484c 100644 --- a/examples/example_sequence_control.jl +++ b/examples/example_sequence_control.jl @@ -13,7 +13,7 @@ using Jutul: expand_to_ministeps, si_unit model_settings = load_model_settings(; from_default_set = "p2d") # Optionally disable ramp-up to see the effect of the initial ramp-up steps on the current and voltage profiles. -# delete!(model_settings, "RampUp") +delete!(model_settings, "RampUp") model = LithiumIonBattery(; model_settings = model_settings) cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") @@ -54,39 +54,15 @@ cycling_protocol = CyclingProtocol( sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) output = solve(sim; info_level = -1, output_substates = true) +states = output.states stored_states = output.jutul_output.states -raw_states, dt, = expand_to_ministeps(stored_states, output.jutul_output.reports[eachindex(stored_states)]) -time = cumsum(dt) -current = [only(state[:Control][:Current]) for state in raw_states] -voltage = [only(state[:Control][:ElectricPotential]) for state in raw_states] +raw_states, _, _ = expand_to_ministeps(stored_states, output.jutul_output.reports[eachindex(stored_states)]) doplot = true if doplot - using GLMakie - - sequence_step_index = [state[:Control][:Controller].step_index for state in raw_states] - sequence_time = [state[:Control][:Controller].time for state in raw_states] - step_boundary_times = Float64[] - - for i in 2:length(sequence_step_index) - if sequence_step_index[i] != sequence_step_index[i - 1] && sequence_step_index[i] <= length(cycling_protocol["Steps"]) - push!(step_boundary_times, sequence_time[i]) - end - end - hour = si_unit("hour") - time_hours = time ./ hour - dt_minutes = dt ./ si_unit("minute") - timestep_time_hours = Float64[] - timestep_values = Float64[] - for i in eachindex(dt) - push!(timestep_time_hours, (time[i] - dt[i]) / hour) - push!(timestep_time_hours, time_hours[i]) - push!(timestep_values, dt_minutes[i]) - push!(timestep_values, dt_minutes[i]) - end - boundary_hours = step_boundary_times ./ hour + using GLMakie fig = Figure(size = (900, 760)) ax_voltage = Axis(fig[1, 1], xlabel = "Time / h", ylabel = "Voltage / V") @@ -94,16 +70,57 @@ if doplot ax_timestep = Axis(fig[3, 1], xlabel = "Time / h", ylabel = "Timestep / min") linkxaxes!(ax_voltage, ax_current, ax_timestep) - scatterlines!(ax_voltage, time_hours, voltage, label = "Voltage", linewidth = 4.0, markersize = 10, markercolor = :black) - scatterlines!(ax_current, time_hours, current, label = "Current", linewidth = 4.0, markersize = 10, markercolor = :black) - lines!(ax_timestep, timestep_time_hours, timestep_values, label = "Timestep", linewidth = 3.0) + hour = si_unit("hour") + min = si_unit("minute") + + function t_mid(times) + return times[2:end] .- diff(times) ./ 2 + end - if !isempty(boundary_hours) - vlines!(ax_voltage, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) - vlines!(ax_current, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) - vlines!(ax_timestep, boundary_hours, color = :black, linestyle = :dash, linewidth = 2) + 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 + + # Plot standard output + scatterlines!(ax_voltage, output.time_series["Time"] ./ si_unit("hour"), output.time_series["Voltage"], label = "output.time_series", linewidth = 3.0, markersize = 8) + scatterlines!(ax_current, output.time_series["Time"] ./ si_unit("hour"), output.time_series["Current"], label = "output.time_series", linewidth = 3.0, markersize = 8) + # lines!(ax_timestep, t_mid(output.time_series["Time"]) ./ si_unit("hour"), diff(output.time_series["Time"]) ./ min, label = "output.time_series", linewidth = 3.0) + 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 = "output.time_series", linewidth = 3.0) + + # Plot expanded output (raw_states) + time = [only(state[:Control][:Controller].time) for state in raw_states] + current = [only(state[:Control][:Current]) for state in raw_states] + voltage = [only(state[:Control][:ElectricPotential]) for state in raw_states] + scatterlines!(ax_voltage, time ./ hour, voltage, label = "expanded", linewidth = 3.0, markersize = 8, marker = :circle) + scatterlines!(ax_current, time ./ hour, current, label = "expanded", linewidth = 3.0, markersize = 8, marker = :circle) + # lines!(ax_timestep, t_mid(time) ./ hour, diff(time) ./ min, label = "expanded", linewidth = 3.0, linestyle = :dash) + expanded_timestep_time, expanded_timestep_vals = t_pw_const(time) + lines!(ax_timestep, expanded_timestep_time ./ hour, expanded_timestep_vals ./ min, label = "expanded", linewidth = 3.0, linestyle = :dash) + + # Plot horizontal lines at step boundaries + sequence_step_index = [state[:Control][:Controller].step_index for state in raw_states] + sequence_time = [state[:Control][:Controller].time for state in raw_states] + step_boundary_times = Float64[] + for i in 2:length(sequence_step_index) + if sequence_step_index[i] != sequence_step_index[i - 1] && sequence_step_index[i] <= length(cycling_protocol["Steps"]) + push!(step_boundary_times, sequence_time[i]) + end end + vlines!(ax_voltage, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) + vlines!(ax_current, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) + vlines!(ax_timestep, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) screen = GLMakie.Screen() GLMakie.display(screen, fig) + end From 6fd08eeff603671327ccd40cd4e755f56bed7eab Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 16:03:55 +0200 Subject: [PATCH 18/27] Add CVCurrentCutoff option to CCCV to allow for absolute current cutoffs --- src/input/formatter.jl | 4 +++ src/input/meta_data/cycling_protocol.jl | 8 +++++ src/input/schemas/get_schema.jl | 1 + .../schemas/lithium_ion/CyclingProtocol.json | 8 +++++ src/input/schemas/parameters_meta_data.json | 8 +++++ src/models/current_and_voltage_boundary.jl | 25 +++++++++---- src/models/formatter.jl | 4 +++ .../intercalation_battery.jl | 12 +++++++ test/test_cycling_protocols.jl | 35 +++++++++++++++++++ 9 files changed, 99 insertions(+), 6 deletions(-) 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 4aca99097..bad8cad33 100644 --- a/src/input/meta_data/cycling_protocol.jl +++ b/src/input/meta_data/cycling_protocol.jl @@ -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, diff --git a/src/input/schemas/get_schema.jl b/src/input/schemas/get_schema.jl index 3e3866245..62acb6a52 100644 --- a/src/input/schemas/get_schema.jl +++ b/src/input/schemas/get_schema.jl @@ -472,6 +472,7 @@ 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"), diff --git a/src/input/schemas/lithium_ion/CyclingProtocol.json b/src/input/schemas/lithium_ion/CyclingProtocol.json index 65dd2b7e1..448589957 100644 --- a/src/input/schemas/lithium_ion/CyclingProtocol.json +++ b/src/input/schemas/lithium_ion/CyclingProtocol.json @@ -107,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, diff --git a/src/input/schemas/parameters_meta_data.json b/src/input/schemas/parameters_meta_data.json index 523bd884d..e143cea9e 100644 --- a/src/input/schemas/parameters_meta_data.json +++ b/src/input/schemas/parameters_meta_data.json @@ -581,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, diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 5dea07c0e..94948a532 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -221,6 +221,7 @@ mutable struct CyclingCVPolicy{R, I} <: AbstractPolicy upperCutoffVoltage::R dIdtLimit::R dEdtLimit::R + cvCurrentCutoff::Union{R, Missing} initialControl::OperationalMode numberOfCycles::I tolerances::Any @@ -242,6 +243,7 @@ function CyclingCVPolicy( use_ramp_up::Bool, rampup_time = zero(lowerCutoffVoltage), current_function = missing, + cv_current_cutoff = missing, ) if initialControl == "charging" @@ -252,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, @@ -266,6 +272,7 @@ function CyclingCVPolicy( upperCutoffVoltage, dIdtLimit, dEdtLimit, + cv_current_cutoff, initialControl, numberOfCycles, tolerances, @@ -1038,13 +1045,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 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 4a4dec6db..6c28f3b1d 100644 --- a/src/models/full_battery_models/intercalation_battery.jl +++ b/src/models/full_battery_models/intercalation_battery.jl @@ -239,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"], @@ -247,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" @@ -311,6 +317,11 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) 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"], @@ -322,6 +333,7 @@ function setup_control_model(input, model_neam, model_peam; T = Float64) 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 diff --git a/test/test_cycling_protocols.jl b/test/test_cycling_protocols.jl index 399b30de7..070bd830b 100644 --- a/test/test_cycling_protocols.jl +++ b/test/test_cycling_protocols.jl @@ -50,6 +50,41 @@ using Test @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 "Crate" begin @test begin From 19a923964beee7ca3711c52c71e00813ee683627 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 16:08:08 +0200 Subject: [PATCH 19/27] Change from OutputSubstrates/OutputSubstates to OutputMinisteps --- .../full_simulation_input/chen_2020.json | 4 +- .../full_simulation_input/chen_2020_p4d.json | 4 +- .../defaults/solver_settings/default.json | 4 +- .../defaults/solver_settings/direct.json | 4 +- .../defaults/solver_settings/iterative.json | 4 +- src/input/meta_data/solver_settings.jl | 6 +- src/input/schemas/get_schema.jl | 2 +- src/simulation/simulation.jl | 56 +++++++++++++++---- 8 files changed, 59 insertions(+), 25 deletions(-) diff --git a/src/input/defaults/full_simulation_input/chen_2020.json b/src/input/defaults/full_simulation_input/chen_2020.json index 24e7ecea9..3125e3b44 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 + "OutputMinisteps": 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..09281f918 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 + "OutputMinisteps": 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..1cdf07b09 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 + "OutputMinisteps": 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..1cdf07b09 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 + "OutputMinisteps": 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..41d60d91d 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 + "OutputMinisteps": false } -} \ No newline at end of file +} diff --git a/src/input/meta_data/solver_settings.jl b/src/input/meta_data/solver_settings.jl index 3453e53d2..857491d8e 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( + "OutputMinisteps" => Dict( "type" => Bool, - "variable_name" => "output_substates", + "variable_name" => "output_ministeps", "category" => "SolverSettings", - "description" => "Store substates (between report steps) as field on each state.", + "description" => "Return output at every accepted solver ministep 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 62acb6a52..8ae37eb10 100644 --- a/src/input/schemas/get_schema.jl +++ b/src/input/schemas/get_schema.jl @@ -756,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"), + "OutputMinisteps" => create_property(parameter_meta, "OutputMinisteps"), ), ), ), diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 5caf8bbd2..ccc6d30eb 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -311,6 +311,7 @@ end info_level = 0, end_report = info_level > -1, include_initial_state = false, + output_ministeps = false, kwargs...) Solves a battery `Simulation` problem by executing the simulation workflow defined in `solve_simulation`. @@ -325,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_ministeps::Bool` (optional): If `true`, the returned BattMo output is sampled at + every accepted solver ministep 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 @@ -353,6 +357,7 @@ function solve( solver_settings = get_default_solver_settings(problem.model), logger = nothing, include_initial_state = false, + output_ministeps::Bool = false, kwargs..., ) @@ -385,6 +390,7 @@ function solve( solver_settings, logger, include_initial_state, + output_ministeps, validate, accept_invalid, kwargs..., @@ -432,6 +438,7 @@ function solve_simulation( solver_settings, logger = nothing, include_initial_state = false, + output_ministeps::Bool = false, validate::Bool = true, accept_invalid::Bool = false, kwargs..., @@ -471,6 +478,12 @@ function solve_simulation( kwargs..., ) + # Jutul stores accepted ministep states using its internal output_substates option. + output_ministeps = output_ministeps || cfg[:output_substates] + if output_ministeps + cfg[:output_substates] = true + end + # Setup hook if given hook = get(kwargs, :hook, nothing) @@ -488,6 +501,11 @@ function solve_simulation( jutul_states, jutul_reports = simulate(state0, simulator, timesteps; forces = forces, config = cfg, kwargs...) + output_states = jutul_states + if output_ministeps + output_states, _, _ = Jutul.expand_to_ministeps(jutul_states, jutul_reports[eachindex(jutul_states)]) + end + # perfom decoupled thermal simulation if needed @@ -497,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_ministeps + 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_ministeps + jutul_reports[i][:ThermalModel] = decoupled_thermal_states[i][:report] + end end @@ -511,7 +536,7 @@ function solve_simulation( jutul_output = ( - states = jutul_states, + states = output_states, reports = jutul_reports, solver_configuration = cfg, multimodel = model.multimodel, @@ -660,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), + "OutputMinisteps" => get(kwargs, :output_ministeps, nothing), ), ) @@ -739,6 +764,7 @@ function solver_configuration( linear_solver_dict = solver_settings["LinearSolver"] output = solver_settings["Output"] verbose = solver_settings["Verbose"] + output_ministeps = get(output, "OutputMinisteps", false) relaxation = non_linear_solver["Relaxation"] if relaxation == "NoRelaxation" @@ -800,7 +826,7 @@ 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_ministeps, ) if !isnothing(rampup_reference_timestep) && !isnothing(rampup_steps) && rampup_steps > 0 && control_ramp_time > 0 @@ -1074,12 +1100,20 @@ function setup_timesteps( 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"] From 12875c98c4d66785945c58bee47e071ecb36a33d Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 16:08:44 +0200 Subject: [PATCH 20/27] Demonstate further the effect of rampup and ministeps output --- examples/example_sequence_control.jl | 227 ++++++++++++++++++++------- 1 file changed, 174 insertions(+), 53 deletions(-) diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl index c0295484c..6d642bb99 100644 --- a/examples/example_sequence_control.jl +++ b/examples/example_sequence_control.jl @@ -1,5 +1,5 @@ using BattMo -using Jutul: expand_to_ministeps, si_unit +using Jutul: si_unit ### # Example: Sequence control with CC, Rest, and CCCV steps @@ -10,15 +10,9 @@ using Jutul: expand_to_ministeps, si_unit # applied to the current-controlled parts of the CC and CCCV steps. ### -model_settings = load_model_settings(; from_default_set = "p2d") - -# Optionally disable ramp-up to see the effect of the initial ramp-up steps on the current and voltage profiles. -delete!(model_settings, "RampUp") - -model = LithiumIonBattery(; model_settings = model_settings) cell_parameters = load_cell_parameters(; from_default_set = "chen_2020") -simulation_settings = load_simulation_settings(; from_default_set = "p2d") -simulation_settings["TimeStepDuration"] = 360.0 +base_simulation_settings = load_simulation_settings(; from_default_set = "p2d") +base_simulation_settings["TimeStepDuration"] = 360.0 cc = Dict( "Protocol" => "CC", @@ -52,12 +46,30 @@ cycling_protocol = CyclingProtocol( ) ) -sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) -output = solve(sim; info_level = -1, output_substates = true) -states = output.states +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 -stored_states = output.jutul_output.states -raw_states, _, _ = expand_to_ministeps(stored_states, output.jutul_output.reports[eachindex(stored_states)]) + model = LithiumIonBattery(; model_settings = model_settings) + sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) + output = solve(sim; info_level = -1, output_ministeps = false) + ministep_output = solve(sim; info_level = -1, output_ministeps = true) + + return ( + label = label, + output = output, + ministep_output = ministep_output, + ) +end + +results = [run_sequence_case(use_rampup) for use_rampup in (true, false)] doplot = true if doplot @@ -65,18 +77,9 @@ if doplot using GLMakie fig = Figure(size = (900, 760)) - ax_voltage = Axis(fig[1, 1], xlabel = "Time / h", ylabel = "Voltage / V") - ax_current = Axis(fig[2, 1], xlabel = "Time / h", ylabel = "Current / A") - ax_timestep = Axis(fig[3, 1], xlabel = "Time / h", ylabel = "Timestep / min") - linkxaxes!(ax_voltage, ax_current, ax_timestep) - hour = si_unit("hour") min = si_unit("minute") - function t_mid(times) - return times[2:end] .- diff(times) ./ 2 - end - function t_pw_const(times) timestep_time = Float64[] timestep_vals = Float64[] @@ -90,37 +93,155 @@ if doplot return timestep_time, timestep_vals end - # Plot standard output - scatterlines!(ax_voltage, output.time_series["Time"] ./ si_unit("hour"), output.time_series["Voltage"], label = "output.time_series", linewidth = 3.0, markersize = 8) - scatterlines!(ax_current, output.time_series["Time"] ./ si_unit("hour"), output.time_series["Current"], label = "output.time_series", linewidth = 3.0, markersize = 8) - # lines!(ax_timestep, t_mid(output.time_series["Time"]) ./ si_unit("hour"), diff(output.time_series["Time"]) ./ min, label = "output.time_series", linewidth = 3.0) - 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 = "output.time_series", linewidth = 3.0) - - # Plot expanded output (raw_states) - time = [only(state[:Control][:Controller].time) for state in raw_states] - current = [only(state[:Control][:Current]) for state in raw_states] - voltage = [only(state[:Control][:ElectricPotential]) for state in raw_states] - scatterlines!(ax_voltage, time ./ hour, voltage, label = "expanded", linewidth = 3.0, markersize = 8, marker = :circle) - scatterlines!(ax_current, time ./ hour, current, label = "expanded", linewidth = 3.0, markersize = 8, marker = :circle) - # lines!(ax_timestep, t_mid(time) ./ hour, diff(time) ./ min, label = "expanded", linewidth = 3.0, linestyle = :dash) - expanded_timestep_time, expanded_timestep_vals = t_pw_const(time) - lines!(ax_timestep, expanded_timestep_time ./ hour, expanded_timestep_vals ./ min, label = "expanded", linewidth = 3.0, linestyle = :dash) - - # Plot horizontal lines at step boundaries - sequence_step_index = [state[:Control][:Controller].step_index for state in raw_states] - sequence_time = [state[:Control][:Controller].time for state in raw_states] - step_boundary_times = Float64[] - for i in 2:length(sequence_step_index) - if sequence_step_index[i] != sequence_step_index[i - 1] && sequence_step_index[i] <= length(cycling_protocol["Steps"]) - push!(step_boundary_times, sequence_time[i]) - end + report_colors = [:dodgerblue3, :darkorange2] + ministep_colors = [:seagreen4, :purple3] + + ax_voltage = Axis( + fig[1, 1], + title = "Voltage: report output and ministep output", + xlabel = "Time / h", + ylabel = "Voltage / V", + ) + ax_current = Axis( + fig[2, 1], + title = "Current: report output and ministep 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 + ministep_output = result.ministep_output + report_color = report_colors[case_index] + ministep_color = ministep_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_ministeps=true shows accepted solver ministeps, including adaptive + # timestep cuts from convergence and control transitions, not only ramp-up steps. + time = ministep_output.time_series["Time"] + scatterlines!( + ax_voltage, + time ./ hour, + ministep_output.time_series["Voltage"], + label = "$label: ministeps", + color = ministep_color, + linewidth = 3.0, + markersize = 8, + marker = :circle, + ) + scatterlines!( + ax_current, + time ./ hour, + ministep_output.time_series["Current"], + label = "$label: ministeps", + color = ministep_color, + linewidth = 3.0, + markersize = 8, + marker = :circle, + ) + ministep_timestep_time, ministep_timestep_vals = t_pw_const(time) + lines!( + ax_timestep, + ministep_timestep_time ./ hour, + ministep_timestep_vals ./ min, + label = "$label: ministeps", + color = ministep_color, + linewidth = 3.0, + ) end - vlines!(ax_voltage, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) - vlines!(ax_current, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) - vlines!(ax_timestep, step_boundary_times ./ hour, color = :black, linestyle = :dash, linewidth = 2) - screen = GLMakie.Screen() - GLMakie.display(screen, fig) + 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].ministep_output, "W rampup: ministeps") + plot_voltage_current!(fig_cases[2, 1], results[2].output, "WO rampup: output") + plot_voltage_current!(fig_cases[2, 2], results[2].ministep_output, "WO rampup: ministeps") + + screen1 = GLMakie.Screen() + screen2 = GLMakie.Screen() + GLMakie.display(screen1, fig) + GLMakie.display(screen2, fig_cases) end From 7cb8714fc09b9f225f00b8e67a2566df09183757 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 16:09:10 +0200 Subject: [PATCH 21/27] add default rampup option --- src/models/current_and_voltage_boundary.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index 94948a532..fc6691a1b 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -240,7 +240,7 @@ 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, From 8d2708d828de3a4eb477f129b325902b7f0e0494 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Fri, 19 Jun 2026 16:10:39 +0200 Subject: [PATCH 22/27] Test for mini-rampup implemented for CC --- test/test_cycling_protocols.jl | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/test_cycling_protocols.jl b/test/test_cycling_protocols.jl index 070bd830b..6ed305f28 100644 --- a/test/test_cycling_protocols.jl +++ b/test/test_cycling_protocols.jl @@ -85,6 +85,33 @@ end @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 From 65eef4cacfc8bf6fc9ff902c928eb9fa7cfb48f7 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Sat, 20 Jun 2026 10:30:13 +0200 Subject: [PATCH 23/27] Change ministeps -> substates to conform with Jutul --- examples/example_sequence_control.jl | 46 +++++++++---------- .../full_simulation_input/chen_2020.json | 2 +- .../full_simulation_input/chen_2020_p4d.json | 2 +- .../defaults/solver_settings/default.json | 2 +- .../defaults/solver_settings/direct.json | 2 +- .../defaults/solver_settings/iterative.json | 2 +- src/input/meta_data/solver_settings.jl | 6 +-- src/input/schemas/get_schema.jl | 2 +- src/simulation/simulation.jl | 30 ++++++------ 9 files changed, 47 insertions(+), 47 deletions(-) diff --git a/examples/example_sequence_control.jl b/examples/example_sequence_control.jl index 6d642bb99..c323f554c 100644 --- a/examples/example_sequence_control.jl +++ b/examples/example_sequence_control.jl @@ -59,13 +59,13 @@ function run_sequence_case(use_rampup::Bool) model = LithiumIonBattery(; model_settings = model_settings) sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) - output = solve(sim; info_level = -1, output_ministeps = false) - ministep_output = solve(sim; info_level = -1, output_ministeps = true) + output = solve(sim; info_level = -1, output_substates = false) + substate_output = solve(sim; info_level = -1, output_substates = true) return ( label = label, output = output, - ministep_output = ministep_output, + substate_output = substate_output, ) end @@ -94,17 +94,17 @@ if doplot end report_colors = [:dodgerblue3, :darkorange2] - ministep_colors = [:seagreen4, :purple3] + substate_colors = [:seagreen4, :purple3] ax_voltage = Axis( fig[1, 1], - title = "Voltage: report output and ministep output", + title = "Voltage: report output and substate output", xlabel = "Time / h", ylabel = "Voltage / V", ) ax_current = Axis( fig[2, 1], - title = "Current: report output and ministep output", + title = "Current: report output and substate output", xlabel = "Time / h", ylabel = "Current / A", ) @@ -119,9 +119,9 @@ if doplot for (case_index, result) in enumerate(results) label = result.label output = result.output - ministep_output = result.ministep_output + substate_output = result.substate_output report_color = report_colors[case_index] - ministep_color = ministep_colors[case_index] + substate_color = substate_colors[case_index] scatter!( ax_voltage, @@ -149,15 +149,15 @@ if doplot linewidth = 3.0, ) - # output_ministeps=true shows accepted solver ministeps, including adaptive + # output_substates=true shows accepted solver substates, including adaptive # timestep cuts from convergence and control transitions, not only ramp-up steps. - time = ministep_output.time_series["Time"] + time = substate_output.time_series["Time"] scatterlines!( ax_voltage, time ./ hour, - ministep_output.time_series["Voltage"], - label = "$label: ministeps", - color = ministep_color, + substate_output.time_series["Voltage"], + label = "$label: substates", + color = substate_color, linewidth = 3.0, markersize = 8, marker = :circle, @@ -165,20 +165,20 @@ if doplot scatterlines!( ax_current, time ./ hour, - ministep_output.time_series["Current"], - label = "$label: ministeps", - color = ministep_color, + substate_output.time_series["Current"], + label = "$label: substates", + color = substate_color, linewidth = 3.0, markersize = 8, marker = :circle, ) - ministep_timestep_time, ministep_timestep_vals = t_pw_const(time) + substate_timestep_time, substate_timestep_vals = t_pw_const(time) lines!( ax_timestep, - ministep_timestep_time ./ hour, - ministep_timestep_vals ./ min, - label = "$label: ministeps", - color = ministep_color, + substate_timestep_time ./ hour, + substate_timestep_vals ./ min, + label = "$label: substates", + color = substate_color, linewidth = 3.0, ) end @@ -235,9 +235,9 @@ if doplot end plot_voltage_current!(fig_cases[1, 1], results[1].output, "W rampup: output") - plot_voltage_current!(fig_cases[1, 2], results[1].ministep_output, "W rampup: ministeps") + 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].ministep_output, "WO rampup: ministeps") + plot_voltage_current!(fig_cases[2, 2], results[2].substate_output, "WO rampup: substates") screen1 = GLMakie.Screen() screen2 = GLMakie.Screen() diff --git a/src/input/defaults/full_simulation_input/chen_2020.json b/src/input/defaults/full_simulation_input/chen_2020.json index 3125e3b44..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, - "OutputMinisteps": false + "OutputSubstates": false } } } 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 09281f918..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, - "OutputMinisteps": false + "OutputSubstates": false } } } diff --git a/src/input/defaults/solver_settings/default.json b/src/input/defaults/solver_settings/default.json index 1cdf07b09..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, - "OutputMinisteps": false + "OutputSubstates": false } } diff --git a/src/input/defaults/solver_settings/direct.json b/src/input/defaults/solver_settings/direct.json index 1cdf07b09..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, - "OutputMinisteps": false + "OutputSubstates": false } } diff --git a/src/input/defaults/solver_settings/iterative.json b/src/input/defaults/solver_settings/iterative.json index 41d60d91d..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, - "OutputMinisteps": false + "OutputSubstates": false } } diff --git a/src/input/meta_data/solver_settings.jl b/src/input/meta_data/solver_settings.jl index 857491d8e..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.", ), - "OutputMinisteps" => Dict( + "OutputSubstates" => Dict( "type" => Bool, - "variable_name" => "output_ministeps", + "variable_name" => "output_substates", "category" => "SolverSettings", - "description" => "Return output at every accepted solver ministep instead of only at report timesteps.", + "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 8ae37eb10..f0c9b51ef 100644 --- a/src/input/schemas/get_schema.jl +++ b/src/input/schemas/get_schema.jl @@ -756,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"), - "OutputMinisteps" => create_property(parameter_meta, "OutputMinisteps"), + "OutputSubstates" => create_property(parameter_meta, "OutputSubstates"), ), ), ), diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index ccc6d30eb..02de4792c 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -311,7 +311,7 @@ end info_level = 0, end_report = info_level > -1, include_initial_state = false, - output_ministeps = false, + output_substates = false, kwargs...) Solves a battery `Simulation` problem by executing the simulation workflow defined in `solve_simulation`. @@ -326,8 +326,8 @@ 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_ministeps::Bool` (optional): If `true`, the returned BattMo output is sampled at - every accepted solver ministep instead of only at report timesteps. The output format is +- `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`. @@ -357,7 +357,7 @@ function solve( solver_settings = get_default_solver_settings(problem.model), logger = nothing, include_initial_state = false, - output_ministeps::Bool = false, + output_substates::Bool = false, kwargs..., ) @@ -390,7 +390,7 @@ function solve( solver_settings, logger, include_initial_state, - output_ministeps, + output_substates, validate, accept_invalid, kwargs..., @@ -438,7 +438,7 @@ function solve_simulation( solver_settings, logger = nothing, include_initial_state = false, - output_ministeps::Bool = false, + output_substates::Bool = false, validate::Bool = true, accept_invalid::Bool = false, kwargs..., @@ -478,9 +478,9 @@ function solve_simulation( kwargs..., ) - # Jutul stores accepted ministep states using its internal output_substates option. - output_ministeps = output_ministeps || cfg[:output_substates] - if output_ministeps + # 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 @@ -502,7 +502,7 @@ function solve_simulation( jutul_states, jutul_reports = simulate(state0, simulator, timesteps; forces = forces, config = cfg, kwargs...) output_states = jutul_states - if output_ministeps + if output_substates output_states, _, _ = Jutul.expand_to_ministeps(jutul_states, jutul_reports[eachindex(jutul_states)]) end @@ -515,7 +515,7 @@ function solve_simulation( thermal_parameters = thermal_sim.parameters thermal_model = thermal_sim.model - thermal_timesteps = if output_ministeps + thermal_timesteps = if output_substates diff(vcat(0.0, [state[:Control][:Controller].time for state in output_states])) else timesteps @@ -526,7 +526,7 @@ function solve_simulation( # 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_ministeps + if !output_substates jutul_reports[i][:ThermalModel] = decoupled_thermal_states[i][:report] end end @@ -685,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), - "OutputMinisteps" => get(kwargs, :output_ministeps, nothing), + "OutputSubstates" => get(kwargs, :output_substates, nothing), ), ) @@ -764,7 +764,7 @@ function solver_configuration( linear_solver_dict = solver_settings["LinearSolver"] output = solver_settings["Output"] verbose = solver_settings["Verbose"] - output_ministeps = get(output, "OutputMinisteps", false) + output_substates_setting = get(output, "OutputSubstates", false) relaxation = non_linear_solver["Relaxation"] if relaxation == "NoRelaxation" @@ -826,7 +826,7 @@ function solver_configuration( output_path = output["OutputPath"] == "" ? nothing : output["OutputPath"], in_memory_reports = output["InMemoryReports"], report_level = output["ReportLevel"], - output_substates = output_ministeps, + output_substates = output_substates_setting, ) if !isnothing(rampup_reference_timestep) && !isnothing(rampup_steps) && rampup_steps > 0 && control_ramp_time > 0 From 134fa35915d05ce4e753ae38be4cbd6314ae779f Mon Sep 17 00:00:00 2001 From: August Johansson Date: Sat, 20 Jun 2026 23:45:08 +0200 Subject: [PATCH 24/27] Fix minor bug regarding using non-post-processed states --- src/simulation/simulation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 02de4792c..fdaaca5c4 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -499,11 +499,11 @@ 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 = jutul_states + output_states = tmp_states if output_substates - output_states, _, _ = Jutul.expand_to_ministeps(jutul_states, jutul_reports[eachindex(jutul_states)]) + output_states, _, _ = Jutul.expand_to_ministeps(tmp_states, jutul_reports[eachindex(tmp_states)]) end @@ -557,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, From 70077f0f5493c5220bb7d0fad843ed531d2053b0 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Sat, 20 Jun 2026 23:45:42 +0200 Subject: [PATCH 25/27] Comment --- src/models/current_and_voltage_boundary.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/models/current_and_voltage_boundary.jl b/src/models/current_and_voltage_boundary.jl index fc6691a1b..352283c9a 100644 --- a/src/models/current_and_voltage_boundary.jl +++ b/src/models/current_and_voltage_boundary.jl @@ -483,6 +483,8 @@ mutable struct SequenceController{R} <: Controller 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 From 4fa77495becfb0663ec9176006ecef79d8b449d4 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Sat, 20 Jun 2026 23:45:53 +0200 Subject: [PATCH 26/27] Update test --- test/test_sequence_protocols.jl | 129 ++++++++++++++++++++++++++++---- 1 file changed, 113 insertions(+), 16 deletions(-) diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl index f462ac1ae..3c3ead2bc 100644 --- a/test/test_sequence_protocols.jl +++ b/test/test_sequence_protocols.jl @@ -1,6 +1,23 @@ 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") @@ -22,6 +39,60 @@ using Test @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 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") @@ -58,49 +129,75 @@ end @test any(abs.(current) .< 1.0e-8) end -@testset "Sequence protocol accepts combined CC Rest and CCCV steps" begin +@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"] = 50 + simulation_settings["TimeStepDuration"] = 20.0 model = LithiumIonBattery(; model_settings) cycling_protocol = CyclingProtocol( Dict( "Protocol" => "Sequence", - "InitialStateOfCharge" => 0.5, + "InitialStateOfCharge" => 1.0, "Steps" => [ Dict( "Protocol" => "CC", "InitialControl" => "discharging", - "DRate" => 0.1, + "DRate" => 1.0, "TotalNumberOfCycles" => 0, - "LowerVoltageLimit" => 3.9, - "UpperVoltageLimit" => 4.1, + "LowerVoltageLimit" => 3.0, + "UpperVoltageLimit" => 4.2, ), Dict( "Protocol" => "Rest", - "Duration" => 20.0, + "Duration" => 120.0, ), Dict( - "Protocol" => "CCCV", + "Protocol" => "CC", "InitialControl" => "charging", "CRate" => 1.0, - "DRate" => 1.0, - "TotalNumberOfCycles" => 1, + "TotalNumberOfCycles" => 0, "LowerVoltageLimit" => 3.0, - "UpperVoltageLimit" => 4.0, - "CurrentChangeLimit" => 1.0e-4, - "VoltageChangeLimit" => 1.0e-4, + "UpperVoltageLimit" => 3.8, ), ], ) ) - sim = Simulation(model, cell_parameters, cycling_protocol; simulation_settings) + 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) - @test !isempty(sim.time_steps) - @test length(cycling_protocol["Steps"]) == 3 + 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 From 8b089911441528b338cbd032244cc4fda8e02e49 Mon Sep 17 00:00:00 2001 From: August Johansson Date: Sat, 20 Jun 2026 23:54:38 +0200 Subject: [PATCH 27/27] Add test of checkbeforesolve necessity --- test/test_sequence_protocols.jl | 71 +++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/test/test_sequence_protocols.jl b/test/test_sequence_protocols.jl index 3c3ead2bc..72427f685 100644 --- a/test/test_sequence_protocols.jl +++ b/test/test_sequence_protocols.jl @@ -93,6 +93,77 @@ end 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")