Degraded Battery Cell Design Optimization¶
This example demonstrates how to optimize battery cell design parameters using ionworkspipeline. We optimize active material volume fractions to maximize discharge capacity while preventing lithium plating by constraining the anode potential at cycle 50 in the batteries’ lifetime.
import ionworkspipeline as iwp
from ionworkspipeline.data_fits.models.metrics import Maximum, Minimum
from ionworkspipeline.data_fits.models.metrics.actions import GreaterThan, Maximize
import matplotlib.pyplot as plt
import numpy as np
import pybamm
from pybamm import Parameter
Configure DFN Model and Optimization Parameters¶
Import Libraries and Dependencies
# Setup model and baseline parameters
model = pybamm.lithium_ion.SPMe(
options={
"SEI": "solvent-diffusion limited",
"SEI porosity change": "false",
"particle mechanics": "swelling only",
"loss of active material": "stress and reaction-driven",
}
)
# Initialize parameters and update LAM rate
params_baseline = pybamm.ParameterValues("Ai2020")
params_baseline.update(
{"Negative electrode LAM constant proportional term [s-1]": 3e-1 / 3600}
)
params_baseline.update(
{"Positive electrode LAM constant proportional term [s-1]": 3e-1 / 3600}
)
# Link porosity to active material fraction
for electrode in ["Positive", "Negative"]:
params_baseline[f"{electrode} electrode porosity"] = 1.0 - Parameter(
f"{electrode} electrode active material volume fraction"
)
# Define optimization parameters
parameters = {
"Positive electrode active material volume fraction": iwp.Parameter(
"Positive electrode active material volume fraction",
initial_value=0.65,
bounds=(0.5, 0.9),
),
"Negative electrode active material volume fraction": iwp.Parameter(
"Negative electrode active material volume fraction",
initial_value=0.65,
bounds=(0.5, 0.85),
),
}
print(f"Model configured with {len(parameters)} optimization parameters")
Model configured with 2 optimization parameters
Define Experiment Protocol and Optimization Metrics¶
# Define experiment protocol
experiment = pybamm.Experiment(
[
(
"Discharge at 2C until 3.00 V",
"Charge at 2C until 4.20 V",
)
]
* 50
)
# Define optimization actions
actions = {
"Maximize discharge capacity [A.h]": Maximize(Maximum("Discharge capacity [A.h]")),
}
# Constraints define bounds that must be satisfied (using GreaterThan/LessThan actions)
constraints = {
"Anode potential minimum constraint [V]": GreaterThan(
Minimum(
"Negative electrode surface potential difference at separator interface [V]"
),
value=0.0,
penalty=150.0,
),
}
Execute Battery Design Optimization¶
Run Validation Simulations¶
# Setup and run optimization
initial_soc = 1.0
params_for_pipeline = {k: v for k, v in params_baseline.items() if k not in parameters}
solver = pybamm.IDAKLUSolver(
options={"num_steps_no_progress": 1000, "t_no_progress": 0.99}
)
objective_options = {
"model": model,
"simulation_kwargs": {
"experiment": experiment,
"solver": solver,
},
}
objective = iwp.objectives.DesignObjective(
options=objective_options,
actions=actions,
constraints=constraints,
save_at_cycles=[50],
validate_against_experiment_steps=False,
parameters={"Initial SOC [%]": 100 * initial_soc},
)
optim = iwp.optimizers.Pints(
max_iterations=10, # Increase for better results
max_unchanged_iterations=30,
max_unchanged_iterations_threshold=1e-4,
log_to_screen=True,
)
design_datafit = iwp.DataFit(
objective, parameters=parameters, options={"seed": 0}, optimizer=optim
)
print("Running optimization...")
results = design_datafit.run(params_for_pipeline)
print(f"\nOptimization completed! Final cost: {results.costs:.6f}")
print("Optimal parameters:")
for name in parameters:
print(f" {name}: {results.parameter_values[name]:.6f}")
Running optimization...
Minimising error measure
Using CMAES
Population size: 6
--------------------------------------------------------------------------------
Iter. Eval. Best Current Time m:s
--------------------------------------------------------------------------------
1 6 1.485628e+02 1.485628e+02 0:09
2 12 1.483528e+02 1.483528e+02 0:18
3 18 -2.233478e+00 -2.233478e+00 0:26
4 24 -2.233478e+00 -2.233478e+00 0:36
5 30 -2.233478e+00 -2.233478e+00 0:48
10 60 -2.463221e+00 -2.463221e+00 1:38
--------------------------------------------------------------------------------
Maximum number of iterations (10) reached.
Optimization complete
Optimization completed! Final cost: -2.463221
Optimal parameters:
Positive electrode active material volume fraction: 0.898269
Negative electrode active material volume fraction: 0.652055
# Helper function to extract metrics from solution
def extract_metrics(solution):
sol_discharge = solution.sub_solutions[0]
discharge_capacity = abs(sol_discharge["Discharge capacity [A.h]"]()[-1])
# Combine charge steps
charge_solutions = solution.sub_solutions[2:]
sol_charge = charge_solutions[0]
for charge_sol in charge_solutions[1:]:
sol_charge = sol_charge + charge_sol
anode_potential_series = sol_charge[
"Negative electrode surface potential difference at separator interface [V]"
]()
min_anode_potential = min(anode_potential_series)
mean_anode_potential = anode_potential_series.mean()
return {
"discharge_capacity": discharge_capacity,
"min_anode_potential": min_anode_potential,
"mean_anode_potential": mean_anode_potential,
"solution": solution,
}
# Setup baseline parameters for comparison
params_optimal = params_baseline.copy()
params_optimal.update(results.parameter_values.copy())
params_baseline_original = params_baseline.copy()
# Run simulations for comparison
sim_baseline = iwp.Simulation(
model=model, parameter_values=params_baseline_original, experiment=experiment
)
sim_optimal = iwp.Simulation(
model=model, parameter_values=params_optimal, experiment=experiment
)
sol_baseline = sim_baseline.solve(initial_soc=initial_soc)
sol_optimal = sim_optimal.solve(initial_soc=initial_soc)
baseline_metrics = extract_metrics(sol_baseline)
optimal_metrics = extract_metrics(sol_optimal)
print("Validation simulations completed")
Validation simulations completed
Analyze Performance Metrics and Parameter Changes¶
# Performance comparison and analysis
capacity_improvement = (
(optimal_metrics["discharge_capacity"] - baseline_metrics["discharge_capacity"])
/ baseline_metrics["discharge_capacity"]
) * 100
min_anode_improvement = (
optimal_metrics["min_anode_potential"] - baseline_metrics["min_anode_potential"]
)
print("Performance Comparison:")
print("=" * 50)
print(f"{'Metric':<35} {'Baseline':<12} {'Optimal':<12} {'Improvement'}")
print("-" * 70)
print(
f"{'Discharge Capacity [A.h]':<35} {baseline_metrics['discharge_capacity']:<12.4f} {optimal_metrics['discharge_capacity']:<12.4f} {capacity_improvement:+.2f}%"
)
print(
f"{'Min Anode Potential [V]':<35} {baseline_metrics['min_anode_potential']:<12.4f} {optimal_metrics['min_anode_potential']:<12.4f} {min_anode_improvement:+.4f}"
)
print(
f"{'Mean Anode Potential [V]':<35} {baseline_metrics['mean_anode_potential']:<12.4f} {optimal_metrics['mean_anode_potential']:<12.4f}"
)
print("=" * 50)
print("\nParameter Changes:")
for name, param in parameters.items():
baseline_val = param.initial_value
optimal_val = results.parameter_values[name]
change = optimal_val - baseline_val
percent_change = (change / baseline_val) * 100
print(f"{name}: {baseline_val:.6f} → {optimal_val:.6f} ({percent_change:+.2f}%)")
Performance Comparison:
==================================================
Metric Baseline Optimal Improvement
----------------------------------------------------------------------
Discharge Capacity [A.h] 2.3748 2.5682 +8.14%
Min Anode Potential [V] -0.0016 0.0032 +0.0048
Mean Anode Potential [V] 0.1774 0.2192
==================================================
Parameter Changes:
Positive electrode active material volume fraction: 0.650000 → 0.898269 (+38.20%)
Negative electrode active material volume fraction: 0.650000 → 0.652055 (+0.32%)
Create Comparative Visualizations¶
# Create comprehensive visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# Voltage comparison
t_baseline = baseline_metrics["solution"].t / 3600
t_optimal = optimal_metrics["solution"].t / 3600
axes[0, 0].plot(
t_baseline,
baseline_metrics["solution"]["Voltage [V]"](baseline_metrics["solution"].t),
label="Baseline",
linewidth=2,
alpha=0.8,
)
axes[0, 0].plot(
t_optimal,
optimal_metrics["solution"]["Voltage [V]"](optimal_metrics["solution"].t),
label="Optimal",
linewidth=2,
alpha=0.8,
linestyle="--",
)
axes[0, 0].set_xlabel("Time [h]")
axes[0, 0].set_ylabel("Voltage [V]")
axes[0, 0].set_title("Cell Voltage: Baseline vs Optimal")
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# Discharge capacity
metrics_names = [
"Discharge\nCapacity @ Cycle 1",
f"Discharge\nCapacity @ Cycle {len(sol_baseline.cycles)}",
]
baseline_values = [
baseline_metrics["solution"].cycles[0]["Discharge capacity [A.h]"]().max(),
baseline_metrics["solution"].cycles[-1]["Discharge capacity [A.h]"]().max(),
]
optimal_values = [
optimal_metrics["solution"].cycles[0]["Discharge capacity [A.h]"]().max(),
optimal_metrics["solution"].cycles[-1]["Discharge capacity [A.h]"]().max(),
]
x = np.arange(len(metrics_names))
width = 0.35
axes[0, 1].bar(
x - width / 2,
baseline_values,
width,
label="Baseline",
alpha=0.8,
color="lightblue",
)
axes[0, 1].bar(
x + width / 2, optimal_values, width, label="Optimal", alpha=0.8, color="orange"
)
for i, (baseline, optimal) in enumerate(
zip(baseline_values, optimal_values, strict=False)
):
y_offset = baseline * 0.02 if baseline > 0 else 0.005
axes[0, 1].text(
i - width / 2,
baseline + y_offset,
f"{baseline:.3f}",
ha="center",
va="bottom",
fontsize=9,
)
axes[0, 1].text(
i + width / 2,
optimal + y_offset,
f"{optimal:.3f}",
ha="center",
va="bottom",
fontsize=9,
)
axes[0, 1].set_xlabel("Capacity [A.h]")
axes[0, 1].set_ylabel("Value")
axes[0, 1].set_title("Degraded Design Optimization")
axes[0, 1].set_xticks(x)
axes[0, 1].set_xticklabels(metrics_names)
axes[0, 1].set_ylim(top=3.4)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# Anode potential during charge
anode_baseline = sol_baseline[
"Negative electrode surface potential difference at separator interface [V]"
](sol_baseline.t)
anode_optimal = sol_optimal[
"Negative electrode surface potential difference at separator interface [V]"
](sol_optimal.t)
axes[1, 0].plot(t_baseline, anode_baseline, label="Baseline", linewidth=2, alpha=0.8)
axes[1, 0].plot(
t_optimal,
anode_optimal,
label="Optimal",
linewidth=2,
alpha=0.8,
linestyle="--",
)
axes[1, 0].axhline(y=0.0, color="r", linestyle=":", alpha=0.7, label="Safety limit")
axes[1, 0].set_xlabel("Charge Time [h]")
axes[1, 0].set_ylabel("Anode - Electrolyte Potential [V]")
axes[1, 0].set_title("Anode Potential During Charge")
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# Performance metrics bar chart
metrics_names = [
"Discharge\nCapacity [A.h]",
"Min Anode\nPotential [V]",
"Mean Anode\nPotential [V]",
]
baseline_values = [
baseline_metrics["discharge_capacity"],
baseline_metrics["min_anode_potential"],
baseline_metrics["mean_anode_potential"],
]
optimal_values = [
optimal_metrics["discharge_capacity"],
optimal_metrics["min_anode_potential"],
optimal_metrics["mean_anode_potential"],
]
x = np.arange(len(metrics_names))
width = 0.35
axes[1, 1].bar(
x - width / 2,
baseline_values,
width,
label="Baseline",
alpha=0.8,
color="lightblue",
)
axes[1, 1].bar(
x + width / 2, optimal_values, width, label="Optimal", alpha=0.8, color="orange"
)
for i, (baseline, optimal) in enumerate(
zip(baseline_values, optimal_values, strict=False)
):
y_offset = baseline * 0.02 if baseline > 0 else 0.005
axes[1, 1].text(
i - width / 2,
baseline + y_offset,
f"{baseline:.3f}",
ha="center",
va="bottom",
fontsize=9,
)
axes[1, 1].text(
i + width / 2,
optimal + y_offset,
f"{optimal:.3f}",
ha="center",
va="bottom",
fontsize=9,
)
axes[1, 1].set_xlabel("Performance Metrics")
axes[1, 1].set_ylabel("Value")
axes[1, 1].set_title("Performance Metrics Comparison")
axes[1, 1].set_xticks(x)
axes[1, 1].set_xticklabels(metrics_names)
axes[1, 1].set_ylim(top=3.0)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
# Parameter comparison bar chart
param_names = list(parameters.keys())
param_labels = [
"Pos. Active\nMaterial [-]",
"Neg. Active\nMaterial [-]",
]
baseline_param_values = [param.initial_value for param in parameters.values()]
optimal_param_values = [results.parameter_values[name] for name in param_names]
x_params = np.arange(len(param_names))
width_params = 0.35
axes[1, 2].bar(
x_params - width_params / 2,
baseline_param_values,
width_params,
label="Baseline",
alpha=0.8,
color="lightblue",
)
axes[1, 2].bar(
x_params + width_params / 2,
optimal_param_values,
width_params,
label="Optimal",
alpha=0.8,
color="orange",
)
# Add value labels on bars
for i, (baseline, optimal) in enumerate(
zip(baseline_param_values, optimal_param_values, strict=False)
):
y_offset_baseline = baseline * 0.02 if baseline > 0 else 0.005
y_offset_optimal = optimal * 0.02 if optimal > 0 else 0.005
# Format based on parameter type
if i < 2: # Volume fractions
baseline_text = f"{baseline:.3f}"
optimal_text = f"{optimal:.3f}"
elif i < 4: # Thicknesses (now in um)
baseline_text = f"{baseline:.1f}"
optimal_text = f"{optimal:.1f}"
else: # Number of electrodes
baseline_text = f"{baseline:.0f}"
optimal_text = f"{optimal:.0f}"
axes[1, 2].text(
i - width_params / 2,
baseline + y_offset_baseline,
baseline_text,
ha="center",
va="bottom",
fontsize=8,
)
axes[1, 2].text(
i + width_params / 2,
optimal + y_offset_optimal,
optimal_text,
ha="center",
va="bottom",
fontsize=8,
)
axes[1, 2].set_xlabel("Optimization Parameters")
axes[1, 2].set_ylabel("Parameter Value")
axes[1, 2].set_title("Parameter Values: Baseline vs Optimal")
axes[1, 2].set_xticks(x_params)
axes[1, 2].set_ylim(top=1.0)
axes[1, 2].set_xticklabels(param_labels, fontsize=9)
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
# Add upper right subplot for summary statistics
param_changes = []
for _i, (name, param) in enumerate(parameters.items()):
baseline_val = param.initial_value
optimal_val = results.parameter_values[name]
percent_change = ((optimal_val - baseline_val) / baseline_val) * 100
param_changes.append(percent_change)
axes[0, 2].barh(
range(len(param_names)),
param_changes,
color=["orange" if x > 0 else "lightblue" for x in param_changes],
alpha=0.7,
)
axes[0, 2].set_yticks(range(len(param_names)))
axes[0, 2].set_yticklabels(param_labels, fontsize=9)
axes[0, 2].set_xlabel("Percentage Change [%]")
axes[0, 2].set_title("Parameter Changes from Baseline")
axes[0, 2].set_xlim(0, 60)
axes[0, 2].axvline(x=0, color="black", linestyle="-", alpha=0.3)
axes[0, 2].grid(True, alpha=0.3)
# Add percentage labels
for i, change in enumerate(param_changes):
axes[0, 2].text(
change + (1 if change >= 0 else -1),
i,
f"{change:+.1f}%",
va="center",
ha="left" if change >= 0 else "right",
fontsize=8,
)
plt.tight_layout()
plt.show()
print("Optimization analysis completed!")
Optimization analysis completed!