diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 000000000..7ef4b7759 --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,31 @@ +name: JossPaperCompilation + +on: + push: + paths: + - papers/joss/** + pull_request: + paths: + - 'papers/joss/**' + workflow_dispatch: +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: papers/joss/paper.md + - name: Upload + uses: actions/upload-artifact@v4 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: papers/joss/paper.pdf diff --git a/.gitignore b/.gitignore index b2123edf4..a2b444c1d 100644 --- a/.gitignore +++ b/.gitignore @@ -323,3 +323,6 @@ results/ # Pycharm *.idea/ + +# JOSS +jats/ diff --git a/CITATION.cff b/CITATION.cff index 1f7be5c7a..fcadd20ac 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,35 +1,62 @@ -cff-version: 1.2.0 -title: 'PyBOP: A Python package for battery model optimisation and parameterisation' -message: >- - If you use this software, please cite the article below. +cff-version: "1.2.0" authors: - - given-names: Brady - family-names: Planden +- family-names: Planden + given-names: Brady + orcid: "https://orcid.org/0000-0002-1082-9125" +- family-names: Courtier + given-names: Nicola E. + orcid: "https://orcid.org/0000-0002-5714-1096" +- family-names: Robinson + given-names: Martin + orcid: "https://orcid.org/0000-0002-1572-6782" +- family-names: Khetarpal + given-names: Agriya + orcid: "https://orcid.org/0000-0002-1112-1786" +- family-names: Planella + given-names: Ferran Brosa + orcid: "https://orcid.org/0000-0001-6363-2812" +- family-names: Howey + given-names: David A. + orcid: "https://orcid.org/0000-0002-0620-3955" +contact: +- family-names: Howey + given-names: David A. + orcid: "https://orcid.org/0000-0002-0620-3955" +doi: 10.5281/zenodo.17711656 +message: If you use this software, please cite our article in the + Journal of Open Source Software. +preferred-citation: + authors: + - family-names: Planden + given-names: Brady orcid: "https://orcid.org/0000-0002-1082-9125" - - given-names: Nicola E. - family-names: Courtier + - family-names: Courtier + given-names: Nicola E. orcid: "https://orcid.org/0000-0002-5714-1096" - - given-names: Martin - family-names: Robinson + - family-names: Robinson + given-names: Martin orcid: "https://orcid.org/0000-0002-1572-6782" - - given-names: Agriya - family-names: Khetarpal + - family-names: Khetarpal + given-names: Agriya orcid: "https://orcid.org/0000-0002-1112-1786" - - given-names: Ferran - family-names: Brosa Planella + - family-names: Planella + given-names: Ferran Brosa orcid: "https://orcid.org/0000-0001-6363-2812" - - given-names: David A. - family-names: Howey + - family-names: Howey + given-names: David A. orcid: "https://orcid.org/0000-0002-0620-3955" - -keywords: - - "python" - - "battery models" - - "parameter inference" - - "optimization" - -journal: "arXiv" -date-released: 2024-12-20 -doi: 10.48550/arXiv.2412.15859 -version: "25.11" # Update this alongside new releases -repository-code: 'https://www.github.com/pybop-team/pybop' + date-published: 2025-12-03 + doi: 10.21105/joss.07874 + issn: 2475-9066 + issue: 116 + journal: Journal of Open Source Software + publisher: + name: Open Journals + start: 7874 + title: "PyBOP: A Python package for battery model optimisation and + parameterisation" + type: article + url: "https://joss.theoj.org/papers/10.21105/joss.07874" + volume: 10 +title: "PyBOP: A Python package for battery model optimisation and + parameterisation" diff --git a/README.md b/README.md index 42a46bb7d..58c94e269 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,7 @@ [![nbviewer](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)](https://nbviewer.org/github/pybop-team/PyBOP/tree/develop/examples/notebooks/) [![Static Badge](https://img.shields.io/badge/https%3A%2F%2Fpybop-team.github.io%2Fpybop-bench%2F?label=Benchmarks)](https://pybop-team.github.io/pybop-bench/) [![Releases](https://img.shields.io/github/v/release/pybop-team/PyBOP?color=gold)](https://github.com/pybop-team/PyBOP/releases) + [![DOI](https://joss.theoj.org/papers/10.21105/joss.07874/status.svg)](https://doi.org/10.21105/joss.07874) [Main Branch Examples](https://github.com/pybop-team/PyBOP/tree/main/examples) [Develop Branch Examples](https://github.com/pybop-team/PyBOP/tree/develop/examples) diff --git a/papers/joss/design_plots.py b/papers/joss/design_plots.py new file mode 100644 index 000000000..c7aec3231 --- /dev/null +++ b/papers/joss/design_plots.py @@ -0,0 +1,164 @@ +# A script to generate design optimisation plots for the JOSS paper. + + +import numpy as np +import pybamm +from pybamm import Parameter + +import pybop +from pybop.plot import PlotlyManager + +go = PlotlyManager().go +np.random.seed(8) +axis_font_size = 24 +tick_font_size = 20 + +# Choose which plots to show and save +create_plot = {} +create_plot["gravimetric"] = True # takes longer +create_plot["prediction"] = True + +# Define model and parameter values +model = pybamm.lithium_ion.SPMe() +pybop.pybamm.add_variable_to_model(model, "Gravimetric energy density [W.h.kg-1]") +parameter_values = pybamm.ParameterValues("Chen2020") +pybop.pybamm.set_formation_concentrations(parameter_values) +parameter_values.update( + { + "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), + "Negative electrode active material density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode porosity": 1.0 + - Parameter("Positive electrode active material volume fraction"), + "Cell mass [kg]": pybop.pybamm.cell_mass(), + }, + check_already_exists=False, +) + +# Fitting parameters +parameter_values.update( + { + "Positive electrode thickness [m]": pybop.Parameter( + initial_value=8.88e-05, + prior=pybop.Gaussian(7.56e-05, 3e-05), + bounds=[50e-06, 120e-06], + transformation=pybop.UnitHyperCube(lower=50e-6, upper=120e-6), + ), + "Positive electrode active material volume fraction": pybop.Parameter( + initial_value=0.42, + prior=pybop.Gaussian(0.58, 0.1), + bounds=[0.3, 0.825], + transformation=pybop.UnitHyperCube(lower=0.3, upper=0.825), + ), + } +) + +# Define test protocol +experiment = pybamm.Experiment(["Discharge at 1C until 2.5 V (10 second period)"]) +Q = parameter_values["Nominal cell capacity [A.h]"] +print(f"The 1C rate is {Q} A.") + +# Generate problem +simulator = pybop.pybamm.Simulator( + model, + parameter_values=parameter_values, + protocol=experiment, + solver=pybamm.IDAKLUSolver(atol=1e-6, rtol=1e-6), +) +cost = pybop.DesignCost(target="Gravimetric energy density [W.h.kg-1]") +problem = pybop.Problem(simulator, cost) + +# Set up the optimiser +options = pybop.PintsOptions(max_iterations=250, max_unchanged_iterations=45) +optim = pybop.NelderMead(problem, options=options) + +# Run optimisation +result = optim.run() +print(result) +print("Estimated parameters:", result.x) +print(f"Initial gravimetric energy density: {problem(result.x0):.1f} W.h.kg-1") +print(f"Optimised gravimetric energy density: {problem(result.x):.1f} W.h.kg-1") + +if create_plot["gravimetric"]: + # Plot the cost landscape with optimisation path + gravimetric_fig = pybop.plot.contour( + result, + steps=25, + show=False, + xaxis=dict( + title=dict( + text="Positive electrode thickness / m", font_size=axis_font_size + ), + tickfont_size=tick_font_size, + exponentformat="power", + ), + yaxis=dict(title_font_size=axis_font_size, tickfont_size=tick_font_size), + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font_size=tick_font_size, + ), + coloraxis_colorbar=dict(tickfont_size=tick_font_size), + margin=dict(t=50), + title=None, + ) + gravimetric_fig.write_image("figures/individual/design_gravimetric.pdf") + +if create_plot["prediction"]: + # Plot the timeseries output + problem.target = "Voltage [V]" + figs = pybop.plot.problem( + problem, + inputs=result.best_inputs, + title=None, + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font_size=tick_font_size - 2, + ), + xaxis=dict( + showline=True, + linewidth=1, + linecolor="black", + mirror=True, + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + ), + yaxis=dict( + showline=True, + linewidth=1, + linecolor="black", + mirror=True, + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + ), + margin=dict(t=60, b=84, r=50, l=15), + show=False, + ) + + prediction_fig = figs[0] + prediction_fig.data[1].update(line=dict(color="#00CC97")) + prediction_fig.data[ + 0 + ].name = f"Initial: {problem(result.x0):.1f} W h kg-1" + prediction_fig.data[ + 1 + ].name = f"Optimised: {problem(result.x):.1f} W h kg-1" + prediction_fig.show() + prediction_fig.write_image("figures/individual/design_prediction.pdf") diff --git a/papers/joss/figures/PyBOP_components.drawio.pdf b/papers/joss/figures/PyBOP_components.drawio.pdf new file mode 100644 index 000000000..efcb9074b Binary files /dev/null and b/papers/joss/figures/PyBOP_components.drawio.pdf differ diff --git a/papers/joss/figures/combined/contour_subplot.pdf b/papers/joss/figures/combined/contour_subplot.pdf new file mode 100644 index 000000000..6cabccc5a Binary files /dev/null and b/papers/joss/figures/combined/contour_subplot.pdf differ diff --git a/papers/joss/figures/combined/converge.pdf b/papers/joss/figures/combined/converge.pdf new file mode 100644 index 000000000..80d7cbe0c Binary files /dev/null and b/papers/joss/figures/combined/converge.pdf differ diff --git a/papers/joss/figures/combined/design.pdf b/papers/joss/figures/combined/design.pdf new file mode 100644 index 000000000..988e7a4d3 Binary files /dev/null and b/papers/joss/figures/combined/design.pdf differ diff --git a/papers/joss/figures/combined/impedance.pdf b/papers/joss/figures/combined/impedance.pdf new file mode 100644 index 000000000..1d9977841 Binary files /dev/null and b/papers/joss/figures/combined/impedance.pdf differ diff --git a/papers/joss/figures/combined/optimisers_parameters.pdf b/papers/joss/figures/combined/optimisers_parameters.pdf new file mode 100644 index 000000000..81c5e3faa Binary files /dev/null and b/papers/joss/figures/combined/optimisers_parameters.pdf differ diff --git a/papers/joss/figures/combined/posteriors.pdf b/papers/joss/figures/combined/posteriors.pdf new file mode 100644 index 000000000..abae4216e Binary files /dev/null and b/papers/joss/figures/combined/posteriors.pdf differ diff --git a/papers/joss/figures/combined/sim-landscape.pdf b/papers/joss/figures/combined/sim-landscape.pdf new file mode 100644 index 000000000..325ea0bf7 Binary files /dev/null and b/papers/joss/figures/combined/sim-landscape.pdf differ diff --git a/papers/joss/figures/individual/convergence_maximising.pdf b/papers/joss/figures/individual/convergence_maximising.pdf new file mode 100644 index 000000000..5f9f986fa Binary files /dev/null and b/papers/joss/figures/individual/convergence_maximising.pdf differ diff --git a/papers/joss/figures/individual/convergence_minimising.pdf b/papers/joss/figures/individual/convergence_minimising.pdf new file mode 100644 index 000000000..641018e6f Binary files /dev/null and b/papers/joss/figures/individual/convergence_minimising.pdf differ diff --git a/papers/joss/figures/individual/design_gravimetric.pdf b/papers/joss/figures/individual/design_gravimetric.pdf new file mode 100644 index 000000000..de0dddbe7 Binary files /dev/null and b/papers/joss/figures/individual/design_gravimetric.pdf differ diff --git a/papers/joss/figures/individual/design_prediction.pdf b/papers/joss/figures/individual/design_prediction.pdf new file mode 100644 index 000000000..15237a3b4 Binary files /dev/null and b/papers/joss/figures/individual/design_prediction.pdf differ diff --git a/papers/joss/figures/individual/evolution_parameters.pdf b/papers/joss/figures/individual/evolution_parameters.pdf new file mode 100644 index 000000000..1676ddcb5 Binary files /dev/null and b/papers/joss/figures/individual/evolution_parameters.pdf differ diff --git a/papers/joss/figures/individual/gradient_parameters.pdf b/papers/joss/figures/individual/gradient_parameters.pdf new file mode 100644 index 000000000..7160aa0b4 Binary files /dev/null and b/papers/joss/figures/individual/gradient_parameters.pdf differ diff --git a/papers/joss/figures/individual/heuristic_parameters.pdf b/papers/joss/figures/individual/heuristic_parameters.pdf new file mode 100644 index 000000000..fae3e3a7d Binary files /dev/null and b/papers/joss/figures/individual/heuristic_parameters.pdf differ diff --git a/papers/joss/figures/individual/impedance_contour.pdf b/papers/joss/figures/individual/impedance_contour.pdf new file mode 100644 index 000000000..e70eddc86 Binary files /dev/null and b/papers/joss/figures/individual/impedance_contour.pdf differ diff --git a/papers/joss/figures/individual/impedance_spectrum.pdf b/papers/joss/figures/individual/impedance_spectrum.pdf new file mode 100644 index 000000000..61279ff52 Binary files /dev/null and b/papers/joss/figures/individual/impedance_spectrum.pdf differ diff --git a/papers/joss/figures/individual/landscape.pdf b/papers/joss/figures/individual/landscape.pdf new file mode 100644 index 000000000..761613af3 Binary files /dev/null and b/papers/joss/figures/individual/landscape.pdf differ diff --git a/papers/joss/figures/individual/simulation.pdf b/papers/joss/figures/individual/simulation.pdf new file mode 100644 index 000000000..8e8f857b6 Binary files /dev/null and b/papers/joss/figures/individual/simulation.pdf differ diff --git a/papers/joss/image_combine.sh b/papers/joss/image_combine.sh new file mode 100755 index 000000000..1e8a27c93 --- /dev/null +++ b/papers/joss/image_combine.sh @@ -0,0 +1,9 @@ +pdfjam --nup 2x1 figures/individual/design_prediction.pdf figures/individual/design_gravimetric.pdf --landscape --outfile figures/combined/design.pdf --papersize '{420px,900px}' + +pdfjam --nup 3x1 figures/individual/gradient_parameters.pdf figures/individual/evolution_parameters.pdf figures/individual/heuristic_parameters.pdf --landscape --outfile figures/combined/optimisers_parameters.pdf --papersize '{900px,1440px}' + +pdfjam --nup 2x1 figures/individual/convergence_minimising.pdf figures/individual/convergence_maximising.pdf --landscape --outfile figures/combined/converge.pdf --papersize '{400px,900px}' --trim '0cm 0cm 0cm 1.5cm' --clip true + +pdfjam --nup 2x1 figures/individual/simulation.pdf figures/individual/landscape.pdf --landscape --outfile figures/combined/sim-landscape.pdf --papersize '{420px,900px}' + +pdfjam --nup 2x1 figures/individual/impedance_spectrum.pdf figures/individual/impedance_contour.pdf --landscape --outfile figures/combined/impedance.pdf --papersize '{420px,900px}' diff --git a/papers/joss/paper.bib b/papers/joss/paper.bib new file mode 100644 index 000000000..2ff12e634 --- /dev/null +++ b/papers/joss/paper.bib @@ -0,0 +1,380 @@ + +@article{Sulzer:2021, + title = {{Python Battery Mathematical Modelling (PyBaMM)}}, + author = {Sulzer, Valentin and Marquis, Scott G. and Timms, Robert and + Robinson, Martin and Chapman, S. Jon}, + doi = {10.5334/jors.309}, + journal = {Journal of Open Research Software}, + publisher = {Software Sustainability Institute}, + volume = {9}, + number = {1}, + pages = {14}, + year = {2021}, + url = {https://github.com/pybamm-team/PyBaMM}, +} + +@article{Clerx:2019, + title={Probabilistic Inference on Noisy Time Series ({PINTS})}, + author={Clerx, Michael and Robinson, Martin and Lambert, Ben and Lei, Chon Lok and + Ghosh, Sanmitra and Mirams, Gary R and Gavaghan, David J}, + journal={Journal of Open Research Software}, + volume={7}, + number={1}, + pages={23}, + year={2019}, + doi={10.5334/jors.252} +} + +@article{SciPy:2020, + author = {Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E. and + Haberland, Matt and Reddy, Tyler and Cournapeau, David and + Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and + Bright, Jonathan and {van der Walt}, St{\'e}fan J. and + Brett, Matthew and Wilson, Joshua and Millman, K. Jarrod and + Mayorov, Nikolay and Nelson, Andrew R. J. and Jones, Eric and + Kern, Robert and Larson, Eric and Carey, C J and + Polat, {\.I}lhan and Feng, Yu and Moore, Eric W. and + {VanderPlas}, Jake and Laxalde, Denis and Perktold, Josef and + Cimrman, Robert and Henriksen, Ian and Quintero, E. A. and + Harris, Charles R. and Archibald, Anne M. and + Ribeiro, Ant{\^o}nio H. and Pedregosa, Fabian and + {van Mulbregt}, Paul and {SciPy 1.0 Contributors}}, + title = {{{SciPy} 1.0: Fundamental Algorithms for Scientific + Computing in Python}}, + journal = {Nature Methods}, + year = {2020}, + volume = {17}, + pages = {261--272}, + url = {https://doi.org/10.1038/s41592-019-0686-2}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2020NatMe..17..261V}, + doi = {10.1038/s41592-019-0686-2}, +} + +@misc{BPX:2023, + title = {Battery Parameter eXchange}, + author = {Korotkin, Ivan and Timms, Robert and Foster, Jamie Foster and + Dickinson, Edmund and Robinson, Martin}, + publisher = {The Faraday Institution}, + year = {2023}, + journal = {GitHub repository}, + url = {https://github.com/FaradayInstitution/BPX}, +} + +@article{Doyle:1993, + author = {Doyle, Marc and Fuller, Thomas F. and Newman, John}, + doi = {10.1149/1.2221597}, + issn = {0013-4651}, + journal = {Journal of The Electrochemical Society}, + number = {6}, + pages = {1526--1533}, + title = {{Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell}}, + volume = {140}, + year = {1993} +} + +@article{Fuller:1994, + doi = {10.1149/1.2054684}, + url = {https://dx.doi.org/10.1149/1.2054684}, + year = {1994}, + month = {jan}, + publisher = {The Electrochemical Society, Inc.}, + volume = {141}, + number = {1}, + pages = {1}, + author = {Thomas F. Fuller and Marc Doyle and John Newman}, + title = {Simulation and Optimization of the Dual Lithium Ion Insertion Cell}, + journal = {Journal of The Electrochemical Society}, +} + +@article{Planella:2022, + author = {Brosa Planella, Ferran and Ai, Weilong and Boyce, Adam M and + Ghosh, Abir and Korotkin, Ivan and Sahu, Smita and Sulzer, Valentin and + Timms, Robert and Tranter, Thomas G and Zyskin, Maxim and + Cooper, Samuel J and Edge, Jacqueline S and Foster, Jamie M and + Marinescu, Monica and Wu, Billy and Richardson, Giles}, + doi = {10.1088/2516-1083/ac7d31}, + issn = {2516-1083}, + journal = {Progress in Energy}, + month = {oct}, + number = {4}, + pages = {042003}, + title = {{A Continuum of Physics-Based Lithium-Ion Battery Models Reviewed}}, + url = {https://iopscience.iop.org/article/10.1088/2516-1083/ac7d31}, + volume = {4}, + year = {2022} +} + +@article{Verbrugge:2017, + doi = {10.1149/2.0341708jes}, + url = {https://dx.doi.org/10.1149/2.0341708jes}, + year = {2017}, + month = {may}, + publisher = {The Electrochemical Society}, + volume = {164}, + number = {11}, + pages = {E3243}, + author = {Mark Verbrugge and Daniel Baker and Brian Koch and Xingcheng Xiao and Wentian Gu}, + title = {Thermodynamic Model for Substitutional Materials: Application to Lithiated Graphite, Spinel Manganese Oxide, Iron Phosphate, and Layered Nickel-Manganese-Cobalt Oxide}, + journal = {Journal of The Electrochemical Society}, +} + +@article{Tranter2022, + doi = {10.21105/joss.04051}, + url = {https://doi.org/10.21105/joss.04051}, + year = {2022}, + publisher = {The Open Journal}, + volume = {7}, + number = {70}, + pages = {4051}, + author = {Thomas G. Tranter and Robert Timms and Valentin Sulzer and + Ferran Brosa Planella and Gavin M. Wiggins and Suryanarayana V. Karra and + Priyanshu Agarwal and Saransh Chopra and Srikanth Allu and + Paul R. Shearing and Dan J. l. Brett}, + title = {{liionpack: A Python package for simulating packs of batteries with PyBaMM}}, + journal = {Journal of Open Source Software} +} + +@article{Wang:2022, + author = {Wang, A. A. and O'Kane, S. E. J. and Brosa Planella, F. and + Houx, J. Le and O'Regan, K. and Zyskin, M. and Edge, J. and + Monroe, C. W. and Cooper, S. J. and Howey, D. A. and Kendrick, E. and + Foster, J. M.}, + journal = {Progress in Energy}, + title = {Review of parameterisation and a novel database {(LiionDB)} for continuum {Li-ion} battery models}, + year = {2022}, + issn = {2516-1083}, + month = {jul}, + number = {3}, + pages = {032004}, + volume = {4}, + doi = {10.1088/2516-1083/ac692c}, + url = {https://iopscience.iop.org/article/10.1088/2516-1083/ac692c}, +} + +@article{Andersson:2022, + author = {Andersson, Malin and Streb, Moritz and Ko, Jing Ying and + {L{\"{o}}fqvist Klass}, Verena and Klett, Matilda and + Ekstr{\"{o}}m, Henrik and Johansson, Mikael and Lindbergh, G{\"{o}}ran}, + journal = {Journal of Power Sources}, + title = {Parametrization of physics-based battery models from input-output data: {A} review of methodology and current research}, + year = {2022}, + issn = {03787753}, + month = {feb}, + number = {November 2021}, + pages = {230859}, + volume = {521}, + doi = {10.1016/j.jpowsour.2021.230859}, + publisher = {Elsevier B.V.}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0378775321013458}, +} + +@article{Miguel:2021, + author = {Miguel, E. and Plett, Gregory L. and Trimboli, M. Scott and + Oca, L. and Iraola, U. and Bekaert, E.}, + journal = {Journal of Energy Storage}, + title = {Review of computational parameter estimation methods for electrochemical models}, + year = {2021}, + issn = {2352152X}, + number = {PB}, + pages = {103388}, + volume = {44}, + doi = {10.1016/j.est.2021.103388}, + publisher = {Elsevier Ltd}, + url = {https://doi.org/10.1016/j.est.2021.103388}, +} + +@article{Chu:2019, + author = {Chu, Zhengyu and Plett, Gregory L. and Trimboli, M. Scott and Ouyang, Minggao}, + journal = {Journal of Energy Storage}, + title = {A control-oriented electrochemical model for lithium-ion battery, {P}art {I}: {L}umped-parameter reduced-order model with constant phase element}, + year = {2019}, + issn = {2352152X}, + number = {August}, + pages = {100828}, + volume = {25}, + doi = {10.1016/j.est.2019.100828}, + publisher = {Elsevier}, + url = {https://doi.org/10.1016/j.est.2019.100828}, +} + +@article{Chen:2020, + author = {Chen, Chang-Hui and {Brosa Planella}, Ferran and O'Regan, Kieran and + Gastol, Dominika and Widanage, W. Dhammika and Kendrick, Emma}, + journal = {Journal of The Electrochemical Society}, + title = {Development of experimental techniques for parameterization of multi-scale lithium-ion battery models}, + year = {2020}, + issn = {0013-4651}, + month = {jan}, + number = {8}, + pages = {080534}, + volume = {167}, + doi = {10.1149/1945-7111/ab9050}, + publisher = {IOP Publishing}, + url = {https://iopscience.iop.org/article/10.1149/1945-7111/ab9050}, +} + +@article{Kirk:2022, + author = {Kirk, Toby L. and Lewis-Douglas, Adam and Howey, David and Please, Colin P. and {Jon Chapman}, S.}, + journal = {Journal of The Electrochemical Society}, + title = {Nonlinear electrochemical impedance spectroscopy for lithium-ion battery model parameterization}, + year = {2023}, + issn = {0013-4651}, + month = {jan}, + number = {1}, + pages = {010514}, + volume = {170}, + doi = {10.1149/1945-7111/acada7}, + url = {https://iopscience.iop.org/article/10.1149/1945-7111/acada7}, +} + +@article{Couto:2023, + author = {Couto, Luis D. and Charkhgard, Mohammad and Karaman, Berke and + Job, Nathalie and Kinnaert, Michel}, + journal = {Energy}, + title = {Lithium-ion battery design optimization based on a dimensionless reduced-order electrochemical model}, + year = {2023}, + issn = {03605442}, + number = {PE}, + pages = {125966}, + volume = {263}, + doi = {10.1016/j.energy.2022.125966}, + publisher = {Elsevier Ltd}, + url = {https://doi.org/10.1016/j.energy.2022.125966}, +} + +@misc{NUTS:2011, + title={{The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo}}, + author={Matthew D. Hoffman and Andrew Gelman}, + year={2011}, + eprint={1111.4246}, + archivePrefix={arXiv}, + primaryClass={stat.CO}, + url={https://arxiv.org/abs/1111.4246}, +} + +@book{Hamiltonian:2011, + title={Handbook of Markov Chain Monte Carlo}, + ISBN={9780429138508}, + url={http://dx.doi.org/10.1201/b10905}, + DOI={10.1201/b10905}, + publisher={Chapman and Hall/CRC}, + author={Brooks, Steve and Gelman, Andrew and Jones, Galin and Meng, Xiao-Li}, + year={2011}, + month=may, +} + +@article{Haario:2001, + title = {An {Adaptive} {Metropolis} {Algorithm}}, + volume = {7}, + issn = {13507265}, + url = {https://www.jstor.org/stable/3318737?origin=crossref}, + doi = {10.2307/3318737}, + number = {2}, + urldate = {2024-11-22}, + journal = {Bernoulli}, + author = {Haario, Heikki and Saksman, Eero and Tamminen, Johanna}, + month = apr, + year = {2001}, + pages = {223}, +} + +@article{DiffEvolution:2006, + title = {A {Markov} {Chain} {Monte} {Carlo} version of the genetic algorithm {Differential} {Evolution}: easy {Bayesian} computing for real parameter spaces}, + volume = {16}, + issn = {1573-1375}, + url = {https://doi.org/10.1007/s11222-006-8769-1}, + doi = {10.1007/s11222-006-8769-1}, + number = {3}, + journal = {Statistics and Computing}, + author = {Braak, Cajo J. F. Ter}, + month = sep, + year = {2006}, + pages = {239--249}, +} + +@article{metropolis:1953, + title = {Equation of {State} {Calculations} by {Fast} {Computing} {Machines}}, + volume = {21}, + issn = {0021-9606, 1089-7690}, + url = {https://pubs.aip.org/jcp/article/21/6/1087/202680/Equation-of-State-Calculations-by-Fast-Computing}, + doi = {10.1063/1.1699114}, + number = {6}, + urldate = {2024-11-22}, + journal = {The Journal of Chemical Physics}, + author = {Metropolis, Nicholas and Rosenbluth, Arianna W. and Rosenbluth, Marshall N. and Teller, Augusta H. and Teller, Edward}, + month = jun, + year = {1953}, + pages = {1087--1092}, + file = {Submitted Version:/Users/engs2510/Zotero/storage/LKZVAGE4/Metropolis et al. - 1953 - Equation of State Calculations by Fast Computing M.pdf:application/pdf}, +} + +@software{pybamm-eis, + author = {Dhoot, Rishit and Timms, Robert and Please, Colin}, + title = {{PyBaMM EIS: Efficient Linear Algebra Methods to Determine Li-ion Battery Behaviour}}, + url = {https://www.github.com/pybamm-team/pybamm-eis}, + year = {2024}, + version = {0.1.4}, +} + +@article{Lu:2021, + author = {Lu, Dongliang and {Scott Trimboli}, M. and Fan, Guodong and Zhang, Ruigang and Plett, Gregory L.}, + journal = {Journal of The Electrochemical Society}, + title = {Implementation of a physics-based model for half-cell open-circuit potential and full-cell open-circuit voltage estimates: {P}art {II}. Processing full-cell data}, + year = {2021}, + issn = {0013-4651}, + number = {7}, + pages = {070533}, + volume = {168}, + doi = {10.1149/1945-7111/ac11a5}, + publisher = {IOP Publishing}, +} + +@misc{Loshchilov:2017, + title = {Decoupled {Weight} {Decay} {Regularization}}, + copyright = {arXiv.org perpetual, non-exclusive license}, + url = {https://arxiv.org/abs/1711.05101}, + doi = {10.48550/ARXIV.1711.05101}, + urldate = {2024-12-13}, + publisher = {arXiv}, + author = {Loshchilov, Ilya and Hutter, Frank}, + year = {2017}, + note = {Version Number: 3}, +} + +@inproceedings{Yang:2009, + address = {Coimbatore, India}, + title = {Cuckoo {Search} via Levy flights}, + isbn = {978-1-4244-5053-4}, + url = {http://ieeexplore.ieee.org/document/5393690/}, + doi = {10.1109/NABIC.2009.5393690}, + urldate = {2024-12-13}, + booktitle = {2009 {World} {Congress} on {Nature} \& {Biologically} {Inspired} {Computing} ({NaBIC})}, + publisher = {IEEE}, + author = {Yang, Xin-She and {Suash Deb}}, + year = {2009}, + pages = {210--214}, +} + +@article{Cauchy:1847, + title={M{\'e}thode g{\'e}n{\'e}rale pour la r{\'e}solution des systemes d’{\'e}quations simultan{\'e}es}, + author={Cauchy, Augustin and others}, + journal={Comp. Rend. Sci. Paris}, + volume={25}, + number={1847}, + pages={536--538}, + year={1847} +} + +@article{Aitio:2020, + title = {Bayesian Parameter Estimation Applied to the Li-ion Battery Single Particle Model with Electrolyte Dynamics}, + journal = {IFAC-PapersOnLine}, + volume = {53}, + number = {2}, + pages = {12497-12504}, + year = {2020}, + note = {21st IFAC World Congress}, + issn = {2405-8963}, + doi = {10.1016/j.ifacol.2020.12.1770}, + url = {https://www.sciencedirect.com/science/article/pii/S2405896320323788}, + author = {Antti Aitio and Scott G. Marquis and Pedro Ascencio and David Howey}, +} diff --git a/papers/joss/paper.md b/papers/joss/paper.md new file mode 100644 index 000000000..50480a875 --- /dev/null +++ b/papers/joss/paper.md @@ -0,0 +1,188 @@ +--- +title: 'PyBOP: A Python package for battery model optimisation and parameterisation' +tags: + - python + - battery + - model + - parameter + - inference + - design optimisation +authors: + - name: Brady Planden + orcid: 0000-0002-1082-9125 + affiliation: 1 + equal-contrib: true + - name: Nicola E. Courtier + affiliation: "1, 2" + orcid: 0000-0002-5714-1096 + equal-contrib: true + - name: Martin Robinson + orcid: 0000-0002-1572-6782 + affiliation: 3 + - name: Agriya Khetarpal + orcid: 0000-0002-1112-1786 + affiliation: 4 + - name: Ferran Brosa Planella + affiliation: "2, 5" + orcid: 0000-0001-6363-2812 + - name: David A. Howey + corresponding: true + affiliation: "1, 2" + orcid: 0000-0002-0620-3955 +affiliations: + - name: Department of Engineering Science, University of Oxford, Oxford, United Kingdom + index: 1 + - name: The Faraday Institution, Harwell Campus, Didcot, United Kingdom + index: 2 + - name: Research Software Engineering Group, University of Oxford, Oxford, United Kingdom + index: 3 + - name: Quansight + index: 4 + - name: Mathematics Institute, University of Warwick, Coventry, United Kingdom + index: 5 +date: 19 December 2024 +bibliography: paper.bib +repository: https://github.com/pybop-team/PyBOP +--- + +# Summary + +The Python Battery Optimisation and Parameterisation (`PyBOP`) package provides methods for estimating and optimising battery model parameters, offering both deterministic and stochastic approaches with example workflows. `PyBOP` enables parameter identification from data for various battery models, including electrochemical and equivalent circuit models from the open-source `PyBaMM` package [@Sulzer:2021]. The same approaches enable design optimisation under user-defined operating conditions across various model structures and design goals. `PyBOP` facilitates optimisation and provides diagnostics to examine optimiser performance and convergence of the cost and parameters. Identified parameters can be used for prediction, online estimation and control, and design optimisation, accelerating battery research and development. + +# Statement of need + +`PyBOP` provides a user-friendly, object-oriented interface for optimising battery model parameters. It leverages the open-source `PyBaMM` package [@Sulzer:2021] to formulate and solve battery models. Together, these tools serve a broad audience including students, engineers, and researchers in academia and industry, enabling advanced applications without specialised knowledge of battery modelling, parameter inference, or software development. `PyBOP` emphasises clear diagnostics and workflows to support users with varying domain expertise, and provides access to numerous optimisation and sampling algorithms. These capabilities are enabled through interfaces to `PINTS` [@Clerx:2019], `SciPy` [@SciPy:2020], and `PyBOP`'s implementations of algorithms including adaptive moment estimation with weight decay (AdamW) [@Loshchilov:2017], gradient descent [@Cauchy:1847], and cuckoo search [@Yang:2009]. + +`PyBOP` complements other lithium-ion battery modelling packages built around `PyBaMM`, such as `liionpack` for battery pack simulation [@Tranter2022] and `pybamm-eis` for fast numerical computation of the electrochemical impedance of any battery model, as well as the battery parameter exchange (BPX) standard [@BPX:2023]. Identified `PyBOP` parameters are easily exported to other packages. + +# Architecture + +`PyBOP` is structured around four core components: a `Simulator`, `Cost`, `Problem`, and `Optimiser`/`Sampler`, as shown in \autoref{fig:classes}. The purpose of the `Simulator` is to generate model predictions. For example, `pybop.pybamm.Simulator` interfaces with `PyBaMM` to efficiently construct, discretise and numerically solve a `PyBaMM` model for candidate parameter values. Custom or built-in `Cost` classes evaluate an error measure, likelihood or design metric for the candidate parameter values and simulation result. Multiple costs can be summed with optional weighting. The `Problem` class coordinates simulator and cost evaluation, and the `Optimiser`/`Sampler` classes perform parameter inference through optimisation algorithms or Monte Carlo sampling. This structure ensures extensibility for new optimisation problems with a consistent interface between models and optimisers. + +![The core `PyBOP` architecture with base class interfaces. Each class provides a direct mapping to a step in the optimisation workflow. Note that 'inputs' means candidate parameter values being optimised, following `PyBaMM` conventions. \label{fig:classes}](figures/PyBOP_components.drawio.pdf){ width=90% } + +The `pybamm.Simulator` object returns a solution with corresponding sensitivities, where possible, to enable gradient-based optimisation. Bayesian inference is provided by sampler classes, with Monte Carlo algorithms provided by `PINTS`. In the typical workflow, the classes in \autoref{fig:classes} are constructed in sequence, from left to right. The optimisation result includes a log of the candidate parameters and corresponding cost values. Beyond convergence information, identifiability metrics are provided through Hessian approximation and Sobol sampling from the `SAlib` package. + +Beyond the core architecture, `PyBOP` provides specialised inference and optimisation features. Parameter inference from electrochemical impedance spectroscopy (EIS) simulations is handled through `pybop.pybamm.EISSimulator`, which discretises and linearises the EIS forward model into sparse mass matrix form with an auto-differentiated Jacobian. The result is returned in the frequency domain and is compatible with the same cost classes as in the time-domain simulations. + +The currently available optimisation algorithms are presented in \autoref{tab:optimisers}. Note that `SciPy` minimize includes several gradient-based and gradient-free methods. Hereafter, point-based parameterisation and design-optimisation tasks are referred to as optimisation tasks. This simplification can be justified by comparing \autoref{eqn:parameterisation} and \autoref{eqn:design}; deterministic parameterisation is an optimisation task to minimise distance-based cost between model output and measured values. + +: Currently supported optimisers classified by optimisation type. \label{tab:optimisers} + +| Gradient-based | Evolutionary | (Meta)heuristic | +|:--------------------------------------------------|:--------------------------------------|:---------------------| +| Weight decayed adaptive moment estimation (AdamW) | Covariance matrix adaptation (CMA-ES) | Particle swarm (PSO) | +| Gradient descent | Exponential natural (xNES) | Nelder-Mead | +| `SciPy` minimize | Separable natural (sNES) | Cuckoo search | +| Improved resilient backpropagation (iRProp-/+) | `SciPy` differential evolution | Simulated annealing | + + +Beyond deterministic optimisers (\autoref{tab:optimisers}), `PyBOP` provides Monte Carlo sampling methods to estimate parameter distributions within a Bayesian framework. These methods estimate posterior parameter distributions that can be used to assess uncertainty and practical identifiability. Individual sampler classes are composed from the `PINTS` library, with a base sampler class implemented for interoperability and direct integration with the `Problem` class. Currently supported samplers are listed in \autoref{tab:samplers}. + +: Sampling methods supported by `PyBOP`, classified according to candidate proposal method. \label{tab:samplers} + +| Gradient-based | Adaptive | Slicing | Other | +|:------------------|:---------------------------|:---------------|:-----------------------------| +| Monomial gamma | Delayed rejection adaptive | Rank shrinking | Metropolis adjusted Langevin | +| No-U-turn | Haario Bardenet | Doubling | Emcee hammer | +| Hamiltonian | Haario | Stepout | Metropolis random walk | +| Relativistic | Rao Blackwell | | Differential evolution | + + +# Background + +## Battery models + +In general, battery models (after spatial discretisation) can be written in the form of a differential-algebraic system of equations, +\begin{equation} +\frac{\mathrm{d} \mathbf{x}}{\mathrm{d} t} = f(t,\mathbf{x},\mathbf{\theta}), +\label{dynamics} +\end{equation} +\begin{equation} +0 = g(t, \mathbf{x}, \mathbf{\theta}), +\label{algebraic} +\end{equation} +\begin{equation} +\mathbf{y}(t) = h(t, \mathbf{x}, \mathbf{\theta}), +\label{output} +\end{equation} +with initial conditions +\begin{equation} +\mathbf{x}(0) = \mathbf{x}_0(\mathbf{\theta}). +\label{initial_conditions} +\end{equation} + +Here, $t$ is time, $\mathbf{x}(t)$ are the (spatially discretised) states, $\mathbf{y}(t)$ are the outputs (e.g., the terminal voltage) and $\mathbf{\theta}$ are the parameters. Here the model input(s) (such as current) are implicitly part of the state vector. + +Common battery models include equivalent circuits (e.g., the Thévenin model), the Doyle–Fuller–Newman (DFN) model [@Doyle:1993; @Fuller:1994] based on porous electrode theory, and its reduced-order variants including the single particle model (SPM) [@Planella:2022] and the multi-species multi-reaction (MSMR) model [@Verbrugge:2017]. Simplified models that retain acceptable predictive accuracy at lower computational cost are widely used, for example in battery management systems, while physics-based models are required to understand the impact of physical parameters on performance. However, different model structures will lead to different parameter estimates from the same dataset for parameters-in-common, such as diffusion time or series resistance. + +# Examples + +## Parameterisation + +Battery model parameterisation is challenging due to the large number of parameters compared to the number of possible measurements [@Miguel:2021; @Wang:2022; @Andersson:2022]. A complete parameterisation often requires stepwise identification of parameter subsets from a variety of excitations and datasets [@Chu:2019; @Chen:2020; @Lu:2021; @Kirk:2022]. Parameter identifiability can be poor for some excitations and datasets, requiring improved experimental design and uncertainty-capable identification methods [@Aitio:2020]. + +A generic data-fitting optimisation problem may be formulated as: +\begin{equation} +\min_{\mathbf{\theta}} ~ \mathcal{L}_{(\hat{\mathbf{y}_i})}(\mathbf{\theta}) ~~~ +\textrm{subject to equations (\ref{dynamics})\textrm{-}(\ref{initial_conditions})} +\label{eqn:parameterisation} +\end{equation} + +where $\mathcal{L} : \mathbf{\theta} \mapsto [0,\infty)$ is a cost function that quantifies the agreement between the model output $\mathbf{y}(t)$ and a sequence of observations $(\hat{\mathbf{y}_i})$ measured at times $t_i$. For gradient-based optimisers, the Jacobian of the cost function with respect to unknown parameters, $\partial \mathcal{L} / \partial \theta$, is computed for step-size and directional information. + +We demonstrate the fitting of synthetic data where the model parameters are known, using `PyBaMM`'s SPM with contact resistance. We target two parameters: the lithium diffusivity in the negative electrode active material particles ("negative particle diffusivity") and the contact resistance, with true values [3.3e-14 $\text{m}^2/\text{s}$, 10 m$\text{\Omega}$]. We generate time-domain data for a one-hour discharge from 100% to 0% state of charge (1C rate) followed by 30 minutes relaxation. The output voltage is corrupted with zero-mean Gaussian noise of amplitude 2 mV (blue dots in \autoref{fig:inference-time-landscape} (left)). Initial states are assumed known, although this is not generally necessary. The `PyBOP` repository contains [example notebooks](https://github.com/pybop-team/PyBOP/tree/develop/examples/notebooks) illustrating similar inference processes. The underlying cost landscape to be explored by the optimiser is shown in \autoref{fig:inference-time-landscape} (right), with the initial position and true values marked. In general, the true values are unknown. + +![A synthetic dataset (left) and cost landscape (right) depicting a time-series parameterisation problem using the root-mean-squared error cost function. \label{fig:inference-time-landscape}](figures/combined/sim-landscape.pdf){ width=100% } + +`PyBOP` can generate and fit EIS data using methods from `pybamm-eis` [@pybamm-eis]. Using `PyBaMM`'s SPM with double-layer capacitance and contact resistance, \autoref{fig:impedance-landscape} shows numerical EIS predictions alongside the cost landscape for the corresponding inference task. At the time of publication, gradient-based optimisation and sampling methods are not available for EIS simulators. + +![Data and model fit (left) and cost landscape (right) for a frequency-domain EIS parameterisation, at 5% SOC, using the root-mean-squared error cost function. \label{fig:impedance-landscape}](figures/combined/impedance.pdf){ width=100% } + +We continue here with time-domain identification (\autoref{fig:inference-time-landscape}), however time- and frequency-domain problems may be combined for improved parameterisation. As gradient information is available for our time-domain example, the choice of distance-based cost function and optimiser is unconstrained. Due to the difference in magnitude between the two parameters, we apply a logarithmic transformation that transforms the search space to allow for a common step size, improving convergence. As a demonstration of `PyBOP`'s parameterisation capabilities, \autoref{fig:convergence-min-max} (left) shows convergence rates for distance-minimising cost functions, while \autoref{fig:convergence-min-max} (right) shows analogous results for likelihood maximisation. Optimisation is performed using `SciPy` minimize with the gradient-based BFGS method. + +![Convergence of the BFGS method for various cost (left) and likelihood (right) functions. \label{fig:convergence-min-max}](figures/combined/converge.pdf){ width=100% } + +Using the same model and parameters, we compare example convergence rates of various algorithms across several categories: gradient-based methods in \autoref{fig:optimiser-inference1} (left), evolutionary strategies in \autoref{fig:optimiser-inference1} (middle) and (meta)heuristics in \autoref{fig:optimiser-inference1} (right) using a mean-squared-error cost. \autoref{fig:optimiser-inference2} shows the optimiser's exploration of the cost landscape, with the three rows showing the gradient-based optimisers (top), evolution strategies (middle), and (meta)heuristics (bottom). Optimiser performance depends on the cost landscape, initial guess or prior for each parameter, and the hyperparameters for each problem. + +![Convergence in the parameter values for optimisation algorithms available in `PyBOP`. \label{fig:optimiser-inference1}](figures/combined/optimisers_parameters.pdf){ width=100% } + +![Cost landscape plots showing the optimisation traces of 12 different optimisers. \label{fig:optimiser-inference2}](figures/combined/contour_subplot.pdf){ width=100% } + +This example parameterisation task can also be approached from a Bayesian perspective, using `PyBOP`'s sampler methods. First, we introduce Bayes' rule, + +\begin{equation} +P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}, +\label{eqn:bayes_theorem} +\end{equation} + +where $P(\theta|D)$ is the posterior parameter distribution, $P(D|\theta)$ is the likelihood function, $P(\theta)$ is the prior parameter distribution, and $P(D)$ is the model evidence, or marginal likelihood, which acts as a normalising constant. For maximum likelihood estimation or maximum a posteriori estimation, one wishes to maximise $P(D|\theta)$ or $P(\theta|D)$, respectively, formulated as an optimisation problem as per \autoref{eqn:parameterisation}. + +To estimate the full posterior parameter distribution, however, one must use sampling or other inference methods to reconstruct $P(\theta|D)$. The posterior distribution provides information about the uncertainty of the identified parameters, e.g., by calculating the variance or other moments. Monte Carlo methods available from the probabilistic inference on noisy time-series (`PINTS`) package include gradient-based methods such as no-u-turn [@NUTS:2011] and Hamiltonian [@Hamiltonian:2011], heuristic methods such as differential evolution [@DiffEvolution:2006], and conventional methods based on random sampling with rejection criteria [@metropolis:1953]. \autoref{fig:posteriors} shows sampled posteriors for the synthetic model using an adaptive covariance-based sampler called Haario Bardenet [@Haario:2001]. + +![Posterior distributions of model parameters and observation noise; shaded areas show the 95th percentile credible interval. \label{fig:posteriors}](figures/combined/posteriors.pdf){ width=100% } + +## Design optimisation + +`PyBOP` supports design optimisation to guide device design development by identifying parameter sensitivities that can unlock improvements in performance. Design workflows are similar to parameterisation workflows, but the aim is to maximise a design metric rather than minimise a distance-based cost function. `PyBOP` performs maximisation by minimising the negative cost. An example design metric is the gravimetric energy (or power) density given by the integral of the discharge energy (or power) normalised by the cell mass. Such metrics are typically quantified for operating conditions such as a 1C discharge, at a given temperature. + +In general, design optimisation can be written as a constrained optimisation problem, +\begin{equation} +\min_{\mathbf{\theta} \in \Omega} ~ -\mathcal{L}(\mathbf{\theta}) ~~~ +\textrm{subject to equations (\ref{dynamics})\textrm{-}(\ref{initial_conditions}),} +\label{eqn:design} +\end{equation} + +where $\mathcal{L} : \mathbf{\theta} \mapsto [0,\infty)$ is a cost function that quantifies the desirability of the design and $\Omega$ is the set of allowable parameter values. + +We consider maximising gravimetric energy density subject to constraints on two of the geometric electrode parameters [@Couto:2023]. We use the `PyBaMM` single particle model with electrolyte (SPMe) to investigate the impact of positive electrode thickness and active material volume fraction on energy density. Since the total volume fraction must sum to unity, the positive electrode porosity is defined relative to the active material volume fraction. The 1C rate can also be optimised (via the nominal capacity parameter) or defined as a function of the parameters for each design. + +![Initial and optimised voltage profiles alongside the gravimetric energy density cost landscape. \label{fig:design_gravimetric}](figures/combined/design.pdf){ width=100% } + +\autoref{fig:design_gravimetric} (left) shows the predicted improvement in the discharge profile between the initial and optimised parameter values for a fixed-rate 1C discharge selected from the initial design and (right) the Nelder-Mead search over the parameter space. + +# Acknowledgements + +We gratefully acknowledge all [contributors](https://github.com/pybop-team/PyBOP/graphs/contributors) to `PyBOP`. This work was supported by the Faraday Institution Multiscale Modelling project (FIRG059), UKRI's Horizon Europe Guarantee (10038031), and EU IntelLiGent project (101069765). + +# References diff --git a/papers/joss/param_plots.py b/papers/joss/param_plots.py new file mode 100644 index 000000000..585fc2087 --- /dev/null +++ b/papers/joss/param_plots.py @@ -0,0 +1,1028 @@ +# A script to generate parameterisation plots for the JOSS paper. + + +import matplotlib.pyplot as plt +import numpy as np +import plotly +import pybamm +from matplotlib.ticker import ScalarFormatter + +import pybop +from pybop.plot import PlotlyManager + +go = PlotlyManager().go +px = PlotlyManager().px +make_subplots = PlotlyManager().make_subplots +plt.rcParams.update({"text.usetex": True}) # Enable LaTeX +np.random.seed(8) # Set random seed for reproducibility +axis_font_size = 24 +tick_font_size = 16 + +# Choose which plots to show and save +create_plot = {} +create_plot["simulation"] = True +create_plot["landscape"] = True +create_plot["minimising"] = True +create_plot["maximising"] = True +create_plot["gradient"] = True +create_plot["evolution"] = True +create_plot["heuristic"] = True +create_plot["posteriors"] = True +create_plot["eis"] = True + +# Define model, solver and parameter values +model = pybamm.lithium_ion.SPM(options={"contact resistance": "true"}) +solver = pybamm.IDAKLUSolver(rtol=1e-7, atol=1e-7) +parameter_values = pybamm.ParameterValues("Chen2020") +parameter_values["Contact resistance [Ohm]"] = 0.01 + +# Define experiment +experiment = pybamm.Experiment( + [ + "Discharge at 1C until 2.8 V (20 second period)", + "Rest for 30 minutes (20 second period)", + ] +) + +# Generate a synthetic dataset with Gaussian noise +sigma = 0.005 +solution = pybamm.Simulation( + model, parameter_values=parameter_values, experiment=experiment, solver=solver +).solve() +values = solution["Voltage [V]"].data +corrupt_values = values + np.random.normal(0, sigma, len(values)) + +if create_plot["simulation"]: + # Plot the data and the simulation + simulation_plot_dict = pybop.plot.StandardPlot( + x=solution["Time [s]"].data, + y=[corrupt_values, solution["Battery open-circuit voltage [V]"].data, values], + trace_names=[ + "Voltage w. noise", + "Open-circuit voltage", + "Voltage", + ], + ) + simulation_plot_dict.traces[0].mode = "markers" + simulation_fig = simulation_plot_dict(show=False) + simulation_fig.update_layout( + width=600, + height=600, + xaxis=dict( + title=dict(text="Time / s", font_size=axis_font_size), + showline=True, + linewidth=1, + linecolor="black", + mirror=True, + tickfont_size=tick_font_size, + ), + yaxis=dict( + title=dict(text="Voltage / V", font_size=axis_font_size), + showline=True, + linewidth=1, + linecolor="black", + mirror=True, + tickfont_size=tick_font_size, + ), + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font_size=tick_font_size, + ), + margin=dict(t=60, b=84, r=50, l=15), + ) + simulation_fig.show() + simulation_fig.write_image("figures/individual/simulation.pdf") + +# Form dataset +dataset = pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": corrupt_values, + } +) + +# Fitting parameters +true_value = [ + parameter_values[p] + for p in ["Contact resistance [Ohm]", "Negative particle diffusivity [m2.s-1]"] +] +initial_value = [0.02, 9e-14] +parameter_values.update( + { + "Contact resistance [Ohm]": pybop.Parameter( + initial_value=initial_value[0], + prior=pybop.Gaussian(0.02, 0.005), + transformation=pybop.ScaledTransformation(coefficient=200), + bounds=[0.005, 0.025], + ), + "Negative particle diffusivity [m2.s-1]": pybop.Parameter( + initial_value=initial_value[1], + prior=pybop.Gaussian(9e-14, 2e-14), + transformation=pybop.LogTransformation(), + bounds=[1.9e-14, 12e-14], + ), + } +) + +# Build the problem +simulator = pybop.pybamm.Simulator( + model, parameter_values=parameter_values, protocol=dataset, solver=solver +) +cost = pybop.SumSquaredError(dataset) +problem = pybop.Problem(simulator, cost) + +if create_plot["landscape"]: + # Plot the cost landscape with the initial and true values + landscape_fig = pybop.plot.contour( + problem, + steps=25, + title=None, + show=False, + xaxis=dict( + title=dict(text="Contact resistance / Ω", font_size=axis_font_size), + tickfont_size=tick_font_size, + ), + yaxis=dict( + title=dict( + text="Negative particle diffusivity / m2 s-1", + font_size=axis_font_size, + standoff=0, + ), + tickfont_size=tick_font_size, + exponentformat="power", + ), + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font_size=tick_font_size, + ), + coloraxis_colorbar=dict(tickfont_size=tick_font_size), + margin=dict(t=50), + ) + landscape_fig.update_traces(colorbar=dict(tickfont_size=tick_font_size)) + landscape_fig.add_trace( + go.Scatter( + x=[initial_value[0]], + y=[initial_value[1]], + mode="markers", + marker=dict( + symbol="x", + color="white", + line_color="black", + line_width=1, + size=16, + showscale=False, + ), + name="Initial values", + ) + ) + landscape_fig.add_trace( + go.Scatter( + x=[true_value[0]], + y=[true_value[1]], + mode="markers", + marker=dict( + symbol="cross", + color="black", + line_color="white", + line_width=1, + size=16, + showscale=False, + ), + name="True values", + ) + ) + landscape_fig.show() + landscape_fig.write_image("figures/individual/landscape.pdf") + +# Categorise the costs +minimising_cost_classes = [ + pybop.Minkowski, # largest + pybop.SumSquaredError, + pybop.SumOfPower, + pybop.RootMeanSquaredError, # smallest +] +maximising_cost_classes = [ + pybop.GaussianLogLikelihood, + pybop.GaussianLogLikelihoodKnownSigma, + pybop.LogPosterior, + pybop.LogPosterior, +] + + +if create_plot["minimising"]: + # Show cost convergence using same optimiser for different cost functions + convergence_traces = [] + for cost in minimising_cost_classes: + # Define keyword arguments for the cost class + kwargs = {} + if cost is pybop.SumOfPower: + kwargs["p"] = 2.5 + + # Define the problem and optimiser + cost = cost(dataset, **kwargs) + problem = pybop.Problem(simulator, cost) + options = pybop.SciPyMinimizeOptions(maxiter=50, method="BFGS", jac=True) + optim = pybop.SciPyMinimize(problem, options=options) + + # Run optimisation + result = optim.run() + print(result) + print("True parameter values:", true_value) + + # Plot convergence + cost_log = result.cost + iteration_numbers = list(range(1, len(cost_log) + 1)) + convergence_plot_dict = pybop.plot.StandardPlot( + x=iteration_numbers, + y=cost_log, + trace_names=[cost.name], + trace_options={"line": {"width": 4, "dash": "dash"}}, + ) + convergence_traces.extend(convergence_plot_dict.traces) + + # Plot minimising convergence traces together + convergence_fig = go.Figure( + data=convergence_traces, + layout=dict( + xaxis=dict( + title=dict(text="Evaluation", font_size=axis_font_size), + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis=dict( + title=dict(text="Cost", font_size=axis_font_size), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + legend=dict( + yanchor="top", + y=0.95, + xanchor="right", + x=0.99, + font_size=tick_font_size, + bordercolor="black", + borderwidth=1, + ), + plot_bgcolor="white", + width=600, + height=600, + ), + ) + convergence_fig.update_traces(dict(mode="markers")) + convergence_fig.show() + convergence_fig.write_image("figures/individual/convergence_minimising.pdf") + + +if create_plot["maximising"]: + ## Do the same for the maximising cost functions + convergence_traces = [] + first_MAP = True + + for cost in maximising_cost_classes: + if cost is pybop.GaussianLogLikelihoodKnownSigma: + cost = cost(dataset, sigma0=sigma) + elif cost is pybop.GaussianLogLikelihood: + cost = cost(dataset, sigma0=4 * sigma) + elif cost is pybop.LogPosterior and first_MAP: + cost = cost( + log_likelihood=pybop.GaussianLogLikelihoodKnownSigma( + dataset, sigma0=sigma + ) + ) + first_MAP = False + elif cost is pybop.LogPosterior: + cost = cost(log_likelihood=pybop.GaussianLogLikelihood(dataset)) + + # Define the problem and optimiser + problem = pybop.Problem(simulator, cost) + options = pybop.SciPyMinimizeOptions(maxiter=50, method="BFGS", jac=True) + optim = pybop.SciPyMinimize(problem, options=options) + + # Run optimisation + result = optim.run() + print(result) + print("True parameter values:", true_value) + + # Plot convergence + cost_log = result.cost + iteration_numbers = list(range(1, len(cost_log) + 1)) + convergence_plot_dict = pybop.plot.StandardPlot( + x=iteration_numbers, + y=cost_log, + trace_names=cost.name + + " " + + ( + cost.log_likelihood.name if isinstance(cost, pybop.LogPosterior) else "" + ), + trace_options={"line": {"width": 4, "dash": "dash"}}, + ) + convergence_traces.extend(convergence_plot_dict.traces) + + # Plot maximising convergence traces together + convergence_fig = go.Figure( + data=convergence_traces, + layout=dict( + xaxis=dict( + title=dict(text="Evaluation", font_size=axis_font_size), + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis=dict( + title=dict(text="Likelihood", font_size=axis_font_size), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + legend=dict( + yanchor="bottom", + y=0.15, + xanchor="right", + x=0.99, + font_size=tick_font_size, + bordercolor="black", + borderwidth=1, + ), + plot_bgcolor="white", + width=600, + height=600, + ), + ) + convergence_fig.update_traces(dict(mode="markers")) + convergence_fig.show() + convergence_fig.write_image("figures/individual/convergence_maximising.pdf") + + +# Categorise the optimisers +gradient_optimiser_classes = [ + pybop.AdamW, + pybop.IRPropMin, + pybop.GradientDescent, + pybop.SciPyMinimize, # with BFGS and jac=True +] +evolution_optimiser_classes = [ + pybop.CMAES, + pybop.XNES, + pybop.SNES, + pybop.SciPyDifferentialEvolution, +] +heuristic_optimiser_classes = [ + pybop.PSO, + pybop.NelderMead, + pybop.CuckooSearch, + pybop.SimulatedAnnealing, +] + +# Define the problem +cost = pybop.RootMeanSquaredError(dataset) +problem = pybop.Problem(simulator, cost) + +# Create subplot figure +max_optims = max( + len(gradient_optimiser_classes), + len(evolution_optimiser_classes), + len(heuristic_optimiser_classes), +) +subplot_contour_fig = make_subplots( + rows=3, + cols=max_optims, + shared_yaxes=True, + shared_xaxes=True, + x_title="Contact resistance / Ω", + y_title="Negative particle diffusivity / m2 s-1", + horizontal_spacing=0.03, + vertical_spacing=0.04, + subplot_titles=[ + cls.__name__ + for cls in [ + *gradient_optimiser_classes, + *evolution_optimiser_classes, + *heuristic_optimiser_classes, + ] + ], +) +bounds = problem.parameters.get_bounds_for_plotly() + +if create_plot["gradient"]: + # Create subplot structure + num_optimisers = len(gradient_optimiser_classes) + parameter_traces = [] + + for i, optimiser in enumerate(gradient_optimiser_classes): + # Set optimiser options + if optimiser is pybop.SciPyMinimize: + options = pybop.SciPyMinimizeOptions(method="BFGS", jac=True, maxiter=50) + else: + options = pybop.PintsOptions( + max_unchanged_iterations=50, max_evaluations=250 + ) + + # Construct the optimiser + optim = optimiser(problem, options=options) + + if optimiser is pybop.AdamW: + optim.optimiser.b1 = 0.85 + optim.optimiser.b2 = 0.9 + optim.optimiser.lam = 0.005 + if optimiser is pybop.GradientDescent: + optim.optimiser.set_learning_rate(eta=[11, 4.5]) + + # Run optimisation + print(optim.name) + result = optim.run() + print(result) + print("True parameter values:", true_value) + + # Plot the parameter traces + parameter_fig = pybop.plot.parameters(result, show=False, title=None) + parameter_fig.update_layout( + legend=dict( + orientation="v", + yanchor="bottom", + y=1.02, + xanchor="left", + x=-0.05, + font_size=tick_font_size, + ) + ) + parameter_fig.update_traces(dict(mode="markers")) + colours = parameter_fig.layout.template.layout.colorway + for j in range(len(parameter_fig.data)): + parameter_fig.data[j].name = optim.name + parameter_fig.data[j].marker.color = colours[i] + if j > 0: + parameter_fig.data[j].showlegend = False + parameter_traces.extend(parameter_fig.data) + + # Collect contour plot data + contour = pybop.plot.contour( + result, + steps=25, + title="", + showlegend=False, + show=False, + margin=dict(l=20, r=20, t=20, b=20), + ) + if i == num_optimisers - 1: + contour.update_traces( + showscale=(i == num_optimisers - 1), + colorbar=dict( + title=dict(text="Cost", font_size=axis_font_size), + tickfont_size=tick_font_size, + ), + selector=dict(type="contour"), + ) + else: + contour.update_traces(showscale=False, selector=dict(type="contour")) + + # Add all traces from the contour figure to the subplot + for trace in contour.data: + subplot_contour_fig.add_trace(trace, row=1, col=i + 1) + + # Plot the parameter traces together + parameter_fig.update_layout( + width=480, + height=910, + plot_bgcolor="white", + xaxis=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis=dict( + title=dict( + text="Contact resistance / Ω", font_size=axis_font_size, standoff=30 + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + range=bounds[0], + ), + legend=dict( + yanchor="bottom", y=1.02, xanchor="left", x=-0.05, font_size=tick_font_size + ), + xaxis2=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis2=dict( + title=dict( + text="Negative particle diffusivity / m2 s-1", + font_size=axis_font_size, + standoff=0, + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + exponentformat="power", + linewidth=1, + linecolor="black", + range=bounds[1], + ), + ) + parameter_fig.data = [] + parameter_fig.add_traces(parameter_traces) + parameter_fig = plotly.subplots.make_subplots(figure=parameter_fig, rows=2, cols=1) + parameter_fig.show() + parameter_fig.write_image("figures/individual/gradient_parameters.pdf") + +if create_plot["evolution"]: + # Create subplot structure + num_optimisers = len(evolution_optimiser_classes) + + parameter_traces = [] + for i, optimiser in enumerate(evolution_optimiser_classes): + # Set optimiser options + if optimiser is pybop.SciPyDifferentialEvolution: + options = pybop.SciPyDifferentialEvolutionOptions(maxiter=50, popsize=3) + else: + options = pybop.PintsOptions( + max_unchanged_iterations=50, max_evaluations=250 + ) + + # Run optimisation + optim = optimiser(problem, options=options) + print(optim.name) + result = optim.run() + print(result) + print("True parameter values:", true_value) + + # Plot the parameter traces + parameter_fig = pybop.plot.parameters(result, show=False, title=None) + parameter_fig.update_traces(dict(mode="markers")) + colours = parameter_fig.layout.template.layout.colorway + for j in range(len(parameter_fig.data)): + parameter_fig.data[j].name = optim.name + parameter_fig.data[j].marker.color = colours[i] + if j > 0: + parameter_fig.data[j].showlegend = False + parameter_traces.extend(parameter_fig.data) + + # Collect contour plot data + contour = pybop.plot.contour( + result, + steps=25, + title="", + showlegend=False, + show=False, + margin=dict(l=20, r=20, t=20, b=20), + ) + contour.update_traces(showscale=False, selector=dict(type="contour")) + # Add all traces from the contour figure to the subplot + for trace in contour.data: + subplot_contour_fig.add_trace(trace, row=2, col=i + 1) + + # Plot the parameter traces together + parameter_fig.update_layout( + width=480, + height=910, + plot_bgcolor="white", + xaxis=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis=dict( + title=dict( + text="Contact resistance / Ω", font_size=axis_font_size, standoff=30 + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + range=bounds[0], + ), + legend=dict( + yanchor="bottom", y=1.02, xanchor="left", x=-0.05, font_size=tick_font_size + ), + xaxis2=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis2=dict( + title=dict( + text="Negative particle diffusivity / m2 s-1", + font_size=axis_font_size, + standoff=0, + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + exponentformat="power", + linewidth=1, + linecolor="black", + range=bounds[1], + ), + ) + parameter_fig.data = [] + parameter_fig.add_traces(parameter_traces) + parameter_fig = plotly.subplots.make_subplots(figure=parameter_fig, rows=2, cols=1) + parameter_fig.show() + parameter_fig.write_image("figures/individual/evolution_parameters.pdf") + + +if create_plot["heuristic"]: + # Create subplot structure + num_optimisers = len(evolution_optimiser_classes) + + parameter_traces = [] + for i, optimiser in enumerate(heuristic_optimiser_classes): + # Set optimiser options + options = pybop.PintsOptions(max_unchanged_iterations=50, max_evaluations=250) + + # Construct the optimiser + optim = optimiser(problem, options=options) + + if isinstance(optim, pybop.SimulatedAnnealing): + optim.optimiser.temperature = 0.1 + + # Run optimisation + print(optim.name) + result = optim.run() + print(result) + print("True parameter values:", true_value) + + # Plot the parameter traces + parameter_fig = pybop.plot.parameters(result, show=False, title=None) + parameter_fig.update_traces(dict(mode="markers")) + colours = parameter_fig.layout.template.layout.colorway + for j in range(len(parameter_fig.data)): + parameter_fig.data[j].name = optim.name + parameter_fig.data[j].marker.color = colours[i] + if j > 0: + parameter_fig.data[j].showlegend = False + parameter_traces.extend(parameter_fig.data) + + # Collect contour plot data + contour = pybop.plot.contour( + result, + steps=25, + title="", + showlegend=False, + show=False, + margin=dict(l=20, r=20, t=20, b=20), + ) + contour.update_traces(showscale=False, selector=dict(type="contour")) + + # Add all traces from the contour figure to the subplot + for trace in contour.data: + subplot_contour_fig.add_trace(trace, row=3, col=i + 1) + + # Update layout to configure the color bar and plot dimensions + subplot_contour_fig.update_xaxes( + dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + range=bounds[0], + ) + ) + subplot_contour_fig.update_yaxes( + dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + range=bounds[1], + ) + ) + subplot_contour_fig.update_layout( + showlegend=False, + height=400 * 3, + width=400 * max_optims, + ) + subplot_contour_fig.update_annotations(font_size=axis_font_size) + subplot_contour_fig.update_annotations( + y=-0.01, font_size=axis_font_size, selector={"text": "Contact resistance / Ω"} + ) + + # Show figure and save image + subplot_contour_fig.show() + subplot_contour_fig.write_image("figures/combined/contour_subplot.pdf") + + # Plot the parameter traces together + parameter_fig.update_layout( + width=480, + height=910, + plot_bgcolor="white", + xaxis=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis=dict( + title=dict( + text="Contact resistance / Ω", font_size=axis_font_size, standoff=30 + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + range=bounds[0], + ), + legend=dict( + yanchor="bottom", y=1.02, xanchor="left", x=-0.05, font_size=tick_font_size + ), + xaxis2=dict( + title_font_size=axis_font_size, + tickfont_size=tick_font_size, + linewidth=1, + linecolor="black", + ), + yaxis2=dict( + title=dict( + text="Negative particle diffusivity / m2 s-1", + font_size=axis_font_size, + standoff=0, + ), + gridcolor=px.colors.qualitative.Pastel2[7], + zerolinecolor=px.colors.qualitative.Pastel2[7], + zerolinewidth=1, + tickfont_size=tick_font_size, + exponentformat="power", + linewidth=1, + linecolor="black", + range=bounds[1], + ), + ) + parameter_fig.data = [] + parameter_fig.add_traces(parameter_traces) + parameter_fig = plotly.subplots.make_subplots(figure=parameter_fig, rows=2, cols=1) + parameter_fig.show() + parameter_fig.write_image("figures/individual/heuristic_parameters.pdf") + + +if create_plot["posteriors"]: + sigma0 = pybop.Parameter( + initial_value=sigma, + prior=pybop.Uniform(1e-8 * sigma, 10 * sigma), + bounds=[1e-8, 10 * sigma], + ) + likelihood = pybop.GaussianLogLikelihood(dataset, sigma0=sigma0) + posterior = pybop.Problem(simulator, pybop.LogPosterior(likelihood)) + + options = pybop.PintsSamplerOptions( + n_chains=5, + max_iterations=3500, + warm_up_iterations=1500, + cov=posterior.parameters.get_sigma0(transformed=True), + ) + sampler = pybop.HaarioBardenetACMC(posterior, options=options) + + result = sampler.run() + print(result) + print("True parameter values:", [true_value, sigma]) + summary = pybop.PosteriorSummary(result.chains) + print(summary.rhat()) + print(summary.effective_sample_size(mixed_chains=True)) + + # Create a grid for subplots + fig = plt.figure(figsize=(15, 6)) + gs = fig.add_gridspec(1, 3, height_ratios=[1]) + + # Function to format axis + def format_axis(ax): + formatter = ScalarFormatter(useMathText=True) + formatter.set_powerlimits((-1, 1)) # Set limits to use scientific notation + formatter.set_scientific(True) # Ensure scientific notation + ax.xaxis.set_major_formatter(formatter) + ax.yaxis.set_major_formatter(formatter) + ax.tick_params(axis="both", which="major", labelsize=tick_font_size) + + # Subplot for parameter 0 + ax1 = fig.add_subplot(gs[0, 0]) + ax1.hist( + summary.all_samples[:, 0], + bins=50, + density=False, + alpha=0.6, + color="tab:blue", + ) + ax1.set_xlabel(r"Contact resistance / $\Omega$", fontsize=18) + ax1.set_ylabel("Density", fontsize=18) + ax1.set_ylim(0, 775) + ax1.tick_params(axis="both", which="major", labelsize=tick_font_size) + ax1.axvspan( + summary.get_summary_statistics()[("ci_lower")][0], + summary.get_summary_statistics()[("ci_upper")][0], + alpha=0.1, + color="tab:blue", + ) + ax1.axvline(x=true_value[0], color="black") + format_axis(ax1) + + # Subplot for parameter 1 + ax2 = fig.add_subplot(gs[0, 1]) + ax2.hist( + summary.all_samples[:, 1], + bins=50, + density=False, + alpha=0.6, + color="tab:red", + label="$R_n$", + ) + ax2.set_xlabel(r"Negative particle diffusivity / m$^2\,$s$^{-1}$", fontsize=18) + ax2.set_ylim(0, 775) + ax2.tick_params(axis="both", which="major", labelsize=tick_font_size) + ax2.axvspan( + summary.get_summary_statistics()[("ci_lower")][1], + summary.get_summary_statistics()[("ci_upper")][1], + alpha=0.1, + color="tab:red", + ) + ax2.axvline(x=true_value[1], color="black") + format_axis(ax2) + + # Subplot for sigma + ax3 = fig.add_subplot(gs[0, 2]) + ax3.hist( + summary.all_samples[:, 2], + bins=50, + density=False, + alpha=0.6, + color="tab:purple", + ) + ax3.set_xlabel("Observation Noise / V", fontsize=18) + ax3.set_ylim(0, 775) + ax3.tick_params(axis="both", which="major", labelsize=tick_font_size) + ax3.axvspan( + summary.get_summary_statistics()[("ci_lower")][2], + summary.get_summary_statistics()[("ci_upper")][2], + alpha=0.1, + color="tab:purple", + ) + ax3.axvline(x=sigma, color="black") + format_axis(ax3) + + # Adjust layout + plt.tight_layout() + plt.savefig("figures/combined/posteriors.pdf") + plt.show() + + +if create_plot["eis"]: + + def noise(sigma, values): + # Generate real part noise + real_noise = np.random.normal(0, sigma, values) + + # Generate imaginary part noise + imag_noise = np.random.normal(0, sigma, values) + + # Combine them into a complex noise + return real_noise + 1j * imag_noise + + var = pybamm.standard_spatial_vars + var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 30, var.r_p: 30} + + # Construct model + model = pybamm.lithium_ion.SPM( + options={"surface form": "differential", "contact resistance": "true"} + ) + + initial_state = {"Initial SoC": 0.05} + n_frequency = 35 + sigma0 = 5e-4 + f_eval = np.logspace(-3.5, 4, n_frequency) + + # Create synthetic data for parameter inference + simulator = pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + f_eval=f_eval, + initial_state=initial_state, + var_pts=var_pts, + ) + solution = simulator.solve( + inputs={ + "Contact resistance [Ohm]": true_value[0], + "Negative particle diffusivity [m2.s-1]": true_value[1], + } + ) + + # Form dataset + dataset = pybop.Dataset( + { + "Frequency [Hz]": f_eval, + "Current function [A]": np.zeros(len(f_eval)), + "Impedance": solution["Impedance"].data + + noise(sigma0, len(solution["Impedance"].data)), + }, + domain="Frequency [Hz]", + ) + + # Build the problem + cost = pybop.RootMeanSquaredError(dataset, target=["Impedance"]) + problem = pybop.Problem(simulator, cost) + + # Set up and run the optimiser + options = pybop.PintsOptions(max_iterations=75, min_iterations=75) + optim = pybop.CMAES(problem, options=options) + result = optim.run() + + parameter_fig = pybop.plot.nyquist( + problem, + result.best_inputs, + title="", + width=600, + height=600, + margin=dict(t=60, b=84, r=50, l=15), + xaxis=dict(title_font_size=axis_font_size, linewidth=1), + yaxis=dict(title_font_size=axis_font_size, linewidth=1), + ) + parameter_fig[0].data[1].update(line=dict(color="#00CC97")) + parameter_fig[0].write_image("figures/individual/impedance_spectrum.pdf") + + landscape_fig = pybop.plot.contour( + problem, + steps=25, + show=False, + xaxis=dict( + title=dict(text="Contact resistance / Ω", font_size=axis_font_size), + tickfont_size=tick_font_size, + ), + yaxis=dict( + title=dict( + text="Negative particle diffusivity / m2 s-1", + font_size=axis_font_size, + standoff=0, + ), + tickfont_size=tick_font_size, + exponentformat="power", + ), + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font_size=tick_font_size, + ), + coloraxis_colorbar=dict(tickfont_size=tick_font_size), + margin=dict(t=50), + title=None, + ) + landscape_fig.add_trace( + go.Scatter( + x=[initial_value[0]], + y=[initial_value[1]], + mode="markers", + marker=dict( + symbol="x", + color="white", + line_color="black", + line_width=1, + size=16, + showscale=False, + ), + name="Initial values", + ) + ) + landscape_fig.add_trace( + go.Scatter( + x=[true_value[0]], + y=[true_value[1]], + mode="markers", + marker=dict( + symbol="cross", + color="black", + line_color="white", + line_width=1, + size=16, + showscale=False, + ), + name="True values", + ) + ) + landscape_fig.show() + landscape_fig.write_image("figures/individual/impedance_contour.pdf") diff --git a/pybop/plot/nyquist.py b/pybop/plot/nyquist.py index 02fc6a94b..80f7eb77a 100644 --- a/pybop/plot/nyquist.py +++ b/pybop/plot/nyquist.py @@ -78,11 +78,7 @@ def nyquist(problem, inputs: Inputs = None, show=True, **layout_kwargs): scaleratio=1, ), legend=dict( - x=0.02, - y=0.98, - bgcolor="rgba(255, 255, 255, 0.5)", - bordercolor="black", - borderwidth=1, + orientation="h", yanchor="bottom", y=1.02, xanchor="right", x=1 ), width=600, height=600, @@ -97,8 +93,8 @@ def nyquist(problem, inputs: Inputs = None, show=True, **layout_kwargs): plot_dict.traces[0].update( mode="lines+markers", - line=dict(color="blue", width=2), - marker=dict(size=8, color="blue", symbol="circle"), + line=dict(color="#00CC96", width=2), + marker=dict(size=8, color="#00CC96", symbol="circle"), ) target_trace = plot_dict.create_trace( @@ -106,27 +102,13 @@ def nyquist(problem, inputs: Inputs = None, show=True, **layout_kwargs): y=-target_output[var].imag, name="Reference", mode="markers", - marker=dict(size=8, color="red", symbol="circle-open"), + marker=dict(size=8, color="#636EFA", symbol="circle-open"), showlegend=True, ) plot_dict.traces.append(target_trace) fig = plot_dict(show=False) - # Add minor gridlines - fig.update_xaxes( - showgrid=True, - gridwidth=1, - gridcolor="lightgray", - minor=dict(showgrid=True, gridwidth=0.5, gridcolor="lightgray"), - ) - fig.update_yaxes( - showgrid=True, - gridwidth=1, - gridcolor="lightgray", - minor=dict(showgrid=True, gridwidth=0.5, gridcolor="lightgray"), - ) - # Overwrite with user-kwargs fig.update_layout(**layout_kwargs) if show: