Skills › Scientific Agent Skills › pybamm › simulate_battery.py
simulate_battery.py
214 lines · skills/pybamm/scripts/simulate_battery.py at 5f2bc03
1"""Run a validated constant-current/rest battery experiment with numerical checks."""
2from __future__ import annotations
4import argparse
5import csv
6import hashlib
7import importlib.metadata
8import json
9import math
10import os
11from pathlib import Path
13# Simulation is local; do not send usage telemetry when the library imports.
14os.environ.setdefault("PYBAMM_DISABLE_TELEMETRY", "true")
17def validate_protocol(protocol):
18 def finite(name, value, low=0, high=math.inf):
19 if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value) or not low < value <= high:
20 raise ValueError(f"{name} must be finite and in ({low}, {high}]")
21 if not isinstance(protocol, dict):
22 raise ValueError("Protocol must be a JSON object")
23 allowed = {"model", "parameter_set", "initial_soc", "temperature_K", "sample_period_s", "steps", "parameter_overrides"}
24 if set(protocol) - allowed:
25 raise ValueError(f"Unsupported protocol fields: {sorted(set(protocol) - allowed)}")
26 if protocol.get("model") not in {"SPM", "DFN"}:
27 raise ValueError("model must be SPM or DFN")
28 if not isinstance(protocol.get("parameter_set"), str) or not protocol["parameter_set"]:
29 raise ValueError("parameter_set is required")
30 finite("initial_soc", protocol.get("initial_soc"), high=1)
31 finite("temperature_K", protocol.get("temperature_K"))
32 finite("sample_period_s", protocol.get("sample_period_s"))
33 steps = protocol.get("steps")
34 if not isinstance(steps, list) or not steps:
35 raise ValueError("steps must be a nonempty list")
36 for index, step in enumerate(steps):
37 if not isinstance(step, dict):
38 raise ValueError(f"step {index}: must be a JSON object")
39 if set(step) - {"kind", "duration_s", "c_rate", "until_voltage_V"}:
40 raise ValueError(f"step {index}: unsupported field")
41 if step.get("kind") not in {"discharge", "charge", "rest"}:
42 raise ValueError(f"step {index}: kind must be discharge, charge or rest")
43 finite(f"step {index} duration_s", step.get("duration_s"))
44 if protocol["sample_period_s"] > step["duration_s"]:
45 raise ValueError("sample_period_s exceeds a step duration")
46 if step["kind"] != "rest":
47 finite(f"step {index} c_rate", step.get("c_rate"))
48 if "until_voltage_V" in step:
49 finite(f"step {index} until_voltage_V", step["until_voltage_V"])
50 elif "c_rate" in step or "until_voltage_V" in step:
51 raise ValueError("Rest steps cannot specify c_rate or voltage cutoff")
52 overrides = protocol.get("parameter_overrides", {})
53 if not isinstance(overrides, dict):
54 raise ValueError("parameter_overrides must be a JSON object")
55 for key, value in overrides.items():
56 if not isinstance(key, str):
57 raise ValueError("Parameter override names must be strings")
58 if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value):
59 raise ValueError(f"Override {key}: only finite numeric existing parameters are supported")
60 return protocol
63def solve_protocol(protocol, mesh_points=20, rtol=1e-6, atol=1e-8):
64 import pybamm
65 import numpy as np
66 validate_protocol(protocol)
67 if not isinstance(mesh_points, int) or mesh_points < 5:
68 raise ValueError("mesh_points must be an integer >=5")
69 if not (0 < rtol < 1 and 0 < atol < 1):
70 raise ValueError("Solver tolerances must be between zero and one")
71 parameters = pybamm.ParameterValues(protocol["parameter_set"])
72 overrides = protocol.get("parameter_overrides", {})
73 unknown = set(overrides) - set(parameters.keys())
74 if unknown:
75 raise ValueError(f"Unknown parameter overrides: {sorted(unknown)}")
76 if {"Ambient temperature [K]", "Initial temperature [K]"} & set(overrides):
77 raise ValueError("Set temperature_K instead of temperature overrides")
78 if "Current function [A]" in overrides:
79 raise ValueError("Set step c_rate instead of Current function [A]; experiment steps replace it")
80 if any(key.startswith("Initial concentration in ") and " electrode [" in key for key in overrides):
81 raise ValueError("Set initial_soc instead of initial electrode concentrations; SOC initialization replaces them")
82 if any(not isinstance(parameters[key], (int, float)) for key in overrides):
83 raise ValueError("Only existing numeric parameters may be overridden; preserve functional parameter laws")
84 parameters.update(overrides)
85 parameters.update({"Ambient temperature [K]": protocol["temperature_K"],
86 "Initial temperature [K]": protocol["temperature_K"]})
87 capacity = float(parameters["Nominal cell capacity [A.h]"])
88 lower, upper = (float(parameters[key]) for key in ("Lower voltage cut-off [V]", "Upper voltage cut-off [V]"))
89 if float(parameters["Number of cells connected in series to make a battery"]) != 1:
90 raise ValueError("Only a single cell is supported; battery and cell voltage differ for series strings")
91 if not (capacity > 0 and 0 < lower < upper):
92 raise ValueError("Invalid nominal capacity or voltage limits in parameter set")
93 steps = []
94 for step in protocol["steps"]:
95 current = 0 if step["kind"] == "rest" else step["c_rate"] * capacity * (1 if step["kind"] == "discharge" else -1)
96 cutoff = step.get("until_voltage_V")
97 if cutoff is not None and not lower <= cutoff <= upper:
98 raise ValueError("Protocol voltage cutoff exceeds the parameter-set limits")
99 steps.append(pybamm.step.current(current, duration=step["duration_s"],
100 termination=f"{cutoff} V" if cutoff is not None else None,
101 period=protocol["sample_period_s"], skip_ok=False))
102 model = {"SPM": pybamm.lithium_ion.SPM, "DFN": pybamm.lithium_ion.DFN}[protocol["model"]](options={"thermal": "isothermal"})
103 mesh = {key: mesh_points for key in ("x_n", "x_s", "x_p", "r_n", "r_p")}
104 simulation = pybamm.Simulation(model, parameter_values=parameters,
105 experiment=pybamm.Experiment([tuple(steps)]), var_pts=mesh,
106 solver=pybamm.IDAKLUSolver(rtol=rtol, atol=atol))
107 solution = simulation.solve(initial_soc=protocol["initial_soc"])
108 segments = [segment for cycle in solution.cycles for segment in cycle.steps]
109 if len(segments) != len(steps) or any(isinstance(s, pybamm.EmptySolution) for s in segments):
110 raise RuntimeError("Protocol ended early or skipped a step; do not interpret as a complete experiment")
111 rows, terminations = [], []
112 for index, segment in enumerate(segments):
113 t = np.asarray(segment["Time [s]"].entries)
114 v = np.asarray(segment["Voltage [V]"].entries)
115 current = np.asarray(segment["Current [A]"].entries)
116 charge = np.asarray(segment["Discharge capacity [A.h]"].entries)
117 requested = protocol["steps"][index]
118 if t[-1] - t[0] < requested["duration_s"] - 1e-5:
119 cutoff = requested.get("until_voltage_V")
120 if cutoff is None or abs(float(v[-1]) - cutoff) > 1e-4:
121 raise RuntimeError(f"Step {index} ended before its duration without reaching the requested cutoff")
122 if not np.isfinite(np.column_stack([t, v, current, charge])).all():
123 raise RuntimeError("Nonfinite simulation output")
124 rows.extend(zip([index] * len(t), t, v, current, charge))
125 terminations.append({"step": index, "kind": protocol["steps"][index]["kind"],
126 "start_s": float(t[0]), "end_s": float(t[-1]), "termination": segment.termination})
127 return {"solution": solution, "segments": segments, "rows": rows, "terminations": terminations,
128 "parameters": simulation.parameter_values.to_json(), "mesh_points": mesh_points, "rtol": rtol, "atol": atol,
129 "nominal_capacity_Ah": capacity, "voltage_limits_V": [lower, upper]}
132def compare_runs(first, second):
133 import numpy as np
134 differences = []
135 for index, (a, b) in enumerate(zip(first["segments"], second["segments"])):
136 duration_a, duration_b = a.t[-1] - a.t[0], b.t[-1] - b.t[0]
137 local = np.linspace(0, min(duration_a, duration_b), 101)
138 va = a["Voltage [V]"](a.t[0] + local)
139 vb = b["Voltage [V]"](b.t[0] + local)
140 differences.append({"step": index, "max_voltage_difference_V": float(np.max(np.abs(va-vb))),
141 "end_time_difference_s": float(abs(duration_a-duration_b))})
142 return {"max_voltage_difference_V": max(d["max_voltage_difference_V"] for d in differences),
143 "steps": differences}
146def compare_measurement(rows, measurement_path, residual_path):
147 import numpy as np
148 with Path(measurement_path).open(newline="") as handle:
149 reader = csv.DictReader(handle)
150 if reader.fieldnames != ["time_s", "voltage_V", "current_A"]:
151 raise ValueError("Measured CSV columns must be time_s,voltage_V,current_A")
152 measured = np.asarray([[float(r[k]) for k in reader.fieldnames] for r in reader])
153 if measured.ndim != 2 or len(measured) < 2 or not np.isfinite(measured).all() or (np.diff(measured[:, 0]) <= 0).any():
154 raise ValueError("Need finite measurements and strictly increasing times")
155 simulated = np.asarray(rows)
156 if measured[0, 0] < simulated[0, 1] or measured[-1, 0] > simulated[-1, 1]:
157 raise ValueError("Measurement times exceed simulated interval; no extrapolation")
158 predicted = np.interp(measured[:, 0], simulated[:, 1], simulated[:, 2])
159 predicted_current = np.interp(measured[:, 0], simulated[:, 1], simulated[:, 3])
160 residual = predicted - measured[:, 1]
161 with Path(residual_path).open("w", newline="") as handle:
162 writer = csv.writer(handle); writer.writerow(["time_s", "predicted_voltage_V", "measured_voltage_V", "residual_V", "current_difference_A"])
163 writer.writerows(zip(measured[:, 0], predicted, measured[:, 1], residual, predicted_current-measured[:, 2]))
164 return {"samples": len(measured), "rmse_V": float(np.sqrt(np.mean(residual**2))),
165 "mae_V": float(np.mean(np.abs(residual))), "bias_V": float(np.mean(residual)),
166 "current_rmse_A": float(np.sqrt(np.mean((predicted_current-measured[:, 2])**2))),
167 "measurement_sha256": hashlib.sha256(Path(measurement_path).read_bytes()).hexdigest(),
168 "comparison_interval_s": [float(measured[0, 0]), float(measured[-1, 0])]}
171def run(protocol_path, output, measured=None, mesh_points=20):
172 import pybamm
173 protocol_path, output = Path(protocol_path), Path(output)
174 protocol = validate_protocol(json.loads(protocol_path.read_text()))
175 if output.exists():
176 raise FileExistsError(f"Refusing to overwrite {output}")
177 baseline = solve_protocol(protocol, mesh_points)
178 tight = solve_protocol(protocol, mesh_points, rtol=1e-8, atol=1e-10)
179 refined = solve_protocol(protocol, mesh_points * 2, rtol=1e-8, atol=1e-10)
180 output.mkdir(parents=True)
181 for name, result in (("curve", baseline), ("tight-tolerance", tight), ("refined-mesh", refined)):
182 with (output / f"{name}.csv").open("w", newline="") as handle:
183 writer = csv.writer(handle)
184 writer.writerow(["step", "time_s", "voltage_V", "current_A", "net_discharge_capacity_Ah"])
185 writer.writerows(result["rows"])
186 (output / "parameters.json").write_text(json.dumps(baseline["parameters"], indent=2) + "\n")
187 report = {"pybamm_version": pybamm.__version__, "solver": "IDAKLUSolver", "protocol": protocol,
188 "versions": {name: importlib.metadata.version(name) for name in ("pybamm", "pybammsolvers", "numpy", "casadi")},
189 "protocol_sha256": hashlib.sha256(protocol_path.read_bytes()).hexdigest(),
190 "parameter_set_description": pybamm.parameter_sets.get_docstring(protocol["parameter_set"]),
191 "parameter_snapshot_scope": "Base simulation parameters after initial_soc; per-step currents are set by the protocol, not the snapshot Current function [A]",
192 "nominal_capacity_Ah": baseline["nominal_capacity_Ah"], "voltage_limits_V": baseline["voltage_limits_V"],
193 "mesh_points": mesh_points, "baseline_tolerances": {"rtol": 1e-6, "atol": 1e-8},
194 "tight_tolerances": {"rtol": 1e-8, "atol": 1e-10}, "refined_mesh_points": mesh_points * 2,
195 "terminations": baseline["terminations"], "solver_sensitivity": compare_runs(baseline, tight),
196 "mesh_sensitivity": compare_runs(tight, refined)}
197 if measured:
198 report["measurement_comparison"] = compare_measurement(baseline["rows"], measured, output / "measurement-residuals.csv")
199 (output / "report.json").write_text(json.dumps(report, indent=2) + "\n")
200 return report
203def main():
204 parser = argparse.ArgumentParser(description=__doc__)
205 parser.add_argument("protocol", type=Path)
206 parser.add_argument("--output", type=Path, required=True)
207 parser.add_argument("--measured", type=Path)
208 parser.add_argument("--mesh-points", type=int, default=20)
209 args = parser.parse_args()
210 print(json.dumps(run(args.protocol, args.output, args.measured, args.mesh_points), indent=2))
213if __name__ == "__main__":
214 main()