Silicon Design Optimization¶
This example demonstrates how to optimize battery cell design parameters using ionworkspipeline. We optimize active material volume fractions and electrode thicknesses to maximize gravimetric energy density while preventing lithium plating by constraining the anode potential.
import copy
from functools import reduce
import ionworkspipeline as iwp
from ionworkspipeline.data_fits.models.metrics import Minimum, PointBased, Voltage
from ionworkspipeline.data_fits.models.metrics.actions import (
GreaterThan,
LessThan,
Maximize,
)
import matplotlib.pyplot as plt
import numpy as np
import pybamm
from pybamm import Parameter
Volume Change Functions¶
Define volume change functions for LiCoO2, graphite, and silicon electrodes based on Ai2020 model.
def lico2_volume_change_Ai2020(sto):
omega = pybamm.Parameter("Positive electrode partial molar volume [m3.mol-1]")
c_s_max = pybamm.Parameter("Maximum concentration in positive electrode [mol.m-3]")
t_change = omega * c_s_max * sto
return t_change
def graphite_volume_change_Ai2020(sto):
p1 = 145.907
p2 = -681.229
p3 = 1334.442
p4 = -1415.710
p5 = 873.906
p6 = -312.528
p7 = 60.641
p8 = -5.706
p9 = 0.386
p10 = -4.966e-05
t_change = (
p1 * sto**9
+ p2 * sto**8
+ p3 * sto**7
+ p4 * sto**6
+ p5 * sto**5
+ p6 * sto**4
+ p7 * sto**3
+ p8 * sto**2
+ p9 * sto
+ p10
)
return t_change
def silicon_volume_change(sto):
return graphite_volume_change_Ai2020(sto) * 30.0
Model Configuration¶
Configure DFN model with particle mechanics, set up parameters, and add custom variables.
# Define DFN model with composite electrode and swelling mechanics
model = pybamm.lithium_ion.DFN(
options={
"particle phases": ("2", "1"),
"open-circuit potential": (("single", "current sigmoid"), "single"),
"particle mechanics": "swelling only",
}
)
submesh_types = model.default_submesh_types
for domain in [
"negative electrode",
"separator",
"positive electrode",
"positive particle",
"negative particle",
]:
submesh_types[domain] = pybamm.SymbolicUniform1DSubMesh
params_baseline = pybamm.ParameterValues("Chen2020_composite")
params_baseline.update(
{
"Positive electrode active material density [kg.m-3]": 5010,
"Positive electrode carbon-binder density [kg.m-3]": 2100,
"Primary: Negative electrode active material density [kg.m-3]": 2170,
"Secondary: Negative electrode active material density [kg.m-3]": 2260,
"Primary: Negative electrode carbon-binder density [kg.m-3]": 2100,
"Secondary: Negative electrode carbon-binder density [kg.m-3]": 2100,
"Positive current collector density [kg.m-3]": 2700,
"Negative current collector density [kg.m-3]": 8960,
"Electrolyte density [kg.m-3]": 1200,
"Separator density [kg.m-3]": 397,
"Cell thickness constraint [m]": 3.045e-3,
"Upper voltage cut-off [V]": 4.5,
"Lower voltage cut-off [V]": 2.5,
"Primary: Maximum concentration in negative electrode [mol.m-3]": 28700,
"Primary: Initial concentration in negative electrode [mol.m-3]": 23000,
"Primary: Negative particle diffusivity [m2.s-1]": 5.5e-14,
"Secondary: Negative particle diffusivity [m2.s-1]": 1.67e-14,
"Secondary: Initial concentration in negative electrode [mol.m-3]": 277000,
"Secondary: Maximum concentration in negative electrode [mol.m-3]": 278000,
"Secondary: Initial hysteresis state in negative electrode": 1,
"Positive electrode volume change": lico2_volume_change_Ai2020,
"Positive electrode partial molar volume [m3.mol-1]": -7.28e-07,
"Positive electrode Young's modulus [Pa]": 375000000000.0,
"Positive electrode Poisson's ratio": 0.2,
"Positive electrode critical stress [Pa]": 375000000.0,
"Positive electrode LAM constant exponential term": 2.0,
"Positive electrode reference concentration for free of deformation [mol.m-3]": 0.0,
"Primary: Negative electrode volume change": graphite_volume_change_Ai2020,
"Primary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0,
"Primary: Negative electrode partial molar volume [m3.mol-1]": 3.1e-06,
"Primary: Negative electrode Young's modulus [Pa]": 15000000000.0,
"Primary: Negative electrode Poisson's ratio": 0.3,
"Secondary: Negative electrode volume change": silicon_volume_change,
"Secondary: Negative electrode partial molar volume [m3.mol-1]": 3.1e-06,
"Secondary: Negative electrode Young's modulus [Pa]": 15000000000.0,
"Secondary: Negative electrode Poisson's ratio": 0.3,
"Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]"
"": 0.0,
}
)
N = Parameter("Number of electrodes connected in parallel to make a cell")
params_baseline.update(
{
"Full cell positive current collector thickness [m]": Parameter(
"Positive current collector thickness [m]"
)
* (N + 1),
"Full cell negative current collector thickness [m]": Parameter(
"Negative current collector thickness [m]"
)
* (N + 1),
"Full cell negative electrode thickness [m]": Parameter(
"Negative electrode thickness [m]"
)
* N,
"Full cell positive electrode thickness [m]": Parameter(
"Positive electrode thickness [m]"
)
* N,
"Cell thickness [m]": (
Parameter("Full cell negative electrode thickness [m]")
+ Parameter("Full cell positive electrode thickness [m]")
+ Parameter("Full cell negative current collector thickness [m]")
+ Parameter("Full cell positive current collector thickness [m]")
+ Parameter("Separator thickness [m]") * N
),
}
)
params_baseline["Positive electrode porosity"] = 1.0 - Parameter(
"Positive electrode active material volume fraction"
)
params_baseline["Negative electrode porosity"] = (
1.0
- Parameter("Primary: Negative electrode active material volume fraction")
- Parameter("Secondary: Negative electrode active material volume fraction")
)
energy_density_var = iwp.models.variables.GravimetricEnergyDensity(
"Gravimetric energy density [W.h.kg-1]", model, params_baseline
)
energy_density_var.add()
model.variables["Cell thickness expansion ratio"] = model.variables[
"Cell thickness change [m]"
] / pybamm.Parameter("Cell thickness [m]")
Optimization Setup¶
Define parameters to optimize, experiment protocol, and constraints.
parameters = {
"Positive electrode active material volume fraction": iwp.Parameter(
"Positive electrode active material volume fraction",
initial_value=0.65,
bounds=(0.5, 0.9),
),
"Primary: Negative electrode active material volume fraction": iwp.Parameter(
"Primary: Negative electrode active material volume fraction",
initial_value=0.65,
bounds=(0.5, 0.85),
),
"Secondary: Negative electrode active material volume fraction": iwp.Parameter(
"Secondary: Negative electrode active material volume fraction",
initial_value=0.02,
bounds=(1e-4, 0.1),
),
"Positive electrode thickness [m]": iwp.Parameter(
"Positive electrode thickness [m]", initial_value=68e-6, bounds=(25e-6, 100e-6)
),
"Negative electrode thickness [m]": iwp.Parameter(
"Negative electrode thickness [m]", initial_value=75e-6, bounds=(30e-6, 105e-6)
),
"Number of electrodes connected in parallel to make a cell": iwp.Parameter(
"Number of electrodes connected in parallel to make a cell",
initial_value=15,
bounds=(1, 50),
),
}
model.variables["Cell thickness [m]"] = pybamm.Parameter("Cell thickness [m]")
print(f"Model configured with {len(parameters)} optimization parameters")
Model configured with 6 optimization parameters
Define Experiment Protocol and Optimization Metrics¶
experiment = pybamm.Experiment(
[
"Discharge at 60A until 3.0 V",
"Rest for 1 hours",
"Charge at 60A until 4.05 V",
"Charge at 30A until 4.20 V",
"Hold at 4.20 V until C/20",
]
)
# Define optimization actions
actions = {
"Gravimetric energy density": Maximize(
Voltage(variable="Gravimetric energy density [W.h.kg-1]", value=3.001)
),
}
# Constraints define bounds that must be satisfied (using GreaterThan/LessThan actions)
constraints = {
"Negative surface potential": GreaterThan(
Minimum(
"Negative electrode surface potential difference at separator interface [V]"
),
value=0.0,
penalty=150.0,
),
"Cell thickness expansion ratio": GreaterThan(
Minimum("Cell thickness expansion ratio"),
value=-0.04,
penalty=1e4,
),
"Cell thickness constraint": LessThan(
PointBased("Cell thickness [m]"),
value=3.045e-3,
),
}
print("Experiment, metrics, and constraints configured")
Experiment, metrics, and constraints configured
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(
on_failure="ignore",
options={
"t_no_progress": 0.99,
"num_steps_no_progress": 1000,
},
)
objective_options = {
"model": model,
"simulation_kwargs": {
"experiment": experiment,
"submesh_types": submesh_types,
"solver": solver,
},
}
objective = iwp.objectives.DesignObjective(
options=objective_options,
actions=actions,
constraints=constraints,
parameters={"Initial SOC [%]": 100 * initial_soc},
)
optim = iwp.optimizers.Pints(
method="CMAES",
max_iterations=10, # Increase for better results
max_unchanged_iterations=100,
max_unchanged_iterations_threshold=1e-5,
# population_size=32, # Uncomment to allow larger exploration per iteration
log_to_screen=True,
)
design_datafit = iwp.DataFit(
objective,
parameters=parameters,
options={"seed": 0},
optimizer=optim,
# parallel=True, # Disabled for CI
)
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: 9
--------------------------------------------------------------------------------
Iter. Eval. Best Current Time m:s
--------------------------------------------------------------------------------
1 9 -3.620839e+02 -3.620839e+02 0:05
2 18 -3.620839e+02 -3.620839e+02 0:12
3 27 -3.620839e+02 -3.620839e+02 0:18
4 36 -3.620839e+02 -3.620839e+02 0:23
5 45 -3.620839e+02 -3.620839e+02 0:29
10 90 -3.808504e+02 -3.808504e+02 0:58
--------------------------------------------------------------------------------
Maximum number of iterations (10) reached.
Optimization complete
Optimization completed! Final cost: -380.850425
Optimal parameters:
Positive electrode active material volume fraction: 0.896622
Primary: Negative electrode active material volume fraction: 0.750268
Secondary: Negative electrode active material volume fraction: 0.024297
Positive electrode thickness [m]: 0.000034
Negative electrode thickness [m]: 0.000071
Number of electrodes connected in parallel to make a cell: 16.511662
Validation and Analysis¶
def extract_metrics(solution):
sol_discharge = solution.sub_solutions[0]
gravimetric_energy_density = abs(
sol_discharge["Gravimetric energy density [W.h.kg-1]"]()[-1]
)
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 {
"gravimetric_energy_density": gravimetric_energy_density,
"min_anode_potential": min_anode_potential,
"mean_anode_potential": mean_anode_potential,
"solution": solution,
}
params_optimal = params_baseline.copy()
params_optimal.update(results.parameter_values.copy())
params_baseline_original = copy.deepcopy(params_baseline)
for name, param in parameters.items():
params_baseline_original[name] = param.initial_value
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 - 0.01)
sol_optimal = sim_optimal.solve(initial_soc=initial_soc - 0.01)
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¶
energy_density_improvement = (
(
optimal_metrics["gravimetric_energy_density"]
- baseline_metrics["gravimetric_energy_density"]
)
/ baseline_metrics["gravimetric_energy_density"]
) * 100
min_anode_improvement = (
optimal_metrics["min_anode_potential"] - baseline_metrics["min_anode_potential"]
)
print("Performance Comparison:")
print("=" * 50)
print(f"{'Metric':<45} {'Baseline':<12} {'Optimal':<12} {'Improvement'}")
print("-" * 80)
print(
f"{'Gravimetric Energy Density [W.h.kg-1]':<45} {baseline_metrics['gravimetric_energy_density']:<12.4f} {optimal_metrics['gravimetric_energy_density']:<12.4f} {energy_density_improvement:+.2f}%"
)
print(
f"{'Min Anode Potential [V]':<45} {baseline_metrics['min_anode_potential']:<12.4f} {optimal_metrics['min_anode_potential']:<12.4f} {min_anode_improvement:+.4f}"
)
print(
f"{'Mean Anode Potential [V]':<45} {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}%)")
baseline_thickness = (
params_baseline_original.evaluate(params_baseline_original["Cell thickness [m]"])
if "Cell thickness [m]" in params_baseline_original
else "N/A"
)
optimal_thickness = params_optimal.evaluate(params_optimal["Cell thickness [m]"])
print(
f"\nCell Thickness: {baseline_thickness} → {optimal_thickness:.6f} m (limit: {params_baseline['Cell thickness constraint [m]']:.6f} m)"
)
Performance Comparison:
==================================================
Metric Baseline Optimal Improvement
--------------------------------------------------------------------------------
Gravimetric Energy Density [W.h.kg-1] 318.4746 374.3864 +17.56%
Min Anode Potential [V] 0.0319 0.0194 -0.0125
Mean Anode Potential [V] 0.0843 0.0672
==================================================
Parameter Changes:
Positive electrode active material volume fraction: 0.650000 → 0.896622 (+37.94%)
Primary: Negative electrode active material volume fraction: 0.650000 → 0.750268 (+15.43%)
Secondary: Negative electrode active material volume fraction: 0.020000 → 0.024297 (+21.49%)
Positive electrode thickness [m]: 0.000068 → 0.000034 (-49.63%)
Negative electrode thickness [m]: 0.000075 → 0.000071 (-4.98%)
Number of electrodes connected in parallel to make a cell: 15.000000 → 16.511662 (+10.08%)
Cell Thickness: 0.0027730000000000003 → 0.002431 m (limit: 0.003045 m)
Visualization¶
def plot_comparison(ax, t_base, y_base, t_opt, y_opt, xlabel, ylabel, title):
ax.plot(t_base, y_base, label="Baseline", linewidth=2, alpha=0.8)
ax.plot(t_opt, y_opt, label="Optimal", linewidth=2, alpha=0.8, linestyle="--")
ax.set(xlabel=xlabel, ylabel=ylabel, title=title)
ax.legend()
ax.grid(True, alpha=0.3)
def add_bar_labels(ax, x_pos, values, fmt, offset=0.02):
for i, val in enumerate(values):
y_off = val * offset if val > 0 else 0.005
ax.text(
x_pos[i], val + y_off, fmt.format(val), ha="center", va="bottom", fontsize=8
)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# Voltage
sol_base, sol_opt = baseline_metrics["solution"], optimal_metrics["solution"]
plot_comparison(
axes[0, 0],
sol_base.t / 3600,
sol_base["Voltage [V]"](sol_base.t),
sol_opt.t / 3600,
sol_opt["Voltage [V]"](sol_opt.t),
"Time [h]",
"Voltage [V]",
"Cell Voltage: Baseline vs Optimal",
)
# Energy density
discharge_base, discharge_opt = sol_base.sub_solutions[0], sol_opt.sub_solutions[0]
plot_comparison(
axes[0, 1],
discharge_base.t / 3600,
abs(discharge_base["Gravimetric energy density [W.h.kg-1]"](discharge_base.t)),
discharge_opt.t / 3600,
abs(discharge_opt["Gravimetric energy density [W.h.kg-1]"](discharge_opt.t)),
"Time [h]",
"Gravimetric Energy Density [W.h.kg-1]",
"Gravimetric Energy Density",
)
# Anode potential
charge_base = reduce(lambda a, b: a + b, sol_base.sub_solutions[2:])
charge_opt = reduce(lambda a, b: a + b, sol_opt.sub_solutions[2:])
anode_var = "Negative electrode surface potential difference at separator interface [V]"
t_charge_base = (charge_base.t - charge_base.t[0]) / 3600
t_charge_opt = (charge_opt.t - charge_opt.t[0]) / 3600
plot_comparison(
axes[1, 0],
t_charge_base,
charge_base[anode_var](charge_base.t),
t_charge_opt,
charge_opt[anode_var](charge_opt.t),
"Charge Time [h]",
"Anode - Electrolyte Potential [V]",
"Anode Potential During Charge",
)
axes[1, 0].axhline(0, color="r", linestyle=":", alpha=0.7, label="Safety limit")
# Performance metrics
metrics_data = [
(
[
"Gravimetric Energy\nDensity [W.h.kg-1]",
"Min Anode\nPotential [V]",
"Mean Anode\nPotential [V]",
],
[
baseline_metrics["gravimetric_energy_density"],
baseline_metrics["min_anode_potential"],
baseline_metrics["mean_anode_potential"],
],
[
optimal_metrics["gravimetric_energy_density"],
optimal_metrics["min_anode_potential"],
optimal_metrics["mean_anode_potential"],
],
)
]
x = np.arange(len(metrics_data[0][0]))
width = 0.35
axes[1, 1].bar(
x - width / 2,
metrics_data[0][1],
width,
label="Baseline",
alpha=0.8,
color="lightblue",
)
axes[1, 1].bar(
x + width / 2, metrics_data[0][2], width, label="Optimal", alpha=0.8, color="orange"
)
add_bar_labels(axes[1, 1], x - width / 2, metrics_data[0][1], "{:.3f}")
add_bar_labels(axes[1, 1], x + width / 2, metrics_data[0][2], "{:.3f}")
axes[1, 1].set(
xlabel="Performance Metrics",
ylabel="Value",
title="Performance Metrics Comparison",
xticks=x,
)
axes[1, 1].set_xticklabels(metrics_data[0][0])
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
# Parameter values
param_labels = [
"Pos. Active\nMaterial [-]",
"Primary: Neg. Active\nMaterial [-]",
"Secondary: Neg. Active\nMaterial [-]",
"Pos. Thickness\n[um]",
"Neg. Thickness\n[um]",
"Number of\nElectrodes [-]",
]
base_vals = [p.initial_value for p in parameters.values()]
opt_vals = [results.parameter_values[n] for n in parameters.keys()]
for i in [3, 4]:
base_vals[i] *= 1e6
opt_vals[i] *= 1e6
x_params = np.arange(len(param_labels))
axes[1, 2].bar(
x_params - width / 2,
base_vals,
width,
label="Baseline",
alpha=0.8,
color="lightblue",
)
axes[1, 2].bar(
x_params + width / 2, opt_vals, width, label="Optimal", alpha=0.8, color="orange"
)
for i, (b, o) in enumerate(zip(base_vals, opt_vals, strict=False)):
fmt = "{:.3f}" if i < 2 else "{:.1f}" if i < 4 else "{:.0f}"
add_bar_labels(axes[1, 2], [x_params[i] - width / 2], [b], fmt)
add_bar_labels(axes[1, 2], [x_params[i] + width / 2], [o], fmt)
axes[1, 2].set(
xlabel="Optimization Parameters",
ylabel="Parameter Value",
title="Parameter Values: Baseline vs Optimal",
xticks=x_params,
ylim=(None, 85),
)
axes[1, 2].set_xticklabels(param_labels, fontsize=9)
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
# Parameter changes
changes = [
100 * (results.parameter_values[n] - p.initial_value) / p.initial_value
for n, p in parameters.items()
]
axes[0, 2].barh(
range(len(changes)),
changes,
alpha=0.7,
color=["orange" if c > 0 else "lightblue" for c in changes],
)
axes[0, 2].set(
xlabel="Percentage Change [%]",
title="Parameter Changes from Baseline",
yticks=range(len(param_labels)),
xlim=(-100, 110),
)
axes[0, 2].set_yticklabels(param_labels, fontsize=9)
axes[0, 2].axvline(0, color="black", linestyle="-", alpha=0.3)
axes[0, 2].grid(True, alpha=0.3)
for i, c in enumerate(changes):
axes[0, 2].text(
c + (1 if c >= 0 else -1),
i,
f"{c:+.1f}%",
va="center",
ha="left" if c >= 0 else "right",
fontsize=8,
)
plt.tight_layout()
plt.show()
print("Optimization analysis completed!")
Optimization analysis completed!
Verify expansion ratio¶
expansion_ratio = np.ptp(
sol_optimal["Cell thickness change [m]"].data
) / params_optimal.evaluate(params_optimal["Cell thickness [m]"])
print(f"Cell thickness expansion percentage: {expansion_ratio * 100.0:.4f}%")
Cell thickness expansion percentage: 3.9028%