diff --git a/examples/data/Baumhofer2014/README.md b/examples/data/Baumhofer2014/README.md new file mode 100644 index 000000000..8774e19b5 --- /dev/null +++ b/examples/data/Baumhofer2014/README.md @@ -0,0 +1,14 @@ +# Degradation data + +This repository contains a subset of the data published in the research paper: +**P. Attia et al., "Review—“Knees” in Lithium-Ion Battery Aging Trajectories,"** *Journal of The Electrochemical Society*, vol. 169, p. 060517, June 2022. [DOI: 10.1149/1945-7111/ac6d13](https://doi.org/10.1149/1945-7111/ac6d13), +which digitised it from the source paper: +**T. Baumhöfer et al., "Production caused variation in capacity aging trend and correlation to initial cell performance,"** *Journal of Power Sources*, vol. 247, p. 332-338, 2014. [DOI: 10.1016/j.jpowsour.2013.08.108](https://doi.org/10.1016/j.jpowsour.2013.08.108) + +## Citation + +If you use this dataset in your work, please cite the original publications as follows: + +```plaintext +P. Attia et al., "Review—“Knees” in Lithium-Ion Battery Aging Trajectories," Journal of The Electrochemical Society, vol. 169, p. 060517, June 2022, doi: 10.1149/1945-7111/ac6d13. +T. Baumhöfer et al., "Production caused variation in capacity aging trend and correlation to initial cell performance," Journal of Power Sources, vol. 247, p. 332-338, 2014, doi: 10.1016/j.jpowsour.2013.08.108. diff --git a/examples/data/Baumhofer2014/baumhofer.mat b/examples/data/Baumhofer2014/baumhofer.mat new file mode 100644 index 000000000..eb2a77314 Binary files /dev/null and b/examples/data/Baumhofer2014/baumhofer.mat differ diff --git a/examples/data/Baumhofer2014/read_dataset.py b/examples/data/Baumhofer2014/read_dataset.py new file mode 100644 index 000000000..a61a8b0a9 --- /dev/null +++ b/examples/data/Baumhofer2014/read_dataset.py @@ -0,0 +1,39 @@ +from pybamm import citations +from scipy.io import loadmat + +from pybop import Datasets + +citations.register("""@article{ + Baumhöfer2014, + title={{Production caused variation in capacity aging trend and correlation to initial cell performance}}, + author={Baumhöfer, T and Brühl, M and Rothgang, S and Sauer, D}, + journal={Journal of Power Sources}, + volume={247}, + pages={332-338}, + year={2014}, + doi={10.1016/j.jpowsour.2013.08.108} +}""") +citations.register("""@article{ + Attia2022, + title={{Review—“Knees” in Lithium-Ion Battery Aging Trajectories}}, + author={Attia, P and Bills, A and Brosa Planella, F and Dechent, P and dos Reis, G and Dubarry, M and Gasper, P and Gilchrist, R and Greenbank, S and Howey, D and Liu, O and Khoo, E and Preger, Y and Soni, A and Sripad, S and Stefanopoulou, A and Sulzer, V}, + journal={Journal of The Electrochemical Society}, + volume={169}, + pages={060517}, + year={2022}, + doi={10.1149/1945-7111/ac6d13} +}""") + +matlab_degradation_data = loadmat("../../data/Baumhofer2014/baumhofer.mat")["lifetime"][ + 0 +][0] +degradation_data = Datasets( + [ + { + "Time [s]": dataset[0][0][0], + "Capacity fade": dataset[0][0][1] / dataset[0][0][1][0], + } + for dataset in matlab_degradation_data + ], + domain="Time [s]", +) diff --git a/examples/data/Gunther2025/README.md b/examples/data/Gunther2025/README.md new file mode 100644 index 000000000..0e4358fe1 --- /dev/null +++ b/examples/data/Gunther2025/README.md @@ -0,0 +1,11 @@ +# High temporal resolution GITT as impedance spectra for various electrode thicknesses + +This repository contains a subset of the data published in the research paper: +**J. Günther et al., "Determination of the Solid State Diffusion Coefficient of Li-ion Battery Single-Phase Materials Using Thin Model Electrodes,"** *Journal of The Electrochemical Society*, vol. 172, p. 030525, March 2025. [DOI: 10.1149/1945-7111/adbfc5](https://doi.org/10.1149/1945-7111/adbfc5) + +## Citation + +If you use this dataset in your work, please cite the original publication as follows: + +```plaintext +J. Günther et al., "Determination of the Solid State Diffusion Coefficient of Li-ion Battery Single-Phase Materials Using Thin Model Electrodes," Journal of The Electrochemical Society, vol. 172, p. 030525, March 2025, doi: 10.1149/1945-7111/adbfc5. diff --git a/examples/data/Gunther2025/drt_finetuning.json b/examples/data/Gunther2025/drt_finetuning.json new file mode 100644 index 000000000..d61b0bf48 --- /dev/null +++ b/examples/data/Gunther2025/drt_finetuning.json @@ -0,0 +1,426 @@ +{ + "monolayer_17_microns": { + "lambda": [ + 1e-7, + 1e-10, + 1e-7, + 1e-7, + 1e-7, + 1e-9, + 1e-7, + 1e-9, + 1e-8, + 1e-7, + 1e-7, + 1e-8, + 1e-7, + 1e-7, + 1e-7, + 1e-9, + 1e-8, + 1e-8, + 1e-7, + 1e-7 + ], + "subsampling": [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1 + ], + "start_frequency": [ + 2e-4, + 1e-4, + 3.7e-4, + 2e-4, + 3e-4, + 1e-4, + 2e-4, + 1e-4, + 1e-4, + 2e-4, + 1e-4, + 2e-4, + 2e-4, + 1e-4, + 3.5e-4, + 2.5e-4, + 2.5e-4, + 3.3e-4, + 2e-4, + 1e-4 + ], + "end_frequency": [ + 10, + 80, + 13, + 50, + 24, + 18, + 28, + 27, + 35, + 40, + 50, + 30, + 38, + 33, + 42, + 36, + 41, + 43, + 40, + 48 + ] + }, + "porous_42_microns": { + "lambda": [ + 1e-7, + 1e-7, + 1e-8, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-9, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-9 + ], + "subsampling": [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1 + ], + "start_frequency": [ + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1.4e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1.7e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4 + ], + "end_frequency": [ + 100, + 100, + 100, + 100, + 100, + 63, + 100, + 100, + 100, + 100, + 100, + 100, + 90, + 100, + 100, + 100, + 100, + 100, + 100, + 74 + ] + }, + "porous_80_microns": { + "lambda": [ + 1e-7, + 1e-7, + 1e-7, + 1e-8, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7, + 1e-7 + ], + "subsampling": [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1 + ], + "start_frequency": [ + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4, + 1e-4 + ], + "end_frequency": [ + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100, + 100 + ] + }, + "18650_LG_3500_MJ1_EIS_anode_discharge": { + "lambda": [ + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7, + 6e-7 + ], + "subsampling": [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1 + ], + "start_frequency": [ + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1, + 1e1 + ], + "end_frequency": [ + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4, + 1e4 + ] + } +} diff --git a/examples/data/Gunther2025/impedance_ocv_alignment.json b/examples/data/Gunther2025/impedance_ocv_alignment.json new file mode 100644 index 000000000..34e3dfc09 --- /dev/null +++ b/examples/data/Gunther2025/impedance_ocv_alignment.json @@ -0,0 +1,284 @@ +{ + "monolayer_17_microns": { + "OCV [V]": [ + 3.621560018178841, + 3.6229118029411365, + 3.6657278211215507, + 3.666774555074636, + 3.7045651160702247, + 3.7054549112369175, + 3.737331367756006, + 3.7378255240682075, + 3.768714668323144, + 3.7693507589594786, + 3.8137058778558295, + 3.815089505533486, + 3.8857812784249113, + 3.8888559124200586, + 3.976965800763847, + 3.980530058887571, + 4.072604760517541, + 4.077257864308728, + 4.175348707141178, + 4.1810006718769746 + ], + "indices": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20 + ], + "Positive electrode SOC [-]": [ + 0.8836622861764264, + 0.8833279812900964, + 0.798401424935284, + 0.7980829975007507, + 0.7131503803880486, + 0.7125810964697968, + 0.6274723181163105, + 0.6266905932570239, + 0.5408726363053408, + 0.5399036578572781, + 0.4539747255530346, + 0.4529617979498503, + 0.3666998205688586, + 0.36555549712711716, + 0.27922773039078513, + 0.27824869952955245, + 0.19226394547893982, + 0.19140665836093645, + 0.10570026261756424, + 0.10488193039919773 + ] + }, + "porous_42_microns": { + "OCV [V]": [ + 3.6140400778604542, + 3.615398884268486, + 3.659232562543624, + 3.6603427657700145, + 3.699105976193917, + 3.7002533400009763, + 3.7322753163854476, + 3.733288825392477, + 3.763783921574187, + 3.765012568487434, + 3.809369764543145, + 3.811459509343397, + 3.888241368563032, + 3.8914558377798354, + 3.990111291865414, + 3.993251389713217, + 4.099423500461939, + 4.102722801182583, + 4.2179125214030915, + 4.222934209922096 + ], + "indices": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20 + ], + "Positive electrode SOC [-]": [ + 0.8875533527990672, + 0.8874754860987171, + 0.800662386066443, + 0.800582521914536, + 0.7138308923161695, + 0.7137399964335671, + 0.6270006758762934, + 0.626906005749099, + 0.5401145004157025, + 0.5400513123486737, + 0.4532353062883149, + 0.4531567155282638, + 0.3663989172497242, + 0.36630981739505175, + 0.2795531848192701, + 0.2794621301857168, + 0.19270148762859085, + 0.1926062254093886, + 0.10589332726948768, + 0.10579895330888094 + ] + }, + "porous_80_microns": { + "OCV [V]": [ + 3.614923365081911, + 3.6162500869813106, + 3.660039936115718, + 3.6614331343070217, + 3.7000277586677313, + 3.701447494832104, + 3.733290779545403, + 3.7344457996195923, + 3.765297953270821, + 3.7666232170142795, + 3.812421802633573, + 3.8148323576358227, + 3.8939963654383987, + 3.8976954739936573, + 3.998060214097357, + 4.001197693185197, + 4.109576273650417, + 4.1132705356666515, + 4.231066838410738, + 4.235770481713148 + ], + "indices": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20 + ], + "Positive electrode SOC [-]": [ + 0.8860897422198584, + 0.8859966464294188, + 0.7993479343820302, + 0.7992638834665555, + 0.7126291600267622, + 0.7125459592508565, + 0.6259218772073403, + 0.6258408323725033, + 0.5392102273559071, + 0.5391244837733075, + 0.4524919173737113, + 0.45240559769927735, + 0.3657868220295677, + 0.3657098738581964, + 0.2790871116415675, + 0.2790092976571204, + 0.1923878300091847, + 0.19231396005310936, + 0.10570635952524471, + 0.10564100567049725 + ] + }, + "18650_LG_3500_MJ1_EIS_anode_discharge": { + "indices": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 15, + 20, + 25, + 30, + 35, + 40, + 45, + 50, + 55, + 60, + 65, + 70, + 75, + 80, + 85, + 90, + 91, + 92, + 93, + 94, + 95, + 96, + 97, + 98, + 99, + 0 + ], + "Negative electrode SOC [-]": [ + 0.009972618970932351, + 0.019940387371970805, + 0.029908155773009256, + 0.039875924174047714, + 0.04984369257508616, + 0.059811460976124615, + 0.06977922937716308, + 0.07974699777820153, + 0.08971476617923997, + 0.09968253458027843, + 0.14952137658547068, + 0.19936021859066297, + 0.24919906059585523, + 0.2990379026010475, + 0.34887674460623974, + 0.398715586611432, + 0.4485544286166243, + 0.49839327062181654, + 0.5482321126270089, + 0.5980709546322011, + 0.6479097966373935, + 0.6977486386425856, + 0.7475874806477779, + 0.7974263226529702, + 0.8472651646581624, + 0.8971040066633547, + 0.9070717750643932, + 0.9170395434654317, + 0.9270073118664701, + 0.9369750802675084, + 0.9469428486685469, + 0.9569106170695854, + 0.9668783854706239, + 0.9768461538716623, + 0.9868139222727007, + 0.9967816922727007 + ] + } +} diff --git a/examples/data/Gunther2025/monolayer_17_microns.parquet b/examples/data/Gunther2025/monolayer_17_microns.parquet new file mode 100644 index 000000000..e9d90f676 Binary files /dev/null and b/examples/data/Gunther2025/monolayer_17_microns.parquet differ diff --git a/examples/data/Gunther2025/parameters.py b/examples/data/Gunther2025/parameters.py new file mode 100644 index 000000000..228a16c95 --- /dev/null +++ b/examples/data/Gunther2025/parameters.py @@ -0,0 +1,674 @@ +from copy import deepcopy +from math import log + +from ep_bolfi.utility.preprocessing import SubstitutionDict + + +def electrode_length(par): + return par["Current collector perpendicular area [m2]"] ** 0.5 + + +def cell_volume(par): + return ( + par["Negative electrode thickness [m]"] + + par["Separator thickness [m]"] + + par["Positive electrode thickness [m]"] + ) * par["Current collector perpendicular area [m2]"] + + +def c_max(par): + return ( + 3600 + * par["Nominal cell capacity [A.h]"] + / ( + (1 - par["Positive electrode porosity"]) + * par["Positive electrode thickness [m]"] + * par["Current collector perpendicular area [m2]"] + * 96485.33212 + ) + ) + + +def c_init(par): + return par["Positive electrode SOC"] * c_max(par) + + +def eff_surf_area(par): + return ( + 3 + * (1 - par["Positive electrode porosity"]) + / par["Positive particle radius [m]"] + ) + + +def sep_brugg(par): + return -log(par["Separator MacMullin number"]) / log(par["Separator porosity"]) + + +def OCV_monolayer_17_microns(SOC): + return ( + (SOC < 0.006015075376884422) + * ((4327.768022488337 * SOC + -70.21489642345898) * SOC + 4.791466455864839) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * ((866.6187041254234 * SOC + -28.576748342249516) * SOC + 4.666238156233664) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * ((205.0113975043605 * SOC + -13.9814916785989) * SOC + 4.585744215588657) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * ((71.99050034119864 * SOC + -8.378571276480159) * SOC + 4.52674461909298) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * ((20.874853114468806 * SOC + -4.687456449102825) * SOC + 4.460109794885175) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * ((6.878569070234789 * SOC + -2.834460974583237) * SOC + 4.398779369292598) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * ((4.41457838775159 * SOC + -2.3599632311958203) * SOC + 4.375935521903028) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * ((1.3772369744866069 * SOC + -1.5922637417263275) * SOC + 4.327425787077231) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * ((-0.13856336977755745 * SOC + -1.0114913143456477) * SOC + 4.271795668361062) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * ((-0.49515533896368424 * SOC + -0.8319445699812036) * SOC + 4.249194883572354) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * ((1.7202387822025003 * SOC + -2.2362817901251333) * SOC + 4.471747048110359) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * ((2.3146820951566838 * SOC + -2.6846473800134163) * SOC + 4.556293252018747) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * ((1.1823279856151316 * SOC + -1.682906618363461) * SOC + 4.334744942815917) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * ((0.5976410928342943 * SOC + -1.0952874767435787) * SOC + 4.187103418838845) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * ((0.34566169395441193 * SOC + -0.809188294535943) * SOC + 4.1058936627797635) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * ( + (-0.05691285989067296 * SOC + -0.30364771234104637) * SOC + + 3.947183134627352 + ) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * ((-0.5616711739030507 * SOC + 0.4567821018514451) * SOC + 3.660781956944817) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * ((-1.1586957303367171 * SOC + 1.4340602986355862) * SOC + 3.2608517264302463) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * ((-3.023761075088487 * SOC + 4.711495579685561) * SOC + 1.8210118232734658) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * ((-4.426761993694981 * SOC + 7.261389822078854) * SOC + 0.6624308556570213) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * ((-15.156394830280078 * SOC + 27.407758210973043) * SOC + -8.794466818351793) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * ((-57.37699159733711 * SOC + 109.22363684862648) * SOC + -48.430559802676726) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * ((-278.97535136616716 * SOC + 545.3091648009504) * SOC + -262.97477832468576) + + (SOC >= 0.9939849246231156) + * ((-3112.093210115854 * SOC + 6177.462047356385) * SOC + -3062.1123075410433) + ) + + +def derivative_OCV_monolayer_17_microns(SOC): + return ( + (SOC < 0.006015075376884422) * (8655.536044976674 * SOC + -70.21489642345898) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * (1733.2374082508468 * SOC + -28.576748342249516) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * (410.022795008721 * SOC + -13.9814916785989) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * (143.98100068239728 * SOC + -8.378571276480159) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * (41.74970622893761 * SOC + -4.687456449102825) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * (13.757138140469579 * SOC + -2.834460974583237) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * (8.82915677550318 * SOC + -2.3599632311958203) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * (2.7544739489732137 * SOC + -1.5922637417263275) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * (-0.2771267395551149 * SOC + -1.0114913143456477) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * (-0.9903106779273685 * SOC + -0.8319445699812036) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * (3.4404775644050005 * SOC + -2.2362817901251333) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * (4.6293641903133675 * SOC + -2.6846473800134163) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * (2.3646559712302633 * SOC + -1.682906618363461) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * (1.1952821856685887 * SOC + -1.0952874767435787) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * (0.6913233879088239 * SOC + -0.809188294535943) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * (-0.11382571978134592 * SOC + -0.30364771234104637) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * (-1.1233423478061013 * SOC + 0.4567821018514451) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * (-2.3173914606734343 * SOC + 1.4340602986355862) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * (-6.047522150176974 * SOC + 4.711495579685561) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * (-8.853523987389963 * SOC + 7.261389822078854) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * (-30.312789660560156 * SOC + 27.407758210973043) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * (-114.75398319467422 * SOC + 109.22363684862648) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * (-557.9507027323343 * SOC + 545.3091648009504) + + (SOC >= 0.9939849246231156) * (-6224.186420231708 * SOC + 6177.462047356385) + ) + + +def OCV_porous_42_microns(SOC): + return ( + (SOC < 0.006015075376884422) + * ((5177.783933267361 * SOC + -83.33697942297869) * SOC + 4.942932375182647) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * ((972.3053138390751 * SOC + -32.74443763950444) * SOC + 4.7907733990147605) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * ((254.3341569552358 * SOC + -16.90577744493163) * SOC + 4.703421994172834) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * ((78.56027100125903 * SOC + -9.502075374246829) * SOC + 4.625459895232533) + + (SOC >= 0.03610552763819096) + * (SOC < 0.05115075376884422) + * ((30.806688928820222 * SOC + -6.0537388195685935) * SOC + 4.5632078898421735) + + (SOC >= 0.05115075376884422) + * (SOC < 0.06619597989949749) + * ((20.462287352613203 * SOC + -4.9954909437476545) * SOC + 4.536142801580912) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * ((8.979508757493477 * SOC + -3.4752633816018488) * SOC + 4.485826325007682) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * ((5.022566084589471 * SOC + -2.713263597124353) * SOC + 4.449141204730351) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * ((1.734345170779477 * SOC + -1.8821533683663318) * SOC + 4.396624644069469) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * ((-0.02434527129776143 * SOC + -1.2083186374788113) * SOC + 4.332080419331714) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * ((-0.5164460376261104 * SOC + -0.960542192330081) * SOC + 4.300891092383221) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * ((2.1073903770014795 * SOC + -2.6237907770519087) * SOC + 4.564474202896157) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * ((2.8048070796405113 * SOC + -3.149825210382005) * SOC + 4.663666057637371) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * ((1.3181215064059302 * SOC + -1.8346239620119604) * SOC + 4.372791787573803) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * ((0.6754236712576471 * SOC + -1.18870294877604) * SOC + 4.2105016986225365) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * ((0.35659474200969044 * SOC + -0.8267023394902253) * SOC + 4.1077471889923345) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * ((-0.03396650940428003 * SOC + -0.336247688204395) * SOC + 3.953772771345399) + + (SOC >= 0.7532613065326633) + * (SOC < 0.7883668341708543) + * ((-0.7185926504641884 * SOC + 0.6951570747978622) * SOC + 3.5653141216738504) + + (SOC >= 0.7883668341708543) + * (SOC < 0.8184572864321608) + * ((-0.3241856717609153 * SOC + 0.07328231244810013) * SOC + 3.81044684049607) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * ((-1.2648799442783911 * SOC + 1.6131184757409756) * SOC + 3.1803017766164885) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * ((-3.01427933808759 * SOC + 4.6872967130098) * SOC + 1.829756574083376) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * ((-4.639934862316295 * SOC + 7.6418561907471485) * SOC + 0.4873101613311519) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * ((-16.08006683322219 * SOC + 29.12228488888377) * SOC + -9.595808159094759) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * ((-62.99048692389624 * SOC + 120.02619271745425) * SOC + -53.63464027002556) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * ((-324.8364689617665 * SOC + 635.3154008067868) * SOC + -307.1452783814493) + + (SOC >= 0.9939849246231156) + * ((-4148.740333427573 * SOC + 8237.12098978099) * SOC + -4085.1853560594755) + ) + + +def derivative_OCV_porous_42_microns(SOC): + return ( + (SOC < 0.006015075376884422) * (10355.567866534722 * SOC + -83.33697942297869) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * (1944.6106276781502 * SOC + -32.74443763950444) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * (508.6683139104716 * SOC + -16.90577744493163) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * (157.12054200251805 * SOC + -9.502075374246829) + + (SOC >= 0.03610552763819096) + * (SOC < 0.05115075376884422) + * (61.613377857640444 * SOC + -6.0537388195685935) + + (SOC >= 0.05115075376884422) + * (SOC < 0.06619597989949749) + * (40.924574705226405 * SOC + -4.9954909437476545) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * (17.959017514986954 * SOC + -3.4752633816018488) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * (10.045132169178942 * SOC + -2.713263597124353) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * (3.468690341558954 * SOC + -1.8821533683663318) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * (-0.04869054259552286 * SOC + -1.2083186374788113) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * (-1.0328920752522208 * SOC + -0.960542192330081) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * (4.214780754002959 * SOC + -2.6237907770519087) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * (5.609614159281023 * SOC + -3.149825210382005) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * (2.6362430128118604 * SOC + -1.8346239620119604) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * (1.3508473425152943 * SOC + -1.18870294877604) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * (0.7131894840193809 * SOC + -0.8267023394902253) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * (-0.06793301880856006 * SOC + -0.336247688204395) + + (SOC >= 0.7532613065326633) + * (SOC < 0.7883668341708543) + * (-1.4371853009283768 * SOC + 0.6951570747978622) + + (SOC >= 0.7883668341708543) + * (SOC < 0.8184572864321608) + * (-0.6483713435218306 * SOC + 0.07328231244810013) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * (-2.5297598885567822 * SOC + 1.6131184757409756) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * (-6.02855867617518 * SOC + 4.6872967130098) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * (-9.27986972463259 * SOC + 7.6418561907471485) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * (-32.16013366644438 * SOC + 29.12228488888377) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * (-125.98097384779248 * SOC + 120.02619271745425) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * (-649.672937923533 * SOC + 635.3154008067868) + + (SOC >= 0.9939849246231156) * (-8297.480666855146 * SOC + 8237.12098978099) + ) + + +def OCV_porous_80_microns(SOC): + return ( + (SOC < 0.006015075376884422) + * ((5252.06550238216 * SOC + -85.540902199796) * SOC + 4.981355681431279) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * ((1078.2838141241373 * SOC + -35.32967927673212) * SOC + 4.8303435361073905) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * ((253.07691926793268 * SOC + -17.12536637010527) * SOC + 4.729945378242956) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * ((83.5435749589742 * SOC + -9.984519676649654) * SOC + 4.654751186051314) + + (SOC >= 0.03610552763819096) + * (SOC < 0.05115075376884422) + * ((31.31914813435651 * SOC + -6.213338704439593) * SOC + 4.586670946640943) + + (SOC >= 0.05115075376884422) + * (SOC < 0.06619597989949749) + * ((21.444391680566696 * SOC + -5.203136232649513) * SOC + 4.560834637695333) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * ((9.357225379113515 * SOC + -3.602892597583832) * SOC + 4.507869789944781) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * ((5.197879662078549 * SOC + -2.8019154791506367) * SOC + 4.469308175456582) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * ((1.880057766203663 * SOC + -1.9633234906144423) * SOC + 4.416318854050315) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * ((-0.021088729886969304 * SOC + -1.234907331112197) * SOC + 4.346546469003318) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * ((-0.36835446293730456 * SOC + -1.0600564169404265) * SOC + 4.32453678068903) + + (SOC >= 0.31694974874371856) + * (SOC < 0.34704020100502514) + * ((1.105221645227175 * SOC + -1.9941555714156038) * SOC + 4.472568026845352) + + (SOC >= 0.34704020100502514) + * (SOC < 0.37713065326633166) + * ((2.8107099460094105 * SOC + -3.1779015768458976) * SOC + 4.677971752676939) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * ((2.7632529419825573 * SOC + -3.142106594984284) * SOC + 4.6712220602305194) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * ((1.4898241747671932 * SOC + -2.015563676626641) * SOC + 4.422072092133504) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * ((0.6760747910667533 * SOC + -1.197733278429098) * SOC + 4.2165891223109355) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * ((0.4031517664107582 * SOC + -0.8878545561732949) * SOC + 4.128629501985841) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * ((-0.0423022776961659 * SOC + -0.3284672460513889) * SOC + 3.9530142130099506) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * ((-0.5695955926386773 * SOC + 0.4659120566277579) * SOC + 3.6538266173006377) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * ((-1.1877907988101697 * SOC + 1.4778447984845116) * SOC + 3.239714754324609) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * ((-3.1857782388037776 * SOC + 4.988860938137805) * SOC + 1.6972583196459254) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * ((-4.599262662306728 * SOC + 7.55780850289193) * SOC + 0.5300202021401219) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * ((-16.67280669139791 * SOC + 30.22755587208667) * SOC + -10.11137565778472) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * ((-64.08420926520012 * SOC + 122.10227712296182) * SOC + -54.620522965688) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * ((-352.042942590193 * SOC + 688.7790177879579) * SOC + -333.4126650747057) + + (SOC >= 0.9939849246231156) + * ((-4527.655695340203 * SOC + 8989.771272383106) * SOC + -4458.943245315109) + ) + + +def derivative_OCV_porous_80_microns(SOC): + return ( + (SOC < 0.006015075376884422) * (10504.13100476432 * SOC + -85.540902199796) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * (2156.5676282482746 * SOC + -35.32967927673212) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * (506.15383853586536 * SOC + -17.12536637010527) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * (167.0871499179484 * SOC + -9.984519676649654) + + (SOC >= 0.03610552763819096) + * (SOC < 0.05115075376884422) + * (62.63829626871302 * SOC + -6.213338704439593) + + (SOC >= 0.05115075376884422) + * (SOC < 0.06619597989949749) + * (42.88878336113339 * SOC + -5.203136232649513) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * (18.71445075822703 * SOC + -3.602892597583832) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * (10.395759324157098 * SOC + -2.8019154791506367) + + (SOC >= 0.12637688442211056) + * (SOC < 0.19157286432160803) + * (3.760115532407326 * SOC + -1.9633234906144423) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * (-0.04217745977393861 * SOC + -1.234907331112197) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * (-0.7367089258746091 * SOC + -1.0600564169404265) + + (SOC >= 0.31694974874371856) + * (SOC < 0.34704020100502514) + * (2.21044329045435 * SOC + -1.9941555714156038) + + (SOC >= 0.34704020100502514) + * (SOC < 0.37713065326633166) + * (5.621419892018821 * SOC + -3.1779015768458976) + + (SOC >= 0.37713065326633166) + * (SOC < 0.44232663316582915) + * (5.526505883965115 * SOC + -3.142106594984284) + + (SOC >= 0.44232663316582915) + * (SOC < 0.5025075376884423) + * (2.9796483495343864 * SOC + -2.015563676626641) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * (1.3521495821335066 * SOC + -1.197733278429098) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * (0.8063035328215165 * SOC + -0.8878545561732949) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * (-0.0846045553923318 * SOC + -0.3284672460513889) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * (-1.1391911852773546 * SOC + 0.4659120566277579) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8786381909547738) + * (-2.3755815976203394 * SOC + 1.4778447984845116) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * (-6.371556477607555 * SOC + 4.988860938137805) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * (-9.198525324613456 * SOC + 7.55780850289193) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * (-33.34561338279582 * SOC + 30.22755587208667) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * (-128.16841853040023 * SOC + 122.10227712296182) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * (-704.085885180386 * SOC + 688.7790177879579) + + (SOC >= 0.9939849246231156) * (-9055.311390680406 * SOC + 8989.771272383106) + ) + + +parameters = SubstitutionDict( + { + ######################################## + # Placeholders for unknown parameters. # + ######################################## + "Positive particle diffusivity [m2.s-1]": 1e-15, + "Positive electrode exchange-current density [A.m-2]": 1.0, + "Positive electrode double-layer capacity [F.m-2]": 0.24, + "Positive electrode Bruggeman coefficient": 1.5, + "Thermodynamic factor": 2, + ################################################### + # Lithium-metal electrode parameters from Xu2019. # + ################################################### + "Exchange-current density for lithium metal electrode [A.m-2]": ( + lambda c_e, c_Li, T: 3.5e-8 * 96458.33212 * c_Li**0.7 + c_e**0.3 + ), + "Lithium metal partial molar volume [m3.mol-1]": 1.3e-05, + ######################################################## + # Parameters taken from the corresponding Excel sheet. # + # Note that these are for MEL; substitute as needed. # + ######################################################## + "Electrolyte diffusivity [m2.s-1]": 5.66e-10, + "Cation transference number": 0.4665, + "Positive electrode thickness [m]": 1.7e-5, + "Positive particle radius [m]": 5.22e-6, + "Current collector perpendicular area [m2]": 0.000254469, + "Positive electrode conductivity [S.m-1]": 0.5, # 100 + "Electrolyte conductivity [S.m-1]": 0.59, + "Positive electrode porosity": 0.75, + "Typical electrolyte concentration [mol.m-3]": 250, + "Initial concentration in electrolyte [mol.m-3]": 250, + "Lower voltage cut-off [V]": 3.4, + "Upper voltage cut-off [V]": 4.3, + "Reference temperature [K]": 298.15, + "Ambient temperature [K]": 298.15, + "Initial temperature [K]": 298.15, + "Nominal cell capacity [A.h]": 7e-4, + "Positive electrode OCP [V]": OCV_monolayer_17_microns, + "Positive electrode OCP derivative by SOC [V]": ( + derivative_OCV_monolayer_17_microns + ), + ##################################### + # Preliminary separator properties. # + ##################################### + "Separator thickness [m]": 176e-6, + "Separator MacMullin number": 3.472, + "Separator porosity": 0.63, + ################################################################### + # Dummy parameters for simulations with a negative electrode. # + # These make the negative electrode have approximately no effect. # + ################################################################### + "Negative particle diffusivity [m2.s-1]": 1e-8, + "Negative electrode thickness [m]": 1e-8, + # This is just the inverse of the thickness for a planar electrode. + "Negative electrode surface area to volume ratio [m-1]": 1e8, + "Negative particle radius [m]": 1e-4, + "Maximum concentration in negative electrode [mol.m-3]": 1, + "Negative electrode conductivity [S.m-1]": 1e3, + "Negative electrode OCP [V]": 1e-6, + "Negative electrode OCP derivative by SOC [V]": 0, + "Negative electrode Bruggeman coefficient (electrolyte)": 1, + "Negative electrode Bruggeman coefficient (electrode)": 1, + "Negative electrode porosity": 0.99, + "Negative electrode exchange-current density [A.m-2]": 1000, + "Initial concentration in negative electrode [mol.m-3]": 1, + "Negative electrode double-layer capacity [F.m-2]": 0.6, + ################################# + # Boilerplate PyBaMM variables. # + ################################# + "Typical current [A]": 1, + "Current function [A]": 0, + "Negative electrode OCP entropic change [V.K-1]": 0, + "Negative electrode OCP entropic change partial derivative by SOC [V.K-1]": 0, + # Not relevant, as we only ever evaluate the model at reference T. + "Positive electrode OCP entropic change [V.K-1]": 0, + "Positive electrode OCP entropic change partial derivative by SOC [V.K-1]": 0, + "Number of electrodes connected in parallel to make a cell": 1, + "Number of cells connected in series to make a battery": 1, + "Negative electrode anodic charge-transfer coefficient": 0.5, + "Negative electrode cathodic charge-transfer coefficient": 0.5, + "Positive electrode anodic charge-transfer coefficient": 0.5, + "Positive electrode cathodic charge-transfer coefficient": 0.5, + "Negative electrode active material volume fraction": 1, + "Positive electrode active material volume fraction": 1, + "Negative electrode electrons in reaction": 1, + "Positive electrode electrons in reaction": 1, + }, + substitutions={ + "Electrode width [m]": electrode_length, + "Electrode height [m]": electrode_length, + "Cell volume [m3]": cell_volume, + "Positive electrode Bruggeman coefficient (electrolyte)": ( + "Positive electrode Bruggeman coefficient" + ), + "Positive electrode Bruggeman coefficient (electrode)": ( + "Positive electrode Bruggeman coefficient" + ), + "Separator Bruggeman coefficient": sep_brugg, + "Separator Bruggeman coefficient (electrolyte)": ( + "Separator Bruggeman coefficient" + ), + "Separator Bruggeman coefficient (electrode)": ( + "Separator Bruggeman coefficient" + ), + "Maximum concentration in positive electrode [mol.m-3]": c_max, + "Initial concentration in positive electrode [mol.m-3]": c_init, + "Positive electrode surface area to volume ratio [m-1]": eff_surf_area, + }, +) + + +def get_parameters_for_cell(cell_name): + completed_parameters = deepcopy(parameters) + match cell_name: + case "monolayer_17_microns": + completed_parameters.update( + { + "Nominal cell capacity [A.h]": 7e-4, + "Positive electrode thickness [m]": 1.7e-5, + "Positive electrode porosity": 0.75, + "Positive electrode OCP [V]": OCV_monolayer_17_microns, + "Positive electrode OCP derivative by SOC [V]": ( + derivative_OCV_monolayer_17_microns + ), + } + ) + case "porous_42_microns": + completed_parameters.update( + { + "Nominal cell capacity [A.h]": 4e-3, + "Positive electrode thickness [m]": 4.2e-5, + "Positive electrode porosity": 0.3102, + "Positive electrode OCP [V]": OCV_porous_42_microns, + "Positive electrode OCP derivative by SOC [V]": ( + derivative_OCV_porous_42_microns + ), + } + ) + case "porous_80_microns": + completed_parameters.update( + { + "Nominal cell capacity [A.h]": 8e-3, + "Positive electrode thickness [m]": 8e-5, + "Positive electrode porosity": 0.271, + "Positive electrode OCP [V]": OCV_porous_80_microns, + "Positive electrode OCP derivative by SOC [V]": ( + derivative_OCV_porous_80_microns + ), + } + ) + case _: + raise ValueError( + "cell_name must be one of 'monolayer_17_microns', " + "'porous_42_microns', or 'porous_80_microns'." + ) + return completed_parameters diff --git a/examples/data/Gunther2025/porous_42_microns.parquet b/examples/data/Gunther2025/porous_42_microns.parquet new file mode 100644 index 000000000..4d1e5f6ab Binary files /dev/null and b/examples/data/Gunther2025/porous_42_microns.parquet differ diff --git a/examples/data/Gunther2025/porous_80_microns.parquet b/examples/data/Gunther2025/porous_80_microns.parquet new file mode 100644 index 000000000..92da3f9c3 Binary files /dev/null and b/examples/data/Gunther2025/porous_80_microns.parquet differ diff --git a/examples/data/Kuhn2026/18650_LG_3500_MJ1_EIS_anode_discharge.parquet b/examples/data/Kuhn2026/18650_LG_3500_MJ1_EIS_anode_discharge.parquet new file mode 100644 index 000000000..82d3bd4d2 Binary files /dev/null and b/examples/data/Kuhn2026/18650_LG_3500_MJ1_EIS_anode_discharge.parquet differ diff --git a/examples/data/Kuhn2026/README.md b/examples/data/Kuhn2026/README.md new file mode 100644 index 000000000..6e6982e4e --- /dev/null +++ b/examples/data/Kuhn2026/README.md @@ -0,0 +1,16 @@ +# LG MJ1 impedance spectra + +This repository contains a subset of the data published in the research paper: +**Y. Kuhn et al., "Workflows and Principles for Collaboration and Communication in Battery Research,"** *Digital Discovery*, Advance Article, 2026. [DOI: 10.1039/D5DD00247H](https://doi.org/10.1039/D5DD00247H) + +## Accessing the full dataset + +The full dataset can be accessed through the following link: +[GITT and impedance data on Zenodo](https://doi.org/10.5281/zenodo.15407848) + +## Citation + +If you use this dataset in your work, please cite the original publication as follows: + +```plaintext +Y. Kuhn et al., "Workflows and Principles for Collaboration and Communication in Battery Research," Digital Discovery, Advance Article, 2026, doi: 10.1039/D5DD00247H. diff --git a/examples/data/Kuhn2026/parameters.py b/examples/data/Kuhn2026/parameters.py new file mode 100644 index 000000000..a7b30b6ee --- /dev/null +++ b/examples/data/Kuhn2026/parameters.py @@ -0,0 +1,1364 @@ +"""!@file +Parameter file for the 18650 LG MJ1 3500 cell with a third reference electrode. +""" + +from copy import deepcopy + +from ep_bolfi.utility.preprocessing import SubstitutionDict +from numpy import exp, log +from pybamm import Parameter + +############################## +# Reaction symmetry factors. # +############################## + +"""! Anodic symmetry factor at Graphite-SiOx. """ +αₙₙ = 0.5 +"""! Cathodic symmetry factor at Graphite-SiOx. """ +αₚₙ = 0.5 +"""! Anodic symmetry factor at NMC. """ +αₙₚ = 0.5 +"""! Cathodic symmetry factor at NMC. """ +αₚₚ = 0.5 + +############# +# Geometry. # +############# + + +def sqrt_area(par): + return par["Current collector perpendicular area [m2]"] ** 0.5 + + +def cell_volume(par): + return par["Current collector perpendicular area [m2]"] * ( + par["Negative electrode thickness [m]"] + + par["Separator thickness [m]"] + + par["Positive electrode thickness [m]"] + ) + + +################################################### +# Maximum and initial active material lithiation. # +################################################### + + +def negative_concentration_correction(par): + return 43.9445 / ( + par["Negative electrode active material volume fraction"] + * par["Negative electrode thickness [m]"] + * par["Current collector perpendicular area [m2]"] + * 96485.33212 + ) + + +def positive_concentration(par): + # Use the areal capacity of 5.26 mAh/cm². Since it applies to the double- + # coated current collector, half it. + return ( + 3600.0 + * 26.3 # areal capacity in C/m² + * par["Current collector perpendicular area [m2]"] + / 96485.33212 # scale to mol instead of C + ) / ( + par["Positive electrode active material volume fraction"] + * par["Positive electrode thickness [m]"] + * par["Current collector perpendicular area [m2]"] + ) + + +def placeholder_initial_negative_concentration(par): + return 0.5 * par["Maximum concentration in negative electrode [mol.m-3]"] + + +def placeholder_initial_positive_concentration(par): + return 0.5 * par["Maximum concentration in positive electrode [mol.m-3]"] + + +#################################### +# Active material volume fraction. # +#################################### + + +def negative_volume_fraction(par): + return 0.93 * (1 - par["Negative electrode porosity"]) + + +def positive_volume_fraction(par): + return 0.92 * (1 - par["Positive electrode porosity"]) + + +######################################################################### +# Bruggeman coefficients, expressed in terms of path-length tortuosity. # +######################################################################### + + +def negative_bruggeman(par): + return 1 - 2 * log(par["Negative electrode tortuosity square"] ** 0.5) / log( + par["Negative electrode porosity"] + ) + + +def positive_bruggeman(par): + return 1 - 2 * log(par["Positive electrode tortuosity square"] ** 0.5) / log( + par["Positive electrode porosity"] + ) + + +###################################################################### +# Standard relation for 1D+1D models to ensure lithium conservation. # +###################################################################### + + +def negative_effective_surface_area(par): + return ( + 3 + * (1 - par["Negative electrode porosity"]) + / par["Negative particle radius [m]"] + ) + + +def positive_effective_surface_area(par): + return ( + 3 + * (1 - par["Positive electrode porosity"]) + / par["Positive particle radius [m]"] + ) + + +################################## +# Open Circuit Potential curves. # +################################## + + +# Discharge profile - Delithiation direction of the Anode +def _18650_LG_3500_MJ1_Anode_GITT_dch_ohne_EIS_2320__extraction(SOC): + return ( + (SOC < 0.006015075376884422) + * ((5008.957495658953 * SOC + -81.88347557587394) * SOC + 1.1056235451272727) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * ((1017.5623642237479 * SOC + -33.86639042684935) * SOC + 0.9612103518524422) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * ((234.26254027875802 * SOC + -16.58656013982172) * SOC + 0.8659107853197145) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * ((81.70919824890336 * SOC + -10.160921381559184) * SOC + 0.7982478405059502) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * ((20.66446523763034 * SOC + -5.752816791750165) * SOC + 0.7186693694562571) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * ((1.6458973973553839 * SOC + -3.23491132280601) * SOC + 0.6353317595507268) + + (SOC >= 0.09628643216080401) + * (SOC < 0.11133165829145729) + * ((0.4299162758333068 * SOC + -3.0007463552735203) * SOC + 0.6240583049203501) + + (SOC >= 0.11133165829145729) + * (SOC < 0.12637688442211056) + * ((43.272535809588476 * SOC + -12.540226111759424) * SOC + 1.1550813551840307) + + (SOC >= 0.12637688442211056) + * (SOC < 0.14643718592964824) + * ((22.28677763548575 * SOC + -7.235996641201524) * SOC + 0.8199153578095064) + + (SOC >= 0.14643718592964824) + * (SOC < 0.1614824120603015) + * ((5.02482586462844 * SOC + -2.1804133592462023) * SOC + 0.44975266328825025) + + (SOC >= 0.1614824120603015) + * (SOC < 0.19157286432160803) + * ((3.107303716296144 * SOC + -1.561121155862704) * SOC + 0.3997502639019974) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * ( + (-0.1348854947984961 * SOC + -0.3188902081786793) * SOC + + 0.28076139350360996 + ) + + (SOC >= 0.2517537688442211) + * (SOC < 0.28685929648241204) + * ((-4.393207820206399 * SOC + 1.8252071805711765) * SOC + 0.01086909431019567) + + (SOC >= 0.28685929648241204) + * (SOC < 0.30190452261306533) + * ((-17.74875373087525 * SOC + 9.487532188717267) * SOC + -1.0881354863179844) + + (SOC >= 0.30190452261306533) + * (SOC < 0.31694974874371856) + * ((30.384751611729655 * SOC + -19.575913715587845) * SOC + 3.299057394296952) + + (SOC >= 0.31694974874371856) + * (SOC < 0.34704020100502514) + * ((3.2542005829125884 * SOC + -2.377871051863501) * SOC + 0.5735997437213527) + + (SOC >= 0.34704020100502514) + * (SOC < 0.37713065326633166) + * ((0.8083175790748953 * SOC + -0.6802315932902943) * SOC + 0.27902517425269924) + + (SOC >= 0.37713065326633166) + * (SOC < 0.5025075376884423) + * ( + (0.05667695122594907 * SOC + -0.11329815128590504) * SOC + + 0.17212118458187764 + ) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * ( + (-0.5279640038688811 * SOC + 0.47427482226713735) * SOC + + 0.024491260505669743 + ) + + (SOC >= 0.5677035175879397) + * (SOC < 0.582748743718593) + * ((-6.403684546701641 * SOC + 7.1456092633268895) * SOC + -1.8691787540919336) + + (SOC >= 0.582748743718593) + * (SOC < 0.5927788944723618) + * ((5.090840704572997 * SOC + -6.251231036317108) * SOC + 2.0343171731161647) + + (SOC >= 0.5927788944723618) + * (SOC < 0.5977939698492463) + * ((-213.7965887884052 * SOC + 253.25246590117195) * SOC + -74.87984011793162) + + (SOC >= 0.5977939698492463) + * (SOC < 0.6028090452261307) + * ((-27.041551689329253 * SOC + 29.970395867572734) * SOC + -8.141502597159956) + + (SOC >= 0.6028090452261307) + * (SOC < 0.607824120603015) + * ((325.22633701847644 * SOC + -394.7301434439819) * SOC + 119.86516070755022) + + (SOC >= 0.607824120603015) + * (SOC < 0.6128391959798994) + * ((-84.21197551394732 * SOC + 103.0028210684236) * SOC + -31.401890017392248) + + (SOC >= 0.6128391959798994) + * (SOC < 0.6278844221105527) + * ((6.337545487311971 * SOC + -7.98177022512953) * SOC + 2.6059638298573695) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * ( + (0.055219990764311166 * SOC + -0.09262159730911579) * SOC + + 0.12922706629572556 + ) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8786381909547738) + * ( + (-0.011749832709947228 * SOC + 0.008269956147848312) * SOC + + 0.09122821460817487 + ) + + (SOC >= 0.8786381909547738) + * (SOC < 0.938819095477387) + * ( + (-0.08775512587896728 * SOC + 0.14183226273387106) * SOC + + 0.03255174288893059 + ) + + (SOC >= 0.938819095477387) + * (SOC < 0.9538643216080401) + * ((-0.583437644148546 * SOC + 1.072544689625488) * SOC + -0.40433355649304303) + + (SOC >= 0.9538643216080401) + * (SOC < 0.9689095477386934) + * ((-0.7849041067989617 * SOC + 1.4568880310711734) * SOC + -0.5876392568193012) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9789396984924623) + * ((-8.958760838734406 * SOC + 17.296343689911737) * SOC + -8.26113916623649) + + (SOC >= 0.9789396984924623) + * (SOC < 0.9839547738693467) + * ((-686.2640629111021 * SOC + 1343.378440086051) * SOC + -657.3383429773813) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * ((-622.0560205052128 * SOC + 1217.0228203938912) * SOC + -595.1742353767207) + + (SOC >= 0.9939849246231156) + * ((-5411.214994336782 * SOC + 10737.726463618066) * SOC + -5326.892181961326) + ) + + +def derivative__18650_LG_3500_MJ1_Anode_GITT_dch_ohne_EIS_2320__extraction(SOC): + return ( + (SOC < 0.006015075376884422) * (10017.914991317906 * SOC + -81.88347557587394) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * (2035.1247284474957 * SOC + -33.86639042684935) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * (468.52508055751605 * SOC + -16.58656013982172) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * (163.41839649780673 * SOC + -10.160921381559184) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * (41.32893047526068 * SOC + -5.752816791750165) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * (3.2917947947107677 * SOC + -3.23491132280601) + + (SOC >= 0.09628643216080401) + * (SOC < 0.11133165829145729) + * (0.8598325516666137 * SOC + -3.0007463552735203) + + (SOC >= 0.11133165829145729) + * (SOC < 0.12637688442211056) + * (86.54507161917695 * SOC + -12.540226111759424) + + (SOC >= 0.12637688442211056) + * (SOC < 0.14643718592964824) + * (44.5735552709715 * SOC + -7.235996641201524) + + (SOC >= 0.14643718592964824) + * (SOC < 0.1614824120603015) + * (10.04965172925688 * SOC + -2.1804133592462023) + + (SOC >= 0.1614824120603015) + * (SOC < 0.19157286432160803) + * (6.214607432592288 * SOC + -1.561121155862704) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * (-0.2697709895969922 * SOC + -0.3188902081786793) + + (SOC >= 0.2517537688442211) + * (SOC < 0.28685929648241204) + * (-8.786415640412798 * SOC + 1.8252071805711765) + + (SOC >= 0.28685929648241204) + * (SOC < 0.30190452261306533) + * (-35.4975074617505 * SOC + 9.487532188717267) + + (SOC >= 0.30190452261306533) + * (SOC < 0.31694974874371856) + * (60.76950322345931 * SOC + -19.575913715587845) + + (SOC >= 0.31694974874371856) + * (SOC < 0.34704020100502514) + * (6.508401165825177 * SOC + -2.377871051863501) + + (SOC >= 0.34704020100502514) + * (SOC < 0.37713065326633166) + * (1.6166351581497906 * SOC + -0.6802315932902943) + + (SOC >= 0.37713065326633166) + * (SOC < 0.5025075376884423) + * (0.11335390245189814 * SOC + -0.11329815128590504) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * (-1.0559280077377622 * SOC + 0.47427482226713735) + + (SOC >= 0.5677035175879397) + * (SOC < 0.582748743718593) + * (-12.807369093403281 * SOC + 7.1456092633268895) + + (SOC >= 0.582748743718593) + * (SOC < 0.5927788944723618) + * (10.181681409145995 * SOC + -6.251231036317108) + + (SOC >= 0.5927788944723618) + * (SOC < 0.5977939698492463) + * (-427.5931775768104 * SOC + 253.25246590117195) + + (SOC >= 0.5977939698492463) + * (SOC < 0.6028090452261307) + * (-54.08310337865851 * SOC + 29.970395867572734) + + (SOC >= 0.6028090452261307) + * (SOC < 0.607824120603015) + * (650.4526740369529 * SOC + -394.7301434439819) + + (SOC >= 0.607824120603015) + * (SOC < 0.6128391959798994) + * (-168.42395102789465 * SOC + 103.0028210684236) + + (SOC >= 0.6128391959798994) + * (SOC < 0.6278844221105527) + * (12.675090974623942 * SOC + -7.98177022512953) + + (SOC >= 0.6278844221105527) + * (SOC < 0.7532613065326633) + * (0.11043998152862233 * SOC + -0.09262159730911579) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8786381909547738) + * (-0.023499665419894455 * SOC + 0.008269956147848312) + + (SOC >= 0.8786381909547738) + * (SOC < 0.938819095477387) + * (-0.17551025175793455 * SOC + 0.14183226273387106) + + (SOC >= 0.938819095477387) + * (SOC < 0.9538643216080401) + * (-1.166875288297092 * SOC + 1.072544689625488) + + (SOC >= 0.9538643216080401) + * (SOC < 0.9689095477386934) + * (-1.5698082135979234 * SOC + 1.4568880310711734) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9789396984924623) + * (-17.917521677468812 * SOC + 17.296343689911737) + + (SOC >= 0.9789396984924623) + * (SOC < 0.9839547738693467) + * (-1372.5281258222042 * SOC + 1343.378440086051) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9939849246231156) + * (-1244.1120410104256 * SOC + 1217.0228203938912) + + (SOC >= 0.9939849246231156) * (-10822.429988673564 * SOC + 10737.726463618066) + ) + + +# Discharge profile - lithiation direction of cathode +def _18650_LG_3500_MJ1_Cathode_GITT_ch_ohne_EIS_2327__extraction(SOC): + return ( + (SOC < 0.006015075376884422) + * ((2943.810072716413 * SOC + -45.47969560495423) * SOC + 4.446944580747672) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * ((634.1988911293374 * SOC + -17.694724907871773) * SOC + 4.363380234203933) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * ((113.15861792600845 * SOC + -6.200419383436554) * SOC + 4.299988272831737) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * ((23.223632167815595 * SOC + -2.4123035511492503) * SOC + 4.26009884204501) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * ((4.859053133897305 * SOC + -1.0861779194034114) * SOC + 4.236158609220656) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * ((0.4366753467279523 * SOC + -0.500690657188386) * SOC + 4.216780157700159) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * ((-0.6822066464878844 * SOC + -0.2852243469172322) * SOC + 4.206406916566717) + + (SOC >= 0.12637688442211056) + * (SOC < 0.1614824120603015) + * ((-4.887988966797366 * SOC + 0.777802985479525) * SOC + 4.139235875404793) + + (SOC >= 0.1614824120603015) + * (SOC < 0.17652763819095477) + * ((-11.823704457352505 * SOC + 3.017795119036691) * SOC + 3.9583762090433083) + + (SOC >= 0.17652763819095477) + * (SOC < 0.19157286432160803) + * ((-3.138277623014801 * SOC + -0.04864065245465099) * SOC + 4.229031541246172) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * ((3.1305873867419223 * SOC + -2.4505295043838373) * SOC + 4.4590999048192685) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * ((0.7365036271785357 * SOC + -1.2450902855862012) * SOC + 4.307362971596806) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * ((-0.6342814441555902 * SOC + -0.3761503177043437) * SOC + 4.1696578193500216) + + (SOC >= 0.37713065326633166) + * (SOC < 0.41223618090452263) + * ((-1.2465376975685558 * SOC + 0.08565088392765574) * SOC + 4.082578124924652) + + (SOC >= 0.41223618090452263) + * (SOC < 0.44232663316582915) + * ((-0.5336389055478321 * SOC + -0.5021144668608031) * SOC + 4.203727196663067) + + (SOC >= 0.44232663316582915) + * (SOC < 0.4724170854271357) + * ((0.8427801315624492 * SOC + -1.7197680638812471) * SOC + 4.473027504629329) + + (SOC >= 0.4724170854271357) + * (SOC < 0.5025075376884423) + * ((2.2869769138767424 * SOC + -3.084294533249931) * SOC + 4.7953403134528685) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * ((1.624528113077531 * SOC + -2.4185235017808964) * SOC + 4.628062832609103) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * ((0.6844382533436715 * SOC + -1.3511388613416102) * SOC + 4.325083825110767) + + (SOC >= 0.6278844221105527) + * (SOC < 0.6930804020100503) + * ((0.12517928333738837 * SOC + -0.6488388709564106) * SOC + 4.1046022133051565) + + (SOC >= 0.6930804020100503) + * (SOC < 0.7532613065326633) + * ((-0.42125376061284214 * SOC + 0.1086051965887691) * SOC + 3.842117393888003) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * ((-1.0456222985791328 * SOC + 1.049230517921501) * SOC + 3.487849064635668) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8485477386934673) + * ((-0.5071542408868481 * SOC + 0.16780430726248596) * SOC + 3.8485539169184904) + + (SOC >= 0.8485477386934673) + * (SOC < 0.8786381909547738) + * ((1.1369428932853225 * SOC + -2.622385503525493) * SOC + 5.032358544153567) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * ((0.6050016717372273 * SOC + -1.687617758334909) * SOC + 4.621697223854881) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * ((-0.904926011592579 * SOC + 1.0566113117188252) * SOC + 3.3748174441032006) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * ((-6.697727169005248 * SOC + 11.933395997483785) * SOC + -1.730849136093184) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * ((-37.178817824355065 * SOC + 71.00023552039784) * SOC + -30.3460615203403) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9889698492462311) + * ((-142.86486270722526 * SOC + 278.9808123081457) * SOC + -132.66780222153466) + + (SOC >= 0.9889698492462311) + * (SOC < 0.9939849246231156) + * ((-315.0640613306459 * SOC + 619.5804433139565) * SOC + -301.08918508613715) + + (SOC >= 0.9939849246231156) + * ((-2264.12836692133 * SOC + 4494.261517070292) * SOC + -2226.7764726042806) + ) + + +def derivative__18650_LG_3500_MJ1_Cathode_GITT_ch_ohne_EIS_2327__extraction(SOC): + return ( + (SOC < 0.006015075376884422) * (5887.620145432826 * SOC + -45.47969560495423) + + (SOC >= 0.006015075376884422) + * (SOC < 0.011030150753768845) + * (1268.3977822586749 * SOC + -17.694724907871773) + + (SOC >= 0.011030150753768845) + * (SOC < 0.02106030150753769) + * (226.3172358520169 * SOC + -6.200419383436554) + + (SOC >= 0.02106030150753769) + * (SOC < 0.03610552763819096) + * (46.44726433563119 * SOC + -2.4123035511492503) + + (SOC >= 0.03610552763819096) + * (SOC < 0.06619597989949749) + * (9.71810626779461 * SOC + -1.0861779194034114) + + (SOC >= 0.06619597989949749) + * (SOC < 0.09628643216080401) + * (0.8733506934559045 * SOC + -0.500690657188386) + + (SOC >= 0.09628643216080401) + * (SOC < 0.12637688442211056) + * (-1.3644132929757689 * SOC + -0.2852243469172322) + + (SOC >= 0.12637688442211056) + * (SOC < 0.1614824120603015) + * (-9.775977933594731 * SOC + 0.777802985479525) + + (SOC >= 0.1614824120603015) + * (SOC < 0.17652763819095477) + * (-23.64740891470501 * SOC + 3.017795119036691) + + (SOC >= 0.17652763819095477) + * (SOC < 0.19157286432160803) + * (-6.276555246029602 * SOC + -0.04864065245465099) + + (SOC >= 0.19157286432160803) + * (SOC < 0.2517537688442211) + * (6.2611747734838445 * SOC + -2.4505295043838373) + + (SOC >= 0.2517537688442211) + * (SOC < 0.31694974874371856) + * (1.4730072543570714 * SOC + -1.2450902855862012) + + (SOC >= 0.31694974874371856) + * (SOC < 0.37713065326633166) + * (-1.2685628883111804 * SOC + -0.3761503177043437) + + (SOC >= 0.37713065326633166) + * (SOC < 0.41223618090452263) + * (-2.4930753951371116 * SOC + 0.08565088392765574) + + (SOC >= 0.41223618090452263) + * (SOC < 0.44232663316582915) + * (-1.0672778110956642 * SOC + -0.5021144668608031) + + (SOC >= 0.44232663316582915) + * (SOC < 0.4724170854271357) + * (1.6855602631248985 * SOC + -1.7197680638812471) + + (SOC >= 0.4724170854271357) + * (SOC < 0.5025075376884423) + * (4.573953827753485 * SOC + -3.084294533249931) + + (SOC >= 0.5025075376884423) + * (SOC < 0.5677035175879397) + * (3.249056226155062 * SOC + -2.4185235017808964) + + (SOC >= 0.5677035175879397) + * (SOC < 0.6278844221105527) + * (1.368876506687343 * SOC + -1.3511388613416102) + + (SOC >= 0.6278844221105527) + * (SOC < 0.6930804020100503) + * (0.25035856667477674 * SOC + -0.6488388709564106) + + (SOC >= 0.6930804020100503) + * (SOC < 0.7532613065326633) + * (-0.8425075212256843 * SOC + 0.1086051965887691) + + (SOC >= 0.7532613065326633) + * (SOC < 0.8184572864321608) + * (-2.0912445971582656 * SOC + 1.049230517921501) + + (SOC >= 0.8184572864321608) + * (SOC < 0.8485477386934673) + * (-1.0143084817736963 * SOC + 0.16780430726248596) + + (SOC >= 0.8485477386934673) + * (SOC < 0.8786381909547738) + * (2.273885786570645 * SOC + -2.622385503525493) + + (SOC >= 0.8786381909547738) + * (SOC < 0.9087286432160804) + * (1.2100033434744546 * SOC + -1.687617758334909) + + (SOC >= 0.9087286432160804) + * (SOC < 0.938819095477387) + * (-1.809852023185158 * SOC + 1.0566113117188252) + + (SOC >= 0.938819095477387) + * (SOC < 0.9689095477386934) + * (-13.395454338010495 * SOC + 11.933395997483785) + + (SOC >= 0.9689095477386934) + * (SOC < 0.9839547738693467) + * (-74.35763564871013 * SOC + 71.00023552039784) + + (SOC >= 0.9839547738693467) + * (SOC < 0.9889698492462311) + * (-285.7297254144505 * SOC + 278.9808123081457) + + (SOC >= 0.9889698492462311) + * (SOC < 0.9939849246231156) + * (-630.1281226612919 * SOC + 619.5804433139565) + + (SOC >= 0.9939849246231156) * (-4528.25673384266 * SOC + 4494.261517070292) + ) + + +##################################################################### +# Exchange-current densities (mean of lithiation and delithiation). # +##################################################################### + +# In order to be pickleable and therefore parellizable, the returned functions +# have to get the parameters via their names rather than directly as numbers. + + +def graphite_exchange_prefactor(cₑ, cₙ, cₙ_max, T): + αₙₙ = Parameter("Anodic symmetry factor at Graphite-SiOx") + αₚₙ = Parameter("Cathodic symmetry factor at Graphite-SiOx") + return Negative_electrode_exchange_current_density_A_per_m2(cₙ / cₙ_max) + + +def graphite_exchange_prefactor_derivative(cₑ, cₙ, cₙ_max, T): + αₙₙ = Parameter("Anodic symmetry factor at Graphite-SiOx") + αₚₙ = Parameter("Cathodic symmetry factor at Graphite-SiOx") + return Negative_electrode_exchange_current_density_A_per_m2(cₙ / cₙ_max) * αₚₙ / cₑ + + +def nmc_exchange_prefactor(cₑ, cₚ, cₚ_max, T): + αₙₚ = Parameter("Anodic symmetry factor at NMC") + αₚₚ = Parameter("Cathodic symmetry factor at NMC") + return Positive_electrode_exchange_current_density_A_per_m2(cₚ / cₚ_max) + + +def nmc_exchange_prefactor_derivative(cₑ, cₚ, cₚ_max, T): + αₙₚ = Parameter("Anodic symmetry factor at NMC") + αₚₚ = Parameter("Cathodic symmetry factor at NMC") + return Positive_electrode_exchange_current_density_A_per_m2(cₚ / cₚ_max) * αₚₚ / cₑ + + +def Negative_electrode_exchange_current_density_A_per_m2(SOC): + return 10 * exp( + (SOC < 0.025) + * ((27.834463836969007 * SOC + 15.96473196936472) * SOC + -2.37212952830029) + + (SOC >= 0.025) + * (SOC < 0.035) + * ((-421.1589223786759 * SOC + 38.41440128014693) * SOC + -2.6527503946850706) + + (SOC >= 0.035) + * (SOC < 0.045) + * ((217.6317365328614 * SOC + -6.300944843660545) * SOC + -1.8702318375184337) + + (SOC >= 0.045) + * (SOC < 0.055) + * ((-164.31938405500932 * SOC + 28.074656009247747) * SOC + -2.6436828567088675) + + (SOC >= 0.055) + * (SOC < 0.065) + * ((12.249126997133317 * SOC + 8.652119793512156) * SOC + -2.1095631107761363) + + (SOC >= 0.065) + * (SOC < 0.07500000000000001) + * ((-265.95640077690587 * SOC + 44.81883840413718) * SOC + -3.2849814656214527) + + (SOC >= 0.07500000000000001) + * (SOC < 0.08499999999999999) + * ((-68.90850892584058 * SOC + 15.261654626477366) * SOC + -2.1765870739592046) + + (SOC >= 0.08499999999999999) + * (SOC < 0.095) + * ((62.94976632331509 * SOC + -7.154252165879257) * SOC + -1.2239110352840648) + + (SOC >= 0.095) + * (SOC < 0.125) + * ((-55.70510160059371 * SOC + 15.39017273966347) * SOC + -2.294771218297342) + + (SOC >= 0.125) + * (SOC < 0.175) + * ((-8.718685210812424 * SOC + 3.6435686422181988) * SOC + -1.5606084622070107) + + (SOC >= 0.175) + * (SOC < 0.225) + * ((-1.3638025119056465 * SOC + 1.0693596976008024) * SOC + -1.3353651795529888) + + (SOC >= 0.225) + * (SOC < 0.275) + * ((-1.664089416211823 * SOC + 1.2044888045385846) * SOC + -1.3505672040834895) + + (SOC >= 0.275) + * (SOC < 0.325) + * ((-5.8675333935170215 * SOC + 3.5163829920564496) * SOC + -1.6684526548671954) + + (SOC >= 0.325) + * (SOC < 0.375) + * ((8.136678418185596 * SOC + -5.586354685550219) * SOC + -0.18925778225610657) + + (SOC >= 0.375) + * (SOC < 0.42500000000000004) + * ((-12.650013231256338 * SOC + 10.003664051531217) * SOC + -3.112386295458876) + + (SOC >= 0.42500000000000004) + * (SOC < 0.475) + * ((7.051634371483175 * SOC + -6.742736410797363) * SOC + 0.4462238027859371) + + (SOC >= 0.475) + * (SOC < 0.525) + * ((13.008911741853609 * SOC + -12.402149912649264) * SOC + 1.7903345094757839) + + (SOC >= 0.525) + * (SOC < 0.575) + * ((-5.123050444802772 * SOC + 6.636410383339978) * SOC + -3.2072875682214033) + + (SOC >= 0.575) + * (SOC < 0.625) + * ((-11.687431501401761 * SOC + 14.185448598428735) * SOC + -5.377636055059426) + + (SOC >= 0.625) + * (SOC < 0.675) + * ((27.71556168605477 * SOC + -35.06829288589191) * SOC + 10.014158158790764) + + (SOC >= 0.675) + * (SOC < 0.7250000000000001) + * ((-26.145992534672075 * SOC + 37.64480531208926) * SOC + -14.526512483027872) + + (SOC >= 0.7250000000000001) + * (SOC < 0.775) + * ((4.67282667236347 * SOC + -7.042482538112154) * SOC + 1.6726293626701505) + + (SOC >= 0.775) + * (SOC < 0.825) + * ((4.586035801125973 * SOC + -6.907956687694082) * SOC + 1.6205005956331178) + + (SOC >= 0.825) + * (SOC < 0.875) + * ((-7.501318516918673 * SOC + 13.036177937079515) * SOC + -6.606454937085999) + + (SOC >= 0.875) + * (SOC < 0.905) + * ((34.06912325716132 * SOC + -59.71209516756062) * SOC + 25.220914546193967) + + (SOC >= 0.905) + * (SOC < 0.915) + * ((-198.116031360747 * SOC + 360.54303469085517) * SOC + -164.9445317147397) + + (SOC >= 0.915) + * (SOC < 0.925) + * ((155.68989428587975 * SOC + -286.92180924247623) * SOC + 131.2706343847599) + + (SOC >= 0.925) + * (SOC < 0.935) + * ((-44.90589213903786 * SOC + 84.18039564362516) * SOC + -40.36413537506178) + + (SOC >= 0.935) + * (SOC < 0.9450000000000001) + * ((-17.448644009964482 * SOC + 32.835341642256026) * SOC + -16.360322629421717) + + (SOC >= 0.9450000000000001) + * (SOC < 0.9550000000000001) + * ((75.85275901431578 * SOC + -143.50431007363295) * SOC + 66.96016280633694) + + (SOC >= 0.9550000000000001) + * (SOC < 0.965) + * ((-214.9570579656629 * SOC + 411.94244035812335) * SOC + -198.2656605248285) + + (SOC >= 0.965) + * (SOC < 0.975) + * ((136.24296138879663 * SOC + -265.87359699598073) * SOC + 128.7805774985277) + + (SOC >= 0.975) + * ((-92.57872949167268 * SOC + 180.32870022093448) * SOC + -88.74304239471894) + ) + + +def derivative_Negative_electrode_exchange_current_density_A_per_m2(SOC): + return ( + 10 + * Negative_electrode_exchange_current_density_A_per_m2(SOC) + * ( + (SOC < 0.025) * (55.668927673938015 * SOC + 15.96473196936472) + + (SOC >= 0.025) + * (SOC < 0.035) + * (-842.3178447573518 * SOC + 38.41440128014693) + + (SOC >= 0.035) + * (SOC < 0.045) + * (435.2634730657228 * SOC + -6.300944843660545) + + (SOC >= 0.045) + * (SOC < 0.055) + * (-328.63876811001865 * SOC + 28.074656009247747) + + (SOC >= 0.055) + * (SOC < 0.065) + * (24.498253994266634 * SOC + 8.652119793512156) + + (SOC >= 0.065) + * (SOC < 0.07500000000000001) + * (-531.9128015538117 * SOC + 44.81883840413718) + + (SOC >= 0.07500000000000001) + * (SOC < 0.08499999999999999) + * (-137.81701785168116 * SOC + 15.261654626477366) + + (SOC >= 0.08499999999999999) + * (SOC < 0.095) + * (125.89953264663018 * SOC + -7.154252165879257) + + (SOC >= 0.095) + * (SOC < 0.125) + * (-111.41020320118741 * SOC + 15.39017273966347) + + (SOC >= 0.125) + * (SOC < 0.175) + * (-17.43737042162485 * SOC + 3.6435686422181988) + + (SOC >= 0.175) + * (SOC < 0.225) + * (-2.727605023811293 * SOC + 1.0693596976008024) + + (SOC >= 0.225) + * (SOC < 0.275) + * (-3.328178832423646 * SOC + 1.2044888045385846) + + (SOC >= 0.275) + * (SOC < 0.325) + * (-11.735066787034043 * SOC + 3.5163829920564496) + + (SOC >= 0.325) + * (SOC < 0.375) + * (16.273356836371192 * SOC + -5.586354685550219) + + (SOC >= 0.375) + * (SOC < 0.42500000000000004) + * (-25.300026462512676 * SOC + 10.003664051531217) + + (SOC >= 0.42500000000000004) + * (SOC < 0.475) + * (14.10326874296635 * SOC + -6.742736410797363) + + (SOC >= 0.475) + * (SOC < 0.525) + * (26.017823483707218 * SOC + -12.402149912649264) + + (SOC >= 0.525) + * (SOC < 0.575) + * (-10.246100889605543 * SOC + 6.636410383339978) + + (SOC >= 0.575) + * (SOC < 0.625) + * (-23.374863002803522 * SOC + 14.185448598428735) + + (SOC >= 0.625) + * (SOC < 0.675) + * (55.43112337210954 * SOC + -35.06829288589191) + + (SOC >= 0.675) + * (SOC < 0.7250000000000001) + * (-52.29198506934415 * SOC + 37.64480531208926) + + (SOC >= 0.7250000000000001) + * (SOC < 0.775) + * (9.34565334472694 * SOC + -7.042482538112154) + + (SOC >= 0.775) + * (SOC < 0.825) + * (9.172071602251947 * SOC + -6.907956687694082) + + (SOC >= 0.825) + * (SOC < 0.875) + * (-15.002637033837345 * SOC + 13.036177937079515) + + (SOC >= 0.875) + * (SOC < 0.905) + * (68.13824651432265 * SOC + -59.71209516756062) + + (SOC >= 0.905) + * (SOC < 0.915) + * (-396.232062721494 * SOC + 360.54303469085517) + + (SOC >= 0.915) + * (SOC < 0.925) + * (311.3797885717595 * SOC + -286.92180924247623) + + (SOC >= 0.925) + * (SOC < 0.935) + * (-89.81178427807572 * SOC + 84.18039564362516) + + (SOC >= 0.935) + * (SOC < 0.9450000000000001) + * (-34.897288019928965 * SOC + 32.835341642256026) + + (SOC >= 0.9450000000000001) + * (SOC < 0.9550000000000001) + * (151.70551802863156 * SOC + -143.50431007363295) + + (SOC >= 0.9550000000000001) + * (SOC < 0.965) + * (-429.9141159313258 * SOC + 411.94244035812335) + + (SOC >= 0.965) + * (SOC < 0.975) + * (272.48592277759326 * SOC + -265.87359699598073) + + (SOC >= 0.975) * (-185.15745898334535 * SOC + 180.32870022093448) + ) + ) + + +def Positive_electrode_exchange_current_density_A_per_m2(SOC): + return 10 * exp( + (SOC < 0.025) + * ((54.0509009047164 * SOC + 3.998665364455661) * SOC + -3.691443849291202) + + (SOC >= 0.025) + * (SOC < 0.035) + * ((273.8859954657528 * SOC + -6.993089363596141) * SOC + -3.554046915190554) + + (SOC >= 0.035) + * (SOC < 0.045) + * ((-143.05889379391738 * SOC + 22.19305288458122) * SOC + -4.06480440453365) + + (SOC >= 0.045) + * (SOC < 0.055) + * ((130.55575990440957 * SOC + -2.432265948268423) * SOC + -3.5107347307945282) + + (SOC >= 0.055) + * (SOC < 0.065) + * ((-754.7160368399291 * SOC + 94.94763169360886) * SOC + -6.1886819159461695) + + (SOC >= 0.065) + * (SOC < 0.07500000000000001) + * ((334.39532447137026 * SOC + -46.63684527685973) * SOC + -1.5871864144059202) + + (SOC >= 0.07500000000000001) + * (SOC < 0.08499999999999999) + * ((-229.28479085109575 * SOC + 37.91517202150908) * SOC + -4.757887063094785) + + (SOC >= 0.08499999999999999) + * (SOC < 0.095) + * ((874.1806496138752 * SOC + -149.67395285753605) * SOC + 3.214650744264631) + + (SOC >= 0.095) + * (SOC < 0.125) + * ((-353.73204821775016 * SOC + 83.62945973047312) * SOC + -7.867261353665796) + + (SOC >= 0.125) + * (SOC < 0.175) + * ((84.72715236061845 * SOC + -25.985340414118866) * SOC + -1.0163363446287885) + + (SOC >= 0.175) + * (SOC < 0.225) + * ((-37.203377576212915 * SOC + 16.690345063772128) * SOC + -4.750458823944243) + + (SOC >= 0.225) + * (SOC < 0.275) + * ((17.518468468908395 * SOC + -7.93448565653253) * SOC + -1.9801653679099687) + + (SOC >= 0.275) + * (SOC < 0.325) + * ((-17.48343037918687 * SOC + 11.316558709919832) * SOC + -4.627183968297189) + + (SOC >= 0.325) + * (SOC < 0.375) + * ((36.811763420237526 * SOC + -23.975317259706003) * SOC + 1.1077458767670336) + + (SOC >= 0.375) + * (SOC < 0.42500000000000004) + * ((-28.853775299433096 * SOC + 25.273836780046963) * SOC + -8.126470505686669) + + (SOC >= 0.42500000000000004) + * (SOC < 0.475) + * ((-15.027367958237164 * SOC + 13.521390540030438) * SOC + -5.629075679683169) + + (SOC >= 0.475) + * (SOC < 0.525) + * ((14.027957971693013 * SOC + -14.081169093403162) * SOC + 0.9265322332573334) + + (SOC >= 0.525) + * (SOC < 0.575) + * ((3.033601250402455 * SOC + -2.53709453604813) * SOC + -2.1037873380483063) + + (SOC >= 0.575) + * (SOC < 0.625) + * ((-5.979243172813028 * SOC + 7.827676550649585) * SOC + -5.083659025473992) + + (SOC >= 0.625) + * (SOC < 0.675) + * ((-0.8724007621615328 * SOC + 1.4441235373352583) * SOC + -3.088798708813215) + + (SOC >= 0.675) + * (SOC < 0.7250000000000001) + * ((-17.167772160183915 * SOC + 23.44287492466549) * SOC + -10.513377302037185) + + (SOC >= 0.7250000000000001) + * (SOC < 0.775) + * ((-17.31420489184609 * SOC + 23.655202385575876) * SOC + -10.590346006617153) + + (SOC >= 0.775) + * (SOC < 0.825) + * ((-3.67823352233529 * SOC + 2.5194467628340362) * SOC + -2.4002407028046946) + + (SOC >= 0.825) + * (SOC < 0.875) + * ((13.257260038885875 * SOC + -25.424117613180897) * SOC + 9.126479602301288) + + (SOC >= 0.875) + * (SOC < 0.905) + * ((-290.7726722841667 * SOC + 506.6282639521596) * SOC + -223.6464373325357) + + (SOC >= 0.905) + * (SOC < 0.915) + * ((1179.0782200890171 * SOC + -2153.8018512433046) * SOC + 980.1981897934111) + + (SOC >= 0.915) + * (SOC < 0.925) + * ((-917.3897065212404 * SOC + 1682.7344544534644) * SOC + -775.0171700628598) + + (SOC >= 0.925) + * (SOC < 0.935) + * ((867.0386453967003 * SOC + -1618.4579965947196) * SOC + 751.7843385469278) + + (SOC >= 0.935) + * (SOC < 0.9450000000000001) + * ((-479.85724574424603 * SOC + 900.237319838845) * SOC + -425.70572188576443) + + (SOC >= 0.9450000000000001) + * (SOC < 0.9550000000000001) + * ((-62.156860698272794 * SOC + 110.78359210195049) * SOC + -52.68883553008527) + + (SOC >= 0.9550000000000001) + * (SOC < 0.965) + * ((-1437.1297827765193 * SOC + 2736.981873271412) * SOC + -1306.6985147884952) + + (SOC >= 0.965) + * (SOC < 0.975) + * ((2251.176102712596 * SOC + -4381.448485722583) * SOC + 2127.944133426101) + + (SOC >= 0.975) + * ((-1191.2980844287886 * SOC + 2331.376179203122) * SOC + -1144.5578907251765) + ) + + +def derivative_Positive_electrode_exchange_current_density_A_per_m2(SOC): + return ( + 10 + * Positive_electrode_exchange_current_density_A_per_m2(SOC) + * ( + (SOC < 0.025) * (108.1018018094328 * SOC + 3.998665364455661) + + (SOC >= 0.025) + * (SOC < 0.035) + * (547.7719909315056 * SOC + -6.993089363596141) + + (SOC >= 0.035) + * (SOC < 0.045) + * (-286.11778758783475 * SOC + 22.19305288458122) + + (SOC >= 0.045) + * (SOC < 0.055) + * (261.11151980881914 * SOC + -2.432265948268423) + + (SOC >= 0.055) + * (SOC < 0.065) + * (-1509.4320736798581 * SOC + 94.94763169360886) + + (SOC >= 0.065) + * (SOC < 0.07500000000000001) + * (668.7906489427405 * SOC + -46.63684527685973) + + (SOC >= 0.07500000000000001) + * (SOC < 0.08499999999999999) + * (-458.5695817021915 * SOC + 37.91517202150908) + + (SOC >= 0.08499999999999999) + * (SOC < 0.095) + * (1748.3612992277504 * SOC + -149.67395285753605) + + (SOC >= 0.095) + * (SOC < 0.125) + * (-707.4640964355003 * SOC + 83.62945973047312) + + (SOC >= 0.125) + * (SOC < 0.175) + * (169.4543047212369 * SOC + -25.985340414118866) + + (SOC >= 0.175) + * (SOC < 0.225) + * (-74.40675515242583 * SOC + 16.690345063772128) + + (SOC >= 0.225) + * (SOC < 0.275) + * (35.03693693781679 * SOC + -7.93448565653253) + + (SOC >= 0.275) + * (SOC < 0.325) + * (-34.96686075837374 * SOC + 11.316558709919832) + + (SOC >= 0.325) + * (SOC < 0.375) + * (73.62352684047505 * SOC + -23.975317259706003) + + (SOC >= 0.375) + * (SOC < 0.42500000000000004) + * (-57.70755059886619 * SOC + 25.273836780046963) + + (SOC >= 0.42500000000000004) + * (SOC < 0.475) + * (-30.05473591647433 * SOC + 13.521390540030438) + + (SOC >= 0.475) + * (SOC < 0.525) + * (28.055915943386026 * SOC + -14.081169093403162) + + (SOC >= 0.525) + * (SOC < 0.575) + * (6.06720250080491 * SOC + -2.53709453604813) + + (SOC >= 0.575) + * (SOC < 0.625) + * (-11.958486345626056 * SOC + 7.827676550649585) + + (SOC >= 0.625) + * (SOC < 0.675) + * (-1.7448015243230657 * SOC + 1.4441235373352583) + + (SOC >= 0.675) + * (SOC < 0.7250000000000001) + * (-34.33554432036783 * SOC + 23.44287492466549) + + (SOC >= 0.7250000000000001) + * (SOC < 0.775) + * (-34.62840978369218 * SOC + 23.655202385575876) + + (SOC >= 0.775) + * (SOC < 0.825) + * (-7.35646704467058 * SOC + 2.5194467628340362) + + (SOC >= 0.825) + * (SOC < 0.875) + * (26.51452007777175 * SOC + -25.424117613180897) + + (SOC >= 0.875) + * (SOC < 0.905) + * (-581.5453445683333 * SOC + 506.6282639521596) + + (SOC >= 0.905) + * (SOC < 0.915) + * (2358.1564401780342 * SOC + -2153.8018512433046) + + (SOC >= 0.915) + * (SOC < 0.925) + * (-1834.7794130424809 * SOC + 1682.7344544534644) + + (SOC >= 0.925) + * (SOC < 0.935) + * (1734.0772907934006 * SOC + -1618.4579965947196) + + (SOC >= 0.935) + * (SOC < 0.9450000000000001) + * (-959.7144914884921 * SOC + 900.237319838845) + + (SOC >= 0.9450000000000001) + * (SOC < 0.9550000000000001) + * (-124.31372139654559 * SOC + 110.78359210195049) + + (SOC >= 0.9550000000000001) + * (SOC < 0.965) + * (-2874.2595655530386 * SOC + 2736.981873271412) + + (SOC >= 0.965) + * (SOC < 0.975) + * (4502.352205425192 * SOC + -4381.448485722583) + + (SOC >= 0.975) * (-2382.596168857577 * SOC + 2331.376179203122) + ) + ) + + +################## +# Diffusivities. # +################## + + +# Discharge profile - Delithiation direction of the Anode +def Negative_electrode_diffusivity__m2_s_1_(Negative_electrode_SOC____, T): + return exp( + (Negative_electrode_SOC____ < 0.09913139478844102) + * ( + (28.968262165814508 * Negative_electrode_SOC____ + 5.625713621193107) + * Negative_electrode_SOC____ + + -33.24419027877485 + ) + + (Negative_electrode_SOC____ >= 0.09913139478844102) + * (Negative_electrode_SOC____ < 0.4956458667127428) + * ( + (-24.866807392848557 * Negative_electrode_SOC____ + 16.29920468895898) + * Negative_electrode_SOC____ + + -33.77322930717966 + ) + + (Negative_electrode_SOC____ >= 0.4956458667127428) + * (Negative_electrode_SOC____ < 0.6939030739440655) + * ( + (-41.70002742886413 * Negative_electrode_SOC____ + 32.985836557593416) + * Negative_electrode_SOC____ + + -37.90855936470258 + ) + + (Negative_electrode_SOC____ >= 0.6939030739440655) + * (Negative_electrode_SOC____ < 0.8921603579383661) + * ( + (183.2186924065055 * Negative_electrode_SOC____ + -279.15774560506054) + * Negative_electrode_SOC____ + + 70.3901362225862 + ) + + (Negative_electrode_SOC____ >= 0.8921603579383661) + * ( + (-301.42626762424334 * Negative_electrode_SOC____ + 585.6042964230546) + * Negative_electrode_SOC____ + + -315.36307025107135 + ) + ) + + +def derivative_Negative_electrode_diffusivity__m2_s_1_(Negative_electrode_SOC____, T): + return Negative_electrode_diffusivity__m2_s_1_(Negative_electrode_SOC____) * ( + (Negative_electrode_SOC____ < 0.09913139478844102) + * (57.936524331629016 * Negative_electrode_SOC____ + 5.625713621193107) + + (Negative_electrode_SOC____ >= 0.09913139478844102) + * (Negative_electrode_SOC____ < 0.4956458667127428) + * (-49.733614785697114 * Negative_electrode_SOC____ + 16.29920468895898) + + (Negative_electrode_SOC____ >= 0.4956458667127428) + * (Negative_electrode_SOC____ < 0.6939030739440655) + * (-83.40005485772826 * Negative_electrode_SOC____ + 32.985836557593416) + + (Negative_electrode_SOC____ >= 0.6939030739440655) + * (Negative_electrode_SOC____ < 0.8921603579383661) + * (366.437384813011 * Negative_electrode_SOC____ + -279.15774560506054) + + (Negative_electrode_SOC____ >= 0.8921603579383661) + * (-602.8525352484867 * Negative_electrode_SOC____ + 585.6042964230546) + ) + + +# Discharge profile - lithiation direction of cathode + + +def Positive_electrode_diffusivity__m2_s_1_(Positive_electrode_SOC____, T): + return exp( + (Positive_electrode_SOC____ < 0.3382277436623123) + * ( + (27.982955482386558 * Positive_electrode_SOC____ + -12.143707202132319) + * Positive_electrode_SOC____ + + -33.0217804711556 + ) + + (Positive_electrode_SOC____ >= 0.3382277436623123) + * (Positive_electrode_SOC____ < 0.6665106367866339) + * ( + (-11.319201599513178 * Positive_electrode_SOC____ + 14.44245261961305) + * Positive_electrode_SOC____ + + -37.51786889573288 + ) + + (Positive_electrode_SOC____ >= 0.6665106367866339) + * (Positive_electrode_SOC____ < 0.8822401369572627) + * ( + (-25.591142520785183 * Positive_electrode_SOC____ + 33.4672534828494) + * Positive_electrode_SOC____ + + -43.85798496478026 + ) + + (Positive_electrode_SOC____ >= 0.8822401369572627) + * ( + (-190.3629702643384 * Positive_electrode_SOC____ + 324.2038932331925) + * Positive_electrode_SOC____ + + -172.10775140069836 + ) + ) + + +def derivative_Positive_electrode_diffusivity__m2_s_1_(Positive_electrode_SOC____, T): + return Positive_electrode_diffusivity__m2_s_1_(Positive_electrode_SOC____) * ( + (Positive_electrode_SOC____ < 0.3382277436623123) + * (55.965910964773116 * Positive_electrode_SOC____ + -12.143707202132319) + + (Positive_electrode_SOC____ >= 0.3382277436623123) + * (Positive_electrode_SOC____ < 0.6665106367866339) + * (-22.638403199026357 * Positive_electrode_SOC____ + 14.44245261961305) + + (Positive_electrode_SOC____ >= 0.6665106367866339) + * (Positive_electrode_SOC____ < 0.8822401369572627) + * (-51.18228504157037 * Positive_electrode_SOC____ + 33.4672534828494) + + (Positive_electrode_SOC____ >= 0.8822401369572627) + * (-380.7259405286768 * Positive_electrode_SOC____ + 324.2038932331925) + ) + + +parameters = SubstitutionDict( + { + ####################################################################### + # Boilerplate values that PyBaMM requires but are not important here. # + ####################################################################### + "Current function [A]": 0, + "Positive electrode cation signed stoichiometry": -1.0, + "Positive electrode OCP entropic change [V.K-1]": 0, + "Positive electrode OCP entropic change partial derivative by SOC [V.K-1]": 0, + "Negative electrode cation signed stoichiometry": -1.0, + "Negative electrode OCP entropic change [V.K-1]": 0, + "Negative electrode OCP entropic change partial derivative by SOC [V.K-1]": 0, + "Number of electrodes connected in parallel to make a cell": 1, + "Number of cells connected in series to make a battery": 1, + ################# + # Placeholders. # + ################# + "Negative electrode double-layer capacity [F.m-2]": 2e-3, + "Positive electrode double-layer capacity [F.m-2]": 2e-3, + # Based on EC:EMC 3:7 weight-% from Landesfeind2019. + "Thermodynamic factor": 1.7, + ######################### + # Geometric properties. # + ######################### + # Circular cut-out with diameter 18 mm for the measurements. + "Current collector perpendicular area [m2]": 2.5447e-4, + ################## + # Cell capacity. # + ################## + # Currently refers to the anode, even though it is oversized. + # Note that it is directly used to calculate the maximum + # lithiation of the negative electrode. + "Nominal cell capacity [A.h]": 0.01202, + "Negative electrode capacity [A.h]": 0.01202, + ################ + # Temperature. # + ################ + "Reference temperature [K]": 298.46, + "Ambient temperature [K]": 298.46, + "Initial temperature [K]": 298.46, + ############################## + # Reaction symmetry factors. # + ############################## + "Anodic symmetry factor at Graphite-SiOx": αₙₙ, + "Cathodic symmetry factor at Graphite-SiOx": αₚₙ, + "Anodic symmetry factor at NMC": αₙₚ, + "Cathodic symmetry factor at NMC": αₚₚ, + ############################################################## + # Commonly used expressions for the Butler-Volmer prefactor. # + ############################################################## + "Negative electrode exchange-current density [A.m-2]": ( + 4e-2 # graphite_exchange_prefactor + ), + "Positive electrode exchange-current density [A.m-2]": (nmc_exchange_prefactor), + ( + "Negative electrode exchange-current density partial derivative " + "by electrolyte concentration [A.m.mol-1]" + ): graphite_exchange_prefactor_derivative, + ( + "Positive electrode exchange-current density partial derivative " + "by electrolyte concentration [A.m.mol-1]" + ): nmc_exchange_prefactor_derivative, + ####################################### + # Graphite-SiOx electrode parameters. # + ####################################### + "Negative electrode thickness [m]": 88e-6, + # Mean value of the range of possible porosity values. + "Negative electrode porosity": 0.26, + # Specific surface area from surface area (98 cm²) and volume (0.022 cm³) + # measurements, then matched to a = 3 * (1 - ε) / R, would give 4.983e-6. + "Negative particle radius [m]": 4.785e-6, + "Negative electrode electrons in reaction": 1.0, + "Negative electrode OCP [V]": ( + _18650_LG_3500_MJ1_Anode_GITT_dch_ohne_EIS_2320__extraction + ), + "Negative electrode OCP derivative by SOC [V]": ( + derivative__18650_LG_3500_MJ1_Anode_GITT_dch_ohne_EIS_2320__extraction + ), + # Lab report gave 4.6 +- 0.9. + "Negative electrode tortuosity square": 1.55, + "Negative electrode conductivity [S.m-1]": 1.4e-4, # 215, adjusted to match highest frequency semicircle + "Negative particle diffusivity [m2.s-1]": Negative_electrode_diffusivity__m2_s_1_, + #################################### + # NMC-830512 electrode properties. # + #################################### + "Positive electrode thickness [m]": 74e-6, + # Mean value of the range of possible porosity values. + "Positive electrode porosity": 0.23, + "Positive particle radius [m]": 2.375e-6, + "Positive electrode electrons in reaction": 1.0, + # With no half-cell data yet, use this literature value for NMC-811 for the + # parameterization of the negative electrode. It should have only a minor + # influence, but grossly wrong values might introduce unnecessary noise. + "Positive electrode OCP [V]": ( + _18650_LG_3500_MJ1_Cathode_GITT_ch_ohne_EIS_2327__extraction + ), + "Positive electrode OCP derivative by SOC [V]": ( + derivative__18650_LG_3500_MJ1_Cathode_GITT_ch_ohne_EIS_2327__extraction + ), + # Lab report gave 2.5 +- 0.1. + "Positive electrode tortuosity square": 2.5, + "Positive electrode conductivity [S.m-1]": 0.25, + "Positive particle diffusivity [m2.s-1]": Positive_electrode_diffusivity__m2_s_1_, + ################################################################# + # Separator (replaced by Whatman GF/A from El-Cell) parameters. # + ################################################################# + "Separator thickness [m]": 260e-6, + "Separator porosity": 0.93, + "Separator Bruggeman coefficient (electrolyte)": 1.0, + "Separator Bruggeman coefficient (electrode)": 1.0, + ################################################################################## + # Electrolyte parameters. 1 M LiPF 6 in EC:EMC:DMC (1:1:1 vol%) with 2 % wt. VC. # + ################################################################################## + "Typical electrolyte concentration [mol.m-3]": 1000.0, + "Initial concentration in electrolyte [mol.m-3]": 1000.0, + "Cation transference number": 0.226, + "Electrolyte diffusivity [m2.s-1]": 2.5e-10, + "Electrolyte conductivity [S.m-1]": 5.5e-3, # 1.026, adjusted to match offset + ################### + # Voltage window. # + ################### + "Lower voltage cut-off [V]": 2.5, + "Upper voltage cut-off [V]": 4.2, + #################################################################### + # SEI properties, mostly taken from Single2019 and Kolzenberg2020. # + #################################################################### + "Charge number of anion in electrolyte salt dissociation": -1, + "Charge number of cation in electrolyte salt dissociation": 1, + "Stoichiometry of anion in electrolyte salt dissociation": 1, + "Stoichiometry of cation in electrolyte salt dissociation": 1, + "Mass density of electrolyte [kg.m-3]": 1.5e3, # Wikipedia, 20 °C + "Mass density of cations in electrolyte [kg.m-3]": 68.53, # Li+PF6- + "Molar mass of electrolyte solvent [kg.mol-1]": 89.08e-3, + "Solvent concentration [mol.m-3]": 13.17e3, + "Partial molar volume of electrolyte solvent [m3.mol-1]": 75.93e-6, + "Inital SEI thickness [m]": 2e-9, + "Reference apparent SEI thickness [m]": 0.05e-9, + "SEI tunneling length [m]": 2.05e-9, + "SEI formation rate constant [A.m-2]": 7.04e-5, + "SEI formation symmetry factor": 0.22, + "SEI molar volume [m3.mol-1]": 1.078e-5, + "SEI thickness [m]": 2e-7, + "SEI ionic conductivity [S.m-1]": 2e-2, + "SEI relative permittivity": 131, + "SEI lithium reference potential [J.mol-1]": 17400, + "Anion transference number in SEI": 1 - 0.063, # t_plus is 0.063 + "SEI porosity": 0.08, + "SEI Bruggeman coefficient": 4.54, + }, + substitutions={ + ########################################################### + # Scalings that are wholly dependent on other parameters. # + ########################################################### + "Typical current [A]": "Nominal cell capacity [A.h]", + "Electrode height [m]": 0.06, + "Electrode width [m]": 1.23494869661675, + "Cell volume [m3]": cell_volume, + ############################################################### + # Given the cell capacity Q, we adjust the maximum lithiation # + # of the electrode via Q = (1 - ε) * L * A * c_max * F. # + ############################################################### + # From the lab report: 1460 mol/m³ Li in Graphite-SiOx at SOC 0, + # 32452 mol/m³ Li in Graphite-SiOx at SOC 100. + # Composition is 96.5% graphite and 3.5% SiOx. With the assumption + # that the SiOx does not contribute to the dynamics at the measured + # SOCs, we assign the values to the graphite. Literature states a + # similar value for it anyways: Ecker2015, Schmalstieg2018 (lliondb.com). + "Maximum concentration in negative electrode [mol.m-3]": ( + negative_concentration_correction + ), + # The provided values give roughly 20000 mol/m³. + "Maximum concentration in positive electrode [mol.m-3]": ( + positive_concentration + ), + # These are just placeholders until the parameterization routine + # changes them. + "Initial concentration in negative electrode [mol.m-3]": ( + placeholder_initial_negative_concentration + ), + "Initial concentration in positive electrode [mol.m-3]": ( + placeholder_initial_positive_concentration + ), + ######################################################### + # Active material as volume fraction of all components. # + ######################################################### + "Negative electrode active material volume fraction": ( + negative_volume_fraction + ), + "Positive electrode active material volume fraction": ( + positive_volume_fraction + ), + ################################################################# + # Conversion from tortuosity τ in the lab report to Bruggeman. # + # The lab report uses the notation N_M = ε / τ # + # instead of N_M = ε / τ² for the MacMullin number. # + # So we first take the square root of the value from there. # + # Formula with τ as the square root: β = 1 - 2 log(τ) / log(ε). # + ################################################################# + "Negative electrode Bruggeman coefficient (electrolyte)": (negative_bruggeman), + "Negative electrode Bruggeman coefficient (electrode)": (negative_bruggeman), + "Positive electrode Bruggeman coefficient (electrolyte)": (positive_bruggeman), + "Positive electrode Bruggeman coefficient (electrode)": (positive_bruggeman), + ##################################################################### + # Match the particle radius as it is used in the derivation of the # + # model equations: just an effective parameter that scales with the # + # active surface area to volume ratio a = 3 * (1 - ε) / R. # + # Note: do not calculate the particle radius here instead, as the # + # solver setup requires it to be known when building the model. # + ##################################################################### + # The provided values give around 226630 1/m. + "Negative electrode surface area to volume ratio [m-1]": ( + negative_effective_surface_area + ), + # The provided values give around 460737 1/m. + "Positive electrode surface area to volume ratio [m-1]": ( + positive_effective_surface_area + ), + ################################################### + # PyBaMM names for the reaction symmetry factors. # + ################################################### + "Negative electrode anodic charge-transfer coefficient": ( + "Anodic symmetry factor at Graphite-SiOx" + ), + "Negative electrode cathodic charge-transfer coefficient": ( + "Cathodic symmetry factor at Graphite-SiOx" + ), + "Positive electrode anodic charge-transfer coefficient": ( + "Anodic symmetry factor at NMC" + ), + "Positive electrode cathodic charge-transfer coefficient": ( + "Cathodic symmetry factor at NMC" + ), + "Open-circuit voltage at 100% SOC [V]": "Upper voltage cut-off [V]", + "Open-circuit voltage at 0% SOC [V]": "Lower voltage cut-off [V]", + }, +) + + +def get_parameters_with_switched_electrodes(): + original_parameters = deepcopy(parameters) + switched_parameters = deepcopy(parameters) + # Switch the parameters between negative and positive electrode. + for k, v in original_parameters.items(): + if "negative" in k: + switched_parameters[k.replace("negative", "positive")] = v + elif "positive" in k: + switched_parameters[k.replace("positive", "negative")] = v + elif "Negative" in k: + switched_parameters[k.replace("Negative", "Positive")] = v + elif "Positive" in k: + switched_parameters[k.replace("Positive", "Negative")] = v + return switched_parameters diff --git a/examples/data/Wycisk2024/README.md b/examples/data/Wycisk2024/README.md new file mode 100644 index 000000000..50d4a6df9 --- /dev/null +++ b/examples/data/Wycisk2024/README.md @@ -0,0 +1,11 @@ +# Galvanostatic Intermittent Titration Technique on graphite with 5% silicon + +This repository contains a subset of the data published in the research paper: +**D. Wycisk et al., "Challenges of open-circuit voltage measurements for silicon-containing Li-Ion cells,"** *Journal of Energy Storage*, vol. 89, p. 111617, June 2024. [DOI: 10.1016/j.est.2024.111617](https://doi.org/10.1016/j.est.2024.111617) + +## Citation + +If you use this dataset in your work, please cite the original publication as follows: + +```plaintext +D. Wycisk et al., "Challenges of open-circuit voltage measurements for silicon-containing Li-Ion cells," Journal of Energy Storage, vol. 89, p. 111617, June 2024, doi: 10.1016/j.est.2024.111617. diff --git a/examples/data/Wycisk2024/gitt_on_graphite_with_5_percent_silicon.parquet b/examples/data/Wycisk2024/gitt_on_graphite_with_5_percent_silicon.parquet new file mode 100644 index 000000000..03999cdb1 Binary files /dev/null and b/examples/data/Wycisk2024/gitt_on_graphite_with_5_percent_silicon.parquet differ diff --git a/examples/data/Wycisk2024/read_dataset.py b/examples/data/Wycisk2024/read_dataset.py new file mode 100644 index 000000000..368094c19 --- /dev/null +++ b/examples/data/Wycisk2024/read_dataset.py @@ -0,0 +1,53 @@ +from pyarrow import parquet +from pybamm import citations + +from pybop import Datasets + + +def read_parquet_cycling(filename): + datafile = parquet.ParquetFile(filename) + indices = [] + timepoints = [] + currents = [] + voltages = [] + for row_group_number in range(datafile.num_row_groups): + row_group = datafile.read_row_group(row_group_number) + indices.append(row_group.column("indices")[0].as_py()) + timepoints.append( + row_group.column("timepoints [s]").combine_chunks().to_numpy().tolist() + ) + currents.append( + row_group.column("currents [A]").combine_chunks().to_numpy().tolist() + ) + voltages.append( + row_group.column("voltages [V]").combine_chunks().to_numpy().tolist() + ) + return Datasets( + [ + { + "Time [s]": t, + "Current function [A]": c, + "Voltage change [V]": v, + "Cycle indices": [i] * len(t), + } + for t, c, v, i in zip(timepoints, currents, voltages, indices, strict=False) + ], + domain="Time [s]", + control_variable="Current function [A]", + ) + + +citations.register("""@article{ + Wycisk2024, + title={{Challenges of open-circuit voltage measurements for silicon-containing Li-Ion cells}}, + author={Wycisk, D and Mertin, G and Oldenburger, M and von Kessel, O and Latz, A}, + journal={Journal of Energy Storage}, + volume={89}, + pages={111617}, + year={2024}, + doi={10.1016/j.est.2024.111617} +}""") + +gitt_on_graphite_with_5_percent_silicon = read_parquet_cycling( + "../../data/Wycisk2024/gitt_on_graphite_with_5_percent_silicon.parquet" +) diff --git a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py index d851937d4..58fa4ba4f 100644 --- a/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py +++ b/examples/scripts/battery_parameterisation/bayesian_feature_fitting.py @@ -5,6 +5,16 @@ import pybop +""" +This example demonstrates how to use EP-BOLFI to parameterise a PyBaMM model using +a "feature"-based cost function. We use the term "feature" to describe a parameter +obtained from fitting either the data or a candidate solution to a simpler model. +Every evaluation of a feature-based cost function runs its own optimisation (based +on the simpler model) to identify the value of the feature. The aim is to minimise +the "feature distance" to identify the parameter values which produce a candidate +solution with a feature value as close as possible to that of the data. +""" + # Define model and parameter values model = pybamm.lithium_ion.SPMe() parameter_values = pybamm.ParameterValues("Chen2020") @@ -115,6 +125,7 @@ optim = pybop.EP_BOLFI(problem, options=options) result = optim.run() + # Plot the optimisation result pybop.plot.convergence(result, yaxis={"type": "log"}) pybop.plot.parameters(result, yaxis={"type": "log"}, yaxis2={"type": "log"}) diff --git a/examples/scripts/battery_parameterisation/bayesian_parameterisation_with_sober.py b/examples/scripts/battery_parameterisation/bayesian_parameterisation_with_sober.py new file mode 100644 index 000000000..41ccabaa6 --- /dev/null +++ b/examples/scripts/battery_parameterisation/bayesian_parameterisation_with_sober.py @@ -0,0 +1,196 @@ +import importlib.util +import sys + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from pybamm import citations, print_citations, DummySolver +from scipy.integrate import quad +from sober import setting_parameters +from torch import device, float64, set_default_dtype, tensor, zeros_like + +import pybop +from pybop.costs.feature_distances import indices_of +from pybop.optimisers.sober_basq_optimiser import ( + SOBER_BASQ_EPLFI, + SOBER_BASQ_EPLFI_Options, +) + +from pybop.models.lithium_ion.silicon_relaxation import SiliconRelaxation + +set_default_dtype(float64) +setting_parameters(device=device("cpu"), dtype=float64) + +seed = 0 + +if __name__ == "__main__": + data_index = 16 + + spec = importlib.util.spec_from_file_location( + "read_dataset", "../../data/Wycisk2024/read_dataset.py" + ) + read_dataset = importlib.util.module_from_spec(spec) + sys.modules["read_dataset"] = read_dataset + spec.loader.exec_module(read_dataset) + measurement = read_dataset.gitt_on_graphite_with_5_percent_silicon + + # pulses = measurement.get_subset(list(range(41, 82, 2))) + relaxations = measurement.get_subset(list(range(42, 83, 2))) + for i in range(len(relaxations)): + relaxations[i]["Time [s]"] = [ + t - relaxations[i]["Time [s]"][0] for t in relaxations[i]["Time [s]"] + ] + relaxations[i]["Voltage change [V]"] = [ + relaxations[i]["Voltage change [V]"][0] - u + for u in relaxations[i]["Voltage change [V]"] + ] + # The first timepoint is set to 0, which messes up the plots. + # The second timepoint is three orders of magnitude smaller than the + # next, which messes up the parameterization. + for i in range(len(relaxations)): + relaxations[i]["Time [s]"] = relaxations[i]["Time [s]"][2:] + relaxations[i]["Current function [A]"] = relaxations[i]["Current function [A]"][ + 2: + ] + relaxations[i]["Voltage change [V]"] = relaxations[i]["Voltage change [V]"][2:] + + t = np.asarray(relaxations[data_index]["Time [s]"]) + short_term_end = indices_of(t, 1e2)[0] + long_term_start = indices_of(t, 1e4)[0] + + dataset = pybop.Dataset( + { + "Time [s]": t, + "Current function [A]": np.asarray( + relaxations[data_index]["Current function [A]"] + ), + "Voltage [V]": np.asarray( + relaxations[data_index]["Voltage change [V]"] + ), + } + ) + len_data = len(t) + + short_term = pybop.MeanSquaredError( + dataset, + "Voltage [V]", + [1] * short_term_end + [0] * (len_data - short_term_end), + ) + mid_term = pybop.MeanSquaredError( + dataset, + "Voltage [V]", + [0] * short_term_end + + [1] * (long_term_start - short_term_end) + + [0] * (len_data - long_term_start), + ) + long_term = pybop.MeanSquaredError( + dataset, + "Voltage [V]", + [0] * long_term_start + [1] * (len_data - long_term_start), + ) + + unknowns = pybop.MultivariateParameters( + { + "Terminal voltage of observed relaxation [V]": pybop.Parameter( + initial_value=0.1, + bounds=[0.01, 0.2], + transformation=pybop.LogTransformation(), + ), + "Logarithmic slope of observed mechanical relaxation [V]": pybop.Parameter( + initial_value=0.01, + bounds=[0.001, 0.2], + transformation=pybop.LogTransformation(), + ), + "Exponential timescale for decay of observed mechanical relaxation [s]": pybop.Parameter( + initial_value=1e5, + bounds=[1e3, 1e7], + transformation=pybop.LogTransformation(), + ), + }, + distribution=pybop.MultivariateUniform( + np.asarray([[0.01, 0.2], [0.001, 0.2], [1e3, 1e7]]) + ), + ) + model = SiliconRelaxation() + simulator = pybop.pybamm.simulator.Simulator(model, model.default_parameter_values) + + + # Override the forced univariate Parameters + simulator.parameters = unknowns + problem = pybop.MetaProblem( + pybop.Problem(simulator, mid_term), pybop.Problem(simulator, long_term) + ) + + # Copy the MultivariateParameters to the meta-problem + problem.parameters = simulator.parameters + options = SOBER_BASQ_EPLFI_Options( + model_initial_samples=128, + seed=seed, + # disable_numpy_mode=True, + # parallelisation=False, + ep_iterations=2, + ep_total_dampening=0.5, + sober_iterations=4, + model_samples_per_iteration=64, + ep_integration_nodes=32, + integration_nodes=512, + batched_input=True, + ) + optim = SOBER_BASQ_EPLFI(problem, options=options) + result = optim.run() + + # Calculate the correlation matrix in addition to the full plot. + raw_taken_samples = result.posterior.distribution.distribution.dataset + raw_mean = np.mean(raw_taken_samples, axis=1) + raw_cov = np.cov(raw_taken_samples) + raw_std = np.var(raw_taken_samples, axis=1) ** 0.5 + raw_corr = (raw_cov / raw_std[:, None]) / raw_std[None, :] + fig_corr, ax_corr = plt.subplots(figsize=(3.75, 3)) + pybop.plot.correlation( + fig_corr, + ax_corr, + raw_corr, + names=["U(t=∞) [V]", "log-slope [V]", "relaxation timescale [s]"], + title="", + entry_color="white", + ) + + # Re-sample the posterior for the predictive posterior. + posterior_resamples = result.posterior.rvs(64, apply_transform=True) + posterior_resamples_pdf = result.posterior.pdf(posterior_resamples) + simulations = simulator.voltage_relaxation(posterior_resamples.T) + fig_pos, ax_pos = plt.subplots(figsize=(3 * 2**0.5, 3), layout="constrained") + norm = matplotlib.colors.Normalize( + posterior_resamples_pdf.min(), posterior_resamples_pdf.max() + ) + cmap = plt.get_cmap("viridis") + for pr, pr_pdf, u in zip( + posterior_resamples.T, posterior_resamples_pdf, simulations.T + ): + ax_pos.semilogx( + t, + u, + color=cmap(norm(pr_pdf)), + lw=0.8, + ls=":", + ) + ax_pos.semilogx( + tensor(t), + tensor(relaxations[data_index]["Voltage change [V]"]), + color="black", + lw=2, + label="experimental data", + ) + fig_pos.colorbar( + matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap), + ax=ax_pos, + label="Posterior PDF from KDE approximation", + ) + ax_pos.set_xlabel("Time since start of relaxation / s") + ax_pos.set_ylabel("Voltage change since start of relaxation / V") + ax_pos.legend() + + print_citations() + + plt.show() + diff --git a/examples/scripts/design_optimisation/bayesian_optimisation.py b/examples/scripts/design_optimisation/bayesian_optimisation.py new file mode 100644 index 000000000..d6efd904c --- /dev/null +++ b/examples/scripts/design_optimisation/bayesian_optimisation.py @@ -0,0 +1,295 @@ +# Miscellaneous imports for plotting, arithmetic, and statistics. + +# In this example, we use parallel processing. + +import matplotlib.pyplot as plt +import numpy as np +import pybamm +from matplotlib.ticker import PercentFormatter +from pybamm import print_citations +from sober import setting_parameters + +# Imports the SOBER interface and ensures that calculations are on CPU. +from torch import device, float64, set_default_dtype + +import pybop +from pybop.optimisers.sober_basq_optimiser import SOBER_BASQ, SOBER_BASQ_Options + +# Sets the precision from 32-bit to 64-bit and forces CPU over GPU calculations. +set_default_dtype(float64) +setting_parameters(dtype=float64, device=device("cpu")) +np.seterr(divide="ignore") + +""" +This example demonstrates how to call SOBER to perform a simple optimisation. +The example considers a heavily simplified operation cycle of a battery used +with a solar panel to provide energy at night. The optimisation task is to find +the right trade-off size between cost for additional capacity and the extended +life-time of the battery. Note: batteries degrade faster at extreme charge states. +""" + + +def day_night_cycle(oversize_factor): + """ + Generates a cycling protocol based on the fraction of unused cell, + and runs the SEI growth model with it. The cycling protocol emulates + a solar panel-coupled battery over 10 years. + + :param oversize_factor: + A number greater than 1, giving the ratio of cell capacity to + minimally required cell capacity. + :returns: + The SEI thickness over time. + """ + reference_current = (1.0 / 2.5) / oversize_factor + # With 1 hour per timestep, each cycle lasts 24 hours: discharge from + # 17:00 to 7:00, rest until 9:00, charge until 15:00, and rest until 17:00. + dt = 3600 + currents = np.asarray( + ( + [reference_current * 14 / 20] * 6 + + [0.0] * 2 + + [-reference_current * 6 / 20] * 14 + + [0.0] * 2 + ) + * 365 + * 10 + ) # roughly 10 years + timepoints = np.asarray(range(len(currents))) * dt + return timepoints, currents + + +def capacity_loss_cutoff(variables): + return 0.6 - variables["Capacity lost to SEI [C]"] / (3600 * variables["Nominal cell capacity [A.h]"]) + + +def set_up_simulator(reference_lifetime=-1): + """ + We need to set up the model twice. First we get the reference lifetime, + then we use it to set up the target relative to that reference. + """ + sei_model = pybamm.lithium_ion.SPM(options={"SEI": "VonKolzenberg2020", "surface form": "algebraic"}) + # Add two variables to the model: one for cutting off at the EOL, + # and one for describing the gain function to optimise. + sei_model.variables["Capacity lost to SEI [C]"] = ( + (sei_model.variables["Negative SEI thickness [m]"] - battery_parameters["Initial SEI thickness [m]"]) + * battery_parameters["Negative electrode surface area to volume ratio [m-1]"] + * battery_parameters["Electrode width [m]"] * battery_parameters["Electrode height [m]"] * battery_parameters["Negative electrode thickness [m]"] + / battery_parameters["SEI partial molar volume [m3.mol-1]"] + * F + ) + capacity_loss_termination = pybamm.step.CustomTermination(name="Capacity loss cut-off", event_function=capacity_loss_cutoff) + if reference_lifetime > 0: + sei_model.variables["Lifetime gained per oversize factor [d]"] = ( + sei_model.variables["Time [d]"] - reference_lifetime + ) / (sei_model.parameters["Oversize factor"] - 1) + protocol = pybamm.Experiment([ + "Charge at 0.28 C for 6 hours", + "Rest for 2 hours", + "Discharge at 0.12 C for 14 hours", + "Rest for 2 hours", + ] * 3650, period="1 hour", termination=[capacity_loss_termination]) + solar_battery_model = pybop.pybamm.Simulator(sei_model, battery_parameters, output_variables=["Voltage [V]"]) + return solar_battery_model + + + +if __name__ == "__main__": + battery_parameters = pybamm.ParameterValues("Marquis2019") + + F = 96485.33212 + # We append parameters from a cell we know the SEI parameters for. + battery_parameters.update( + { + "Negative electrode surface area to volume ratio [m-1]": 3 * (1 - battery_parameters["Negative electrode porosity"]) / battery_parameters["Negative particle radius [m]"], + "Electrolyte diffusivity [m2.s-1]": 2.8e-10, + "Electrode width [m]": (1 / 24) ** 0.5, + "Electrode height [m]": (1 / 24) ** 0.5, + "Nominal cell capacity [A.h]": 0.05, + "SEI reaction exchange current density [A.m-2]": 1e-5, + "SEI lithium ion conductivity [S.m-1]": 1e-8, + "Initial concentration in negative electrode [mol.m-3]": ( + 0.11 + * battery_parameters[ + "Maximum concentration in negative electrode [mol.m-3]" + ] + ), + "Initial concentration in positive electrode [mol.m-3]": ( + 0.83 + * battery_parameters[ + "Maximum concentration in positive electrode [mol.m-3]" + ] + ), + # "SEI growth transfer coefficient": 0.22, # currently not supported in PyBaMM + "Inital SEI thickness [m]": 2e-9, + "Tunneling distance for electrons [m]": 2.05e-9, # tunneling within SEI + "SEI partial molar volume [m3.mol-1]": 1.078e-5, + "SEI lithium interstitial diffusivity [m2.s-1]": 1e-15, + # "SEI relative permittivity": 131, # not required for growth + "SEI open-circuit potential [V]": 17400 / F, + # "Anion transference number in SEI": 1 - 0.063, # t_plus is 0.063, not required for growth + "SEI porosity": 0.1, + # "SEI Bruggeman coefficient": 4.54, # not required for growth + "Relative capacity cut-off for End-Of-Life": 0.4, + # Change cell capacity by adjusting cross-section area. + "Electrode height [m]": battery_parameters["Electrode height [m]"] * pybamm.Parameter("Oversize factor"), + "Oversize factor": "[input]", + }, + check_already_exists=False, + ) + + pybop_prior = pybop.MultivariateParameters( + {"Adjustment factor": pybop.Parameter(initial_value=1.1, bounds=[1.0, 1.2])}, + distribution=pybop.MultivariateUniform(np.asarray([[1.0, 1.2 - 1.0]])), # bug work-around; change to actual [1.0, 1.2] bounds when PR862 is merged + ) + # In this simple example, we first plot the whole target function. + # A note, as this is a common point of confusion: this plot solves + # the optimization already, as we can see the optimum point. + # The point of this example is to show the application of the + # optimiser in an easily verifiable example. As soon as many + # variables are to be optimised at once, say 5 or more, or the + # landscape of the target function becomes much more complex, + # we use to optimiser to get away with much sparser evaluations. + reference_lifetime = set_up_simulator().solve()["Time [d]"].entries[-1] + simulator = set_up_simulator(reference_lifetime) + cost = DesignCost(target="Time [d]") + oversize_factors = np.asarray([1 + 0.002 * i for i in range(1, 101)]) + eol_days = [reference_lifetime] + for oversize_factor in oversize_factors: + eol_days.append(cost.evaluate(solar_battery_model.solve(inputs={"Oversize factor": oversize_factor}))) + oversize_factors = np.append([1.0], oversize_factors) + eol_days = np.asarray(eol_days) + + fig_kde, ax_kde = plt.subplots(figsize=(3 * 2**0.5, 3), layout="constrained") + ax_eol = ax_kde.twinx() + eol_plot = ax_eol.plot( + oversize_factors, eol_days, label="Time to End-Of-Life / d" + )[0] + gain = [float("NaN")] + list( + (eol_days[1:] - eol_days[0]) / (oversize_factors[1:] - 1) + ) + gain_plot = ax_eol.plot( + oversize_factors, gain, label="Gain per extra capacity / d" + )[0] + ax_kde.set_xlabel("Oversize factor") + ax_eol.set_ylabel("Time / d") + + # We have seen that the optimum ratio of battery lifetime gained to + # battery oversizing is achieved at ~ 10.3% oversizing fraction. + # Now we showcase how to obtain this result with SOBER instead. + # We utilise the same interface as for the parameterization, and for + # most of its arguments we refer you to its documentation. + # In the special case of optimization, we set the 'data' argument to + # a suitably-shaped 0,such that "target function - data" is just + # the target function, and set 'maximize' to True. + # Since it is needed for EOL gain calculation, we also pass the + # reference EOL, which will be passed on to the target function. + """ + sober_wrapper = SoberWrapper( + calculate_eol_gain, + tensor([0]), + model_initial_samples=16, + bounds=tensor([[1.0], [1.2]]), + prior='Uniform', + maximize=True, + seed=0, + names=["Oversize factor"], + true_optimum=tensor([1.103]), + offset=eol_indices[0] / 24 + ) + """ + + # We now invoke SOBER to explore the target function efficiently. + # Its settings, for which we refer you to its documentation, are + # best found by trial-and-error and a coarse initial guess about + # the complexity of the target function. Its results are stored + # in the interface instance, which we will access later. + """ + sober_wrapper.run_SOBER( + sober_iterations=7, + model_samples_per_iteration=16, + visualizations=False, + verbose=True + ) + """ + + # We now invoke BASQ to assess the quality with which SOBER has + # explored the target function. Its settings, for which we refer + # you to its documentation, are best found by trial-and-error and + # a coarse initial guess about the complexity of the target + # function. We get five return values: + # 1. samples from the probability distribution that SOBER generated + # as a (faster) surrogate to the original target function, + # 2. the optimal point in terms of the maximum value of the + # surrogate, called the Maximum A Posteriori (MAP) point, + # 3. the optimal point that has been evaluated on the original + # target function, + # 4. the SOBER approximation quality criterion, the expected log + # marginal likelihood (lower is better, scale is relative), + # 5. the SOBER approximation quality criterion quality criterion, + # i.e., the self-assessment about the accuracy with which the + # expected log marginal likelihood was calculated, expressed + # in terms of the variance of the log marginal likelihood. + """ + ( + taken_samples, + MAP, + best_observed, + log_expected_marginal_likelihood, + log_approx_variance_marginal_likelihood + ) = sober_wrapper.run_BASQ( + integration_nodes=128, + visualizations=False, + verbose=True + ) + """ + + cost = pybop.DesignCost("EOL [d]") + cost._target_data = np.asarray([0]) + pybop_problem = pybop.Problem(solar_battery_model, cost) + pybop_problem.parameters = pybop_prior + pybop_options = SOBER_BASQ_Options( + model_initial_samples=16, + maximise=True, + sober_iterations=7, + model_samples_per_iteration=16, + integration_nodes=128, + ) + sober_basq_wrapper = SOBER_BASQ(pybop_problem, pybop_options) + pybop_result = sober_basq_wrapper.run() + + # sober_wrapper = sober_basq_wrapper.optim + + # We have seen visualizations form 'run_SOBER' and 'run_BASQ' about + # their internal states. For a more intuitive visualization, we + # employ the so-called predictive posterior. Rather than showing + # the probability distribution of the model parameter values, we + # plot the model realizations for a representative sample of it. + # We will use the samples SOBER took and the samples BASQ generated + # for plotting a so-called Kernel Density Estimate (KDE). If you are + # not familiar with this, think of it as a spruced up histogram. + eval_kde = np.linspace(1.0, 1.2, 101) + post_approx = pybop_result.posterior.pdf(eval_kde) + post_norm = sum(post_approx) * (1.2 - 1.0) + post_approx /= post_norm + kde_plot = ax_kde.plot( + eval_kde, post_approx, label="Posterior for optimal sizing", ls="--" + )[0] + ax_kde.set_xlabel("Oversize factor") + ax_kde.set_ylabel("Posterior probability density") + ax_kde.xaxis.set_major_formatter(PercentFormatter(xmax=1.0)) + plots_for_legend = [eol_plot, gain_plot, kde_plot] + fig_kde.legend( + plots_for_legend, + [p.get_label() for p in plots_for_legend], + loc="outside lower center", + # borderpad=0.3, + # handlelength=0.5, + # handletextpad=0.3, + # borderaxespad=0.0, + # columnspacing=0.5, + ) + + print_citations() + + plt.show() diff --git a/examples/scripts/model_selection/counting_kneepoints.py b/examples/scripts/model_selection/counting_kneepoints.py new file mode 100644 index 000000000..3e669edc5 --- /dev/null +++ b/examples/scripts/model_selection/counting_kneepoints.py @@ -0,0 +1,156 @@ +import importlib.util +import sys + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import torch +from pybamm import print_citations +from scipy.stats import gaussian_kde +from torch import tensor + +import pybop +from pybop.optimisers.sober_basq_optimiser import SOBER_BASQ, SOBER_BASQ_Options + +from pybop.models.kneepoints import KneepointModel +from pybop.plot.predictive import predictive + +import plotly.graph_objects as go + +# For the machine this was tested on; comment out for default behaviour. +import plotly.io as pio +pio.renderers.default = "browser" + + +""" +class MeanSquaredErrorPyTorch(pybop.costs.error_measures.ErrorMeasure): + + def __call__(self, r: torch.Tensor) -> torch.Tensor: + e = torch.sum(r**2, axis=1)**0.5 + return e +""" + + +if __name__ == "__main__": + data_index = 0 + + spec = importlib.util.spec_from_file_location( + "read_dataset", "../../data/Baumhofer2014/read_dataset.py" + ) + read_dataset = importlib.util.module_from_spec(spec) + sys.modules["read_dataset"] = read_dataset + spec.loader.exec_module(read_dataset) + measurements = read_dataset.degradation_data + + fig, ax = plt.subplots(figsize=(2.4 * 2**0.5, 2.4), constrained_layout=True) + for m in measurements: + ax.plot(m["Time [s]"], m["Capacity fade"]) + + # Cast non-standard dtypes into NumPy floats to avoid PyTorch errors. + t = np.ndarray.astype(measurements[data_index]["Time [s]"].T[0], np.float64)[1:] + t[5] = t[4] + 1 + data = np.ndarray.astype( + measurements[data_index]["Capacity fade"].T[0], np.float64 + )[1:] + dataset = pybop.Dataset({"Time [s]": t, "Capacity fade": data}) + """ + ax.plot( + t.cpu(), + kneepoint_model( + torch.tensor([[0.0002, 1500, 0.0008]]), + t + )[0].cpu(), + color='black', + lw=2, + label='model' + ) + """ + ax.set_xlabel("Full cycles") + ax.set_ylabel("Relative remaining capacity") + # ax.legend() + plt.show() + + for n_kneepoints, mean, bounds, names in zip( + (1, 2), + ( + np.array([0.0002, 1500, 0.0008]), + np.array([0.0002, 1500, 0.0008, 2000, 0.001]), + ), + ( + np.array([[0.00001, 0.001], [200, 2000], [0.00001, 0.01]]), + np.array( + [ + [0.00001, 0.001], + [200, 2000], + [0.00001, 0.01], + [700, 2500], + [0.00001, 0.01], + ] + ), + ), + ( + [ + "1st degr. rate [Capacity/Cycle]", + "1st kneepoint [Cycle]", + "2nd degr. rate [Capacity/Cycle]", + ], + [ + "1st degr. rate [Capacity/Cycle]", + "1st kneepoint [Cycle]", + "2nd degr. rate [Capacity/Cycle]", + "2nd kneepoint [Cycle]", + "3rd degr. rate [Capacity/Cycle]", + ], + ), + strict=False, + ): + initial_values = np.exp(0.5 * (np.log(bounds.T[0]) + np.log(bounds.T[1]))) + pybop_prior = pybop.MultivariateParameters( + { + n: pybop.Parameter( + initial_value=i, bounds=b, transformation=pybop.LogTransformation() + ) + for n, i, b in zip(names, initial_values, bounds, strict=False) + }, + distribution=pybop.MultivariateGaussian(mean=mean, bounds=bounds), + ) + simulator = KneepointModel(pybop_prior, tensor(t), n_kneepoints) + # Override the forced univariate Parameters + simulator.parameters = pybop_prior + cost = pybop.MeanSquaredError(dataset, "Capacity fade") + cost._target_data = np.asarray([0]) + pybop_problem = pybop.Problem(simulator, cost) + # Copy the MultivariateParameters to the meta-problem + pybop_problem.parameters = simulator.parameters + pybop_options = SOBER_BASQ_Options( + model_initial_samples=256, + sober_iterations=12, + model_samples_per_iteration=64, + integration_nodes=256, + batched_input=True, + ) + sober_basq_wrapper = SOBER_BASQ(pybop_problem, pybop_options) + pybop_result = sober_basq_wrapper.run() + kneepoint_pdf_x_eval = np.linspace(np.log(t[1]), np.log(t[-1]), 201) + kneepoint_pdf_y = np.zeros_like(kneepoint_pdf_x_eval) + for i in range(n_kneepoints): + kneepoint_pdf_y += pybop_result.posterior.distribution.distribution.marginal(1 + 2 * i).pdf(kneepoint_pdf_x_eval) + kneepoint_pdf_x_plot = np.exp(kneepoint_pdf_x_eval) + fig = predictive( + pybop_result, + simulator=lambda p: simulator(tensor(p.T)).T.detach().cpu().numpy(), + number_of_traces=64, + dataset_y="Capacity fade", + data_legend_entry="degradation data", + rvs_legend_entry="candidate fits", + pdf_plot=[kneepoint_pdf_x_plot, kneepoint_pdf_y], + pdf_label="PDF for knee points", + show=True, + xaxis_title="Cycles", + yaxis_title="Capacity", + yaxis2={'title': "PDF for knee points", 'overlaying': 'y', 'side': 'right'}, + yaxis_range=[-0.1, 1.1], + title=str(n_kneepoints) + "-knee point model", + ) + print_citations() + diff --git a/examples/scripts/model_selection/impedance_data_plot.py b/examples/scripts/model_selection/impedance_data_plot.py new file mode 100644 index 000000000..dc0e01f9f --- /dev/null +++ b/examples/scripts/model_selection/impedance_data_plot.py @@ -0,0 +1,84 @@ +import json + +import matplotlib.pyplot as plt +from ep_bolfi.utility.dataset_formatting import read_parquet_table +from ep_bolfi.utility.visualization import nyquist_plot + +x_lims = { + "monolayer_17_microns": (0, 300), + "porous_42_microns": (0, 75), + "porous_80_microns": (0, 75), +} +y_lims = { + "monolayer_17_microns": (0, 200), + "porous_42_microns": (0, 50), + "porous_80_microns": (0, 50), +} + +with open("../../data/Gunther2025/impedance_ocv_alignment.json") as f: + alignment = json.load(f) + +data = {} +""" +for filename, plotname in [ + ("monolayer_17_microns", "17 µm"), + ("porous_42_microns", "42 µm"), + ("porous_80_microns", "80 µm") +]: + soc = alignment[filename]["Positive electrode SOC [-]"] + ocv = alignment[filename]["OCV [V]"] + data[filename] = read_parquet_table(filename + ".parquet", 'impedance') + for (direction, index) in (("delithiation", 0), ("lithiation", 1)): + fig, ax = plt.subplots(figsize=(5, 3)) + nyquist_plot( + fig, + ax, + data[filename].frequencies[index::2], + data[filename].complex_impedances[index::2], + title_text="NMC " + plotname + " - " + direction + " direction", + legend_text=['{:3.3g} V'.format(o) for o in ocv[index::2]] + ) + ax.set_xlim(x_lims[filename]) + ax.set_ylim(y_lims[filename]) + plt.draw() + legend_position = ax.get_legend().get_bbox_to_anchor().transformed(ax.transAxes.inverted()) + legend_position.x0 -= 1.25 + legend_position.x1 -= 1.25 + legend_position.y0 += 0.15 + legend_position.y1 += 0.15 + ax.get_legend().set_bbox_to_anchor(legend_position, transform=ax.transAxes) + fig.tight_layout() + fig.savefig(filename + "_" + direction + ".pdf", bbox_inches='tight', pad_inches=0.0) +""" +filename = "18650_LG_3500_MJ1_EIS_anode_discharge" +soc = alignment[filename]["Negative electrode SOC [-]"] +data[filename] = read_parquet_table( + "../../data/Kuhn2026/" + filename + ".parquet", "impedance" +) +fig, ax = plt.subplots(figsize=(5, 3)) +nyquist_plot( + fig, + ax, + data[filename].frequencies, + data[filename].complex_impedances, + title_text="graphite - delithiation", + legend_text=[f"{int(round(100 * s)):d} %" for s in soc], +) +handles = ax.get_legend().legend_handles +labels = [t._text for t in ax.get_legend().texts] +ax.legend(handles, labels, fontsize=6, ncols=2) +ax.set_xlim((0, 5)) +ax.set_ylim((0, 5)) +plt.draw() +legend_position = ( + ax.get_legend().get_bbox_to_anchor().transformed(ax.transAxes.inverted()) +) +legend_position.x0 -= 1.25 +legend_position.x1 -= 1.25 +legend_position.y0 += 0.15 +legend_position.y1 += 0.15 +ax.get_legend().set_bbox_to_anchor(legend_position, transform=ax.transAxes) +fig.tight_layout() +fig.savefig(filename + "_delithiation.pdf", bbox_inches="tight", pad_inches=0.0) + +plt.show() diff --git a/examples/scripts/model_selection/impedance_model_selection.py b/examples/scripts/model_selection/impedance_model_selection.py new file mode 100644 index 000000000..9129ccee7 --- /dev/null +++ b/examples/scripts/model_selection/impedance_model_selection.py @@ -0,0 +1,755 @@ +import importlib.util +import json +import sys +from contextlib import redirect_stdout +from multiprocessing import Pool + +import matplotlib.pyplot as plt +import numpy as np +import pybamm +import pybammeis +from botorch.acquisition.active_learning import qNegIntegratedPosteriorVariance +from ep_bolfi.models.solversetup import spectral_mesh_pts_and_method +from ep_bolfi.utility.dataset_formatting import read_parquet_table +from ep_bolfi.utility.fitting_functions import fit_drt +from ep_bolfi.utility.preprocessing import find_occurrences +from ep_bolfi.utility.visualization import ( + interactive_impedance_model, + nyquist_plot, +) +from scipy import constants +from scipy.stats import gaussian_kde +from sober import SoberWrapper +from torch import exp, log, pi, tensor, zeros + +spec_gunther = importlib.util.spec_from_file_location( + "gunther", "../../data/Gunther2025/parameters.py" +) +gunther = importlib.util.module_from_spec(spec_gunther) +sys.modules["gunther"] = gunther +spec_gunther.loader.exec_module(gunther) +spec_kuhn = importlib.util.spec_from_file_location( + "kuhn", "../../data/Kuhn2026/parameters.py" +) +kuhn = importlib.util.module_from_spec(spec_kuhn) +sys.modules["kuhn"] = kuhn +spec_kuhn.loader.exec_module(kuhn) + +parameter_ranges = { + "pouch_SPM": { + "Positive particle diffusivity [m2.s-1]": (1e-16, 1e-14), + "Positive electrode double-layer capacity [F.m-2]": (0.1, 0.5), + "Positive electrode exchange-current density [A.m-2]": (0.1, 1.0), + }, + "pouch_DFN": { + "Positive electrode Bruggeman coefficient": (1, 2.5), + "Positive electrode conductivity [S.m-1]": (0.01, 10), + "Electrolyte diffusivity [m2.s-1]": (1e-10, 1e-9), + "Positive particle diffusivity [m2.s-1]": (1e-16, 1e-14), + "Positive electrode double-layer capacity [F.m-2]": (0.1, 0.5), + "Positive electrode exchange-current density [A.m-2]": (0.1, 1.0), + }, + "LG_MJ1_variant_1": { + "Positive electrode double-layer capacity [F.m-2]": (7e-4, 7e-3), + "Positive electrode exchange-current density [A.m-2]": (2e-2, 8e-2), + "SEI relative permittivity": (100.0, 200.0), + "SEI Bruggeman coefficient": (4.3, 4.8), + }, + "LG_MJ1_variant_2": { + "Positive electrode double-layer capacity [F.m-2]": (1e-2, 1e-1), + "Positive electrode exchange-current density [A.m-2]": (2e-2, 8e-2), + "SEI relative permittivity": (10.0, 50.0), + "SEI Bruggeman coefficient": (4.3, 4.8), + }, +} + +no_transform = (lambda x: x, lambda x: x) +log_transform = (lambda x: log(x), lambda x: exp(x)) +parameter_transforms = { + "pouch_SPM": { + "Positive particle diffusivity [m2.s-1]": log_transform, + "Positive electrode double-layer capacity [F.m-2]": log_transform, + "Positive electrode exchange-current density [A.m-2]": log_transform, + }, + "pouch_DFN": { + "Positive electrode Bruggeman coefficient": no_transform, + "Positive electrode conductivity [S.m-1]": log_transform, + "Electrolyte diffusivity [m2.s-1]": log_transform, + "Positive particle diffusivity [m2.s-1]": log_transform, + "Positive electrode double-layer capacity [F.m-2]": log_transform, + "Positive electrode exchange-current density [A.m-2]": log_transform, + }, + "LG_MJ1_variant_1": { + "Positive electrode double-layer capacity [F.m-2]": log_transform, + "Positive electrode exchange-current density [A.m-2]": log_transform, + "SEI relative permittivity": log_transform, + "SEI Bruggeman coefficient": no_transform, + }, + "LG_MJ1_variant_2": { + "Positive electrode double-layer capacity [F.m-2]": log_transform, + "Positive electrode exchange-current density [A.m-2]": log_transform, + "SEI relative permittivity": log_transform, + "SEI Bruggeman coefficient": no_transform, + }, +} +np_log_transform = (lambda x: np.log(x), lambda x: np.exp(x)) +parameter_transforms_numpy = { + "pouch_SPM": { + "Positive particle diffusivity [m2.s-1]": np_log_transform, + "Positive electrode double-layer capacity [F.m-2]": np_log_transform, + "Positive electrode exchange-current density [A.m-2]": np_log_transform, + }, + "pouch_DFN": { + "Positive electrode Bruggeman coefficient": no_transform, + "Positive electrode conductivity [S.m-1]": np_log_transform, + "Electrolyte diffusivity [m2.s-1]": np_log_transform, + "Positive particle diffusivity [m2.s-1]": np_log_transform, + "Positive electrode double-layer capacity [F.m-2]": np_log_transform, + "Positive electrode exchange-current density [A.m-2]": np_log_transform, + }, + "LG_MJ1_variant_1": { + "Positive electrode double-layer capacity [F.m-2]": np_log_transform, + "Positive electrode exchange-current density [A.m-2]": np_log_transform, + "SEI relative permittivity": np_log_transform, + "SEI Bruggeman coefficient": no_transform, + }, + "LG_MJ1_variant_2": { + "Positive electrode double-layer capacity [F.m-2]": np_log_transform, + "Positive electrode exchange-current density [A.m-2]": np_log_transform, + "SEI relative permittivity": np_log_transform, + "SEI Bruggeman coefficient": no_transform, + }, +} + + +def Z_SEI(f, parameters): + """ + Impedance of the Solid-Electrolyte Interphase (Single2019). + + :param f: + An array of the frequencies to evaluate. + :param parameters: + A dictionary of the model parameters. + :returns: + The evaluated impedances as an array. + """ + + εₙ = parameters["Negative electrode porosity"] + βₙ = parameters["Negative electrode Bruggeman coefficient (electrolyte)"] + Dₑ = parameters["Electrolyte diffusivity [m2.s-1]"] + ε_SEI = parameters["SEI porosity"] + β_SEI = parameters["SEI Bruggeman coefficient"] + nₚ = parameters["Stoichiometry of cation in electrolyte salt dissociation"] + nₙ = parameters["Stoichiometry of anion in electrolyte salt dissociation"] + zₚ_salt = parameters["Charge number of cation in electrolyte salt dissociation"] + zₙ_salt = parameters["Charge number of anion in electrolyte salt dissociation"] + ρₑ = parameters["Mass density of electrolyte [kg.m-3]"] + ρₑ_plus = parameters["Mass density of cations in electrolyte [kg.m-3]"] + M_N = parameters["Molar mass of electrolyte solvent [kg.mol-1]"] + c_N = parameters["Solvent concentration [mol.m-3]"] + ρ_N = M_N * c_N + v_N = parameters["Partial molar volume of electrolyte solvent [m3.mol-1]"] + tilde_ρ_N = M_N / v_N + one_plus_dlnf_dlnc = parameters["Thermodynamic factor"] + Lₙ = parameters["Negative electrode thickness [m]"] + Lₛ = parameters["Separator thickness [m]"] + Lₚ = parameters["Positive electrode thickness [m]"] + L_electrolyte_for_SEI_model = (((Lₙ + Lₚ) / 2) + Lₛ) / 2 + L_SEI = parameters["SEI thickness [m]"] + t_SEI_minus = parameters["Anion transference number in SEI"] + permittivity_SEI = parameters["SEI relative permittivity"] + F = constants.physical_constants["Faraday constant"][0] + ɛ_0 = constants.physical_constants["vacuum electric permittivity"][0] + C_SEI = ɛ_0 * permittivity_SEI / L_SEI + κ_SEI = parameters["SEI ionic conductivity [S.m-1]"] + R_SEI = L_SEI / (ε_SEI**β_SEI * κ_SEI) + + # Notations refer to Single2019. + k_electrolyte = (1 - 1j) * np.sqrt(εₙ ** (-βₙ) * 2 * pi * f / (2 * Dₑ)) + k_SEI = (1 - 1j) * np.sqrt(ε_SEI ** (1 - β_SEI) * 2 * pi * f / (2 * Dₑ)) + Theta = ( + -(nₚ + nₙ) + / (nₚ * nₙ) + / (zₚ_salt * zₙ_salt * F**2) + * ρₑ**2 + / (ρ_N * tilde_ρ_N) + * one_plus_dlnf_dlnc + ) + Psi = 1 - ε_SEI ** ((1 + β_SEI) / 2) * np.tan( + k_electrolyte * L_electrolyte_for_SEI_model + ) * np.tan(k_SEI * L_SEI) + Z_D_SEI = ( + L_SEI + * Theta + / (Dₑ * ε_SEI**β_SEI) + * (t_SEI_minus - ρₑ_plus / ρₑ) ** 2 + * np.tan(k_SEI * L_SEI) + / (Psi * k_SEI * L_SEI) + ) + + return 1 / (2 * pi * 1j * f * C_SEI + 1 / (R_SEI + Z_D_SEI)) + + +def preprocess_data(data_index, electrode, cell_name): + with open("../../data/Gunther2025/impedance_ocv_alignment.json") as f: + ocv_alignment = json.load(f)[cell_name] + raw_data_index = ocv_alignment["indices"].index(data_index) + soc = ocv_alignment[electrode.capitalize() + " electrode SOC [-]"][raw_data_index] + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge": + directory = "../../Kuhn2026/" + else: + directory = "../../Gunther2025/" + data = read_parquet_table(directory + cell_name + ".parquet", "impedance") + + frequencies = tensor(data.frequencies[raw_data_index]) + impedances = tensor(data.complex_impedances[raw_data_index]) + + return frequencies, impedances, soc + + +########################################################## +# Simulator used for visualization and parameterization. # +########################################################## + + +def composed_model( + parameters, + angular_frequencies, + working_electrode="both", + electrolyte=True, + sei=False, + three_electrode=None, + reference_electrode_location=0.5, +): + # "surface form": "differential" enables surface capacitance. + model_options = { + "surface form": "differential", + "working electrode": working_electrode, + } + if electrolyte: + model = pybamm.lithium_ion.DFN(options=model_options) + else: + model = pybamm.lithium_ion.SPM(options=model_options) + discretization = { + "order_s_n": 10, + "order_s_p": 10, + "order_e": 10, + "volumes_e_n": 1, + "volumes_e_s": 1, + "volumes_e_p": 1, + "halfcell": False, + } + eis_sim = pybammeis.EISSimulation( + model, + pybamm.ParameterValues(dict(parameters)), + None, # geometry + *spectral_mesh_pts_and_method(**discretization), + three_electrode=three_electrode, + reference_electrode_location=reference_electrode_location, + ) + impedance = eis_sim.solve(angular_frequencies) + # Put the SEI analytic model in series to DFN with capacitance. + if sei: + impedance += Z_SEI(angular_frequencies, parameters) + return impedance + + +def simulator( + trial, + angular_frequencies, + cell_name, + soc, + electrolyte=True, + variable_parameters=None, +): + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge": + parameters = kuhn.get_parameters_with_switched_electrodes() + sei = True + working_electrode = "both" + three_electrode = "positive" + reference_electrode_location = 0.5 + else: + parameters = gunther.get_parameters_for_cell(cell_name) + sei = False + working_electrode = "positive" + three_electrode = "positive" + reference_electrode_location = 1.0 + + if variable_parameters is None: + variable_parameters = [] + + parameters["Positive electrode SOC"] = soc + for i, key in enumerate(variable_parameters): + parameters[key] = trial[i].item() + + return composed_model( + parameters, + angular_frequencies.detach().numpy(), + working_electrode, + electrolyte, + sei, + three_electrode, + reference_electrode_location, + ) + + +def drt_simulator( + trial, + angular_frequencies, + cell_name, + soc, + lambda_value=-2.0, + num_data_peaks=None, + electrolyte=True, + variable_parameters=None, +): + if variable_parameters is None: + variable_parameters = [] + + parallel_arguments = [ + (t, angular_frequencies, cell_name, soc, electrolyte, variable_parameters) + for t in trial + ] + with Pool() as p: + impedances = p.starmap(simulator, parallel_arguments) + parallel_arguments = [(angular_frequencies, i, lambda_value) for i in impedances] + with Pool() as p: + drt_results = p.starmap(fit_drt, parallel_arguments) + results = [] + for drt_tau, drt_resistance, _ in drt_results: + num_data_peaks = num_data_peaks or len(drt_tau) + results.append( + tensor( + [ + [[dr, dt]] * num_data_peaks + for dr, dt in zip(drt_tau, drt_resistance, strict=False) + ] + ).log() + ) + return results + + +def manual_model_assessment( + frequencies, + impedances, + parameters, + unknowns, + transform_unknowns, + working_electrode, + electrolyte, + sei, + three_electrode, + reference_electrode_location, + lambda_value, +): + interactive_impedance_model( + frequencies.detach().numpy(), + impedances.detach().numpy(), + parameters, + unknowns=unknowns, + transform_unknowns=transform_unknowns, + model=lambda par, freq: composed_model( + par, + freq / 1j, + working_electrode, + electrolyte, + sei, + three_electrode, + reference_electrode_location, + ), + lambda_value=lambda_value, + ) + plt.show() + + +def automated_model_assessment( + frequencies, + impedances, + mean, + bounds, + transforms, + cell_name, + soc, + electrolyte, + variable_parameters, + lambda_value, +): + drt_tau, drt_resistance, drt = fit_drt( + frequencies, impedances, lambda_value or -2.0 + ) + drt_data = tensor(list(zip(drt_tau, drt_resistance, strict=False))).log() + num_data_peaks = len(drt_tau) + sober_wrapper = SoberWrapper( + model=drt_simulator, + data=drt_data, + model_initial_samples=48, + mean=mean, + bounds=bounds, + prior="Uniform", + use_bolfi=False, + transforms=transforms, + disable_numpy_mode=True, + parallelization=False, + visualizations=False, + names=variable_parameters, + angular_frequencies=frequencies, + cell_name=cell_name, + soc=soc, + lambda_value=drt.lambda_value, + num_data_peaks=num_data_peaks, + electrolyte=electrolyte, + variable_parameters=variable_parameters, + ) + + # Re-define the distance function to counter the effect of less + # simulator DRT peaks automatically lessening the error. + # These are two versions to choose from. Taking the distance to all + # data peaks is more stable, but has weak optima. Taking only the + # distance to the nearest data peak leads to better convergence. + + def distance_function_all(observations): + # Distance of every simulator peak to every data peak. + num_sim_peaks = observations.shape[1] + return ((observations - sober_wrapper.data) * sober_wrapper.weights).view( + observations.shape[0], -1 + ).norm(dim=1).to( + device=sober_wrapper.tm.device, dtype=sober_wrapper.tm.dtype + ) / (num_sim_peaks / num_data_peaks) ** 0.5 + + def distance_function_nearest(observations): + # Distance of every simulator peak to its nearest data peak. + num_sim_peaks = observations.shape[1] + peak_distances = ( + (observations - sober_wrapper.data) * sober_wrapper.weights + ).norm(dim=3) / (num_sim_peaks / num_data_peaks) ** 0.5 + distances = zeros(observations.shape[0]).to( + device=sober_wrapper.tm.device, dtype=sober_wrapper.tm.dtype + ) + for i in range(len(distances)): + pd_entry = peak_distances[i].tolist() + for _ in range(num_sim_peaks): + nearest = tensor(pd_entry).argmin() + distances[i] += pd_entry[nearest // num_data_peaks][ + nearest % num_data_peaks + ] + pd_entry.pop(nearest // num_data_peaks) + return distances + + sober_wrapper.distance_function = distance_function_nearest + + sober_wrapper.run_SOBER( + sober_iterations=1, + model_samples_per_iteration=48, + acquisition_function=None, + visualizations=False, + verbose=True, + ) + + for _ in range(18): + raw_taken_samples = sober_wrapper.run_BASQ( + integration_nodes=3 ** len(variable_parameters), + visualizations=False, + return_raw_samples=True, + verbose=True, + )[0] + acquisition_function = qNegIntegratedPosteriorVariance( + sober_wrapper.surrogate_model, raw_taken_samples + ) + sober_wrapper.run_SOBER( + sober_iterations=1, + model_samples_per_iteration=48, + acquisition_function=lambda x: acquisition_function(x.unsqueeze(1)), + visualizations=False, + verbose=True, + ) + basq_results = sober_wrapper.run_BASQ( + integration_nodes=3 ** len(variable_parameters), + visualizations=False, + return_raw_samples=True, + verbose=True, + ) + return sober_wrapper, basq_results + + +def visualize_automated_assessment( + sober_wrapper, + basq_results, + frequencies, + cell_name, + prior_name, + soc, + lithiation, + electrolyte, + variable_parameters, +): + ( + raw_taken_samples, + MAP, + best_observed, + log_expected_marginal_likelihood, + log_approx_variance_marginal_likelihood, + ) = basq_results + filename = ( + cell_name + + "_" + + prior_name + + "_soc_" + + f"{soc:.3g}" + + "_" + + ("lithiation" if lithiation else "delithiation") + ) + with open(filename + "_posterior.json", "w") as f: + json.dump( + { + "variable_parameters": variable_parameters, + "diagonalization": sober_wrapper.diagonalization.tolist(), + "bounds": sober_wrapper.bounds.tolist(), + "diag_order": sober_wrapper.diag_order, + "raw_taken_samples": raw_taken_samples.numpy().tolist(), + }, + f, + ) + kde_posterior = gaussian_kde(raw_taken_samples.numpy().T) + raw_posterior_resamples = kde_posterior.resample(32) + posterior_resamples_pdf = kde_posterior(raw_posterior_resamples) + alphas = posterior_resamples_pdf * (0.5 / posterior_resamples_pdf.max()) + posterior_resamples = sober_wrapper.reverse_transform( + sober_wrapper.denormalize_input(tensor(raw_posterior_resamples.T)) + ) + parallel_arguments = [ + (trial, frequencies, cell_name, soc, electrolyte, variable_parameters) + for trial in posterior_resamples + ] + with Pool() as p: + simulations = p.starmap(simulator, parallel_arguments) + + fig_imp, ax_imp = plt.subplots(figsize=(4 * 2**0.5, 4)) + if electrolyte: + title_text = cell_name + " DFN fit" + else: + title_text = cell_name + " SPM fit" + nyquist_plot( + fig_imp, + ax_imp, + frequencies.detach().numpy(), + impedances.detach().numpy(), + ls=":", + lw=2, + title_text=title_text, + legend_text="experiment", + ) + nyquist_plot( + fig_imp, + ax_imp, + frequencies.detach().numpy(), + simulator(MAP, frequencies, cell_name, soc, electrolyte, variable_parameters), + title_text=title_text, + legend_text="optimal fit", + add_frequency_colorbar=False, + ) + for simulated_impedances, alpha in zip(simulations, alphas, strict=False): + nyquist_plot( + fig_imp, + ax_imp, + frequencies.detach().numpy(), + simulated_impedances, + lw=1, + alpha=alpha, + title_text=title_text, + legend_text=None, + add_frequency_colorbar=False, + ) + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge": + ax_imp.set_xlim(0.8, 4.8) + ax_imp.set_ylim(0.0, 2.0) + """ + elif cell_name == "monolayer_17_microns": + ax_imp.set_xlim(0, 125) + ax_imp.set_ylim(0, 83) + elif cell_name == "porous_42_microns": + ax_imp.set_xlim(0, 12.9) + ax_imp.set_ylim(0, 18.7) + elif cell_name == "porous_80_microns": + ax_imp.set_xlim(0, 12.9) + ax_imp.set_ylim(0, 12.9) + """ + fig_imp.savefig(filename + ".pdf", bbox_inches="tight", pad_inches=0.0) + plt.show() + + +if __name__ == "__main__": + ########################################################## + # Select which cell to parameterize on which data point. # + ########################################################## + + # One of 'monolayer_17_microns', 'porous_42_microns', 'porous_80_microns', + # or '18650_LG_3500_MJ1_EIS_anode_discharge' (recommendation: soc_index 3). + cell_name = "porous_80_microns" + # Choose the parameter prior. + prior_name = "pouch_SPM" + + soc_index = 6 if "microns" in cell_name else 3 + # Select the (de-)lithiation impedance measurement at that point. + # (Not applicable for the delithiation-only LG-MJ1 measurement.) + lithiation = True + + ###################################### + # Read in model parameters and data. # + ###################################### + + electrolyte = False if "SPM" in prior_name else True + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge": + parameters = kuhn.get_parameters_with_switched_electrodes() + sei = True + working_electrode = "both" + three_electrode = "positive" + reference_electrode_location = 0.5 + else: + parameters = gunther.get_parameters_for_cell(cell_name) + sei = False + working_electrode = "positive" + three_electrode = "positive" + reference_electrode_location = 1.0 + if "variant_2" in prior_name: + parameters.update( + { + "Positive electrode double-layer capacity [F.m-2]": 3e-2, + "SEI relative permittivity": 20.0, + } + ) + + parameter_range = parameter_ranges[prior_name] + parameter_transform = parameter_transforms[prior_name] + parameter_transform_numpy = parameter_transforms_numpy[prior_name] + + variable_parameters = list(parameter_range.keys()) + data_index = ( + soc_index + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge" + else soc_index * 2 - 1 + lithiation + ) + electrode = ( + "negative" + if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge" + else "positive" + ) + location = 0.5 if cell_name == "18650_LG_3500_MJ1_EIS_anode_discharge" else 1.0 + frequencies, impedances, soc = preprocess_data(data_index, electrode, cell_name) + parameters[electrode.capitalize() + " electrode SOC"] = soc + + with open("impedance_ocv_alignment.json") as f: + ocv_alignment = json.load(f)[cell_name] + raw_data_index = ocv_alignment["indices"].index(data_index) + with open("drt_finetuning.json") as f: + drt_settings = json.load(f) + lambda_value = drt_settings[cell_name]["lambda"][raw_data_index] + subsampling = drt_settings[cell_name]["subsampling"][raw_data_index] + start_frequency = drt_settings[cell_name]["start_frequency"][raw_data_index] + end_frequency = drt_settings[cell_name]["end_frequency"][raw_data_index] + """ + drt_finetuning = interactive_drt_finetuning( + frequencies.detach().cpu().numpy(), + impedances.detach().cpu().numpy(), + lambda_value, + subsampling, + start_frequency, + end_frequency, + ) + lambda_value = drt_finetuning["lambda"] + subsampling = drt_finetuning["subsampling"] + start_frequency = drt_finetuning["start"] + end_frequency = drt_finetuning["end"] + """ + + start = find_occurrences(frequencies, start_frequency)[0] + end = find_occurrences(frequencies, end_frequency)[0] + optimizer_frequencies = frequencies[start : end + 1 : subsampling] + optimizer_impedances = impedances[start : end + 1 : subsampling] + + unknowns = {key: parameter_range[key] for key in variable_parameters} + transform_unknowns = { + key: (parameter_transform_numpy[key][1], parameter_transform_numpy[key][0]) + for key in variable_parameters + } + + # Make sure that no function parameters try to get passed as numbers. + for key in variable_parameters: + if "function" in str(type(parameters[key])): + forward = transform_unknowns[key][0] + backward = transform_unknowns[key][1] + parameters[key] = forward( + 0.5 * (backward(unknowns[key][0]) + backward(unknowns[key][1])) + ) + + mean = tensor([parameters[key] for key in variable_parameters]) + bounds = tensor( + [ + [parameter_range[key][0] for key in variable_parameters], + [parameter_range[key][1] for key in variable_parameters], + ] + ) + transforms = [parameter_transform[key] for key in variable_parameters] + + ###################################### + # Perform a manual model assessment. # + ###################################### + + """ + manual_model_assessment( + optimizer_frequencies, + optimizer_impedances, + parameters, + unknowns, + transform_unknowns, + working_electrode, + electrolyte, + sei, + three_electrode, + reference_electrode_location, + lambda_value, + ) + """ + + ############################################ + # Perform the model assessment with SOBER. # + ############################################ + + filename = ( + cell_name + + "_" + + prior_name + + "_soc_" + + f"{soc:.3g}" + + "_" + + ("lithiation" if lithiation else "delithiation") + + ".log" + ) + with open(filename, "w") as f: + with redirect_stdout(f): + sober_wrapper, basq_results = automated_model_assessment( + optimizer_frequencies, + optimizer_impedances, + mean, + bounds, + transforms, + cell_name, + soc, + electrolyte, + variable_parameters, + lambda_value, + ) + # from pandas import DataFrame + # from seaborn import pairplot + # pairplot(DataFrame(sober_wrapper.tm.numpy(sober_wrapper.X_all), columns=sober_wrapper.names)) + # plt.show() + visualize_automated_assessment( + sober_wrapper, + basq_results, + frequencies, + cell_name, + prior_name, + soc, + lithiation, + electrolyte, + variable_parameters, + ) diff --git a/examples/scripts/model_selection/visualize_large_variance_evidences.py b/examples/scripts/model_selection/visualize_large_variance_evidences.py new file mode 100644 index 000000000..468f552ce --- /dev/null +++ b/examples/scripts/model_selection/visualize_large_variance_evidences.py @@ -0,0 +1,90 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.stats import lognorm + +""" +Suppose a normally distributed random variable ``X`` has mean ``mu`` +and standard deviation ``sigma``. Then ``Y = exp(X)`` is lognormally +distributed with ``s = sigma`` and ``scale = exp(mu)``. +""" + + +def calculate_distribution_parameters(log_evidence, log_variance): + sigma = np.log(1 + np.exp(log_variance) / np.exp(log_evidence) ** 2) ** 0.5 + mu = log_evidence - sigma**2 / 2 + mean = np.exp(log_evidence) + lower = mean / np.exp(sigma) + upper = mean * np.exp(sigma) + print(mean, "within one confidence interval: [", lower, ",", upper, "]") + return mu, sigma + + +def plot_evidence_distribution(log_evidences, log_variances, labels, title): + fig, ax = plt.subplots(figsize=(3**0.5, 3)) + log_evidences = np.array(log_evidences) + log_variances = np.array(log_variances) + min_extent = float("inf") + max_extent = -float("inf") + for le, lv, label in zip(log_evidences, log_variances, labels, strict=False): + mu, sigma = calculate_distribution_parameters(le, lv) + lower = mu - 2 * sigma + upper = mu + 2 * sigma + evidence_eval = np.logspace(lower, upper, 200, base=np.exp(1)) + min_extent = min(lower, min_extent) + max_extent = max(upper, max_extent) + distribution = lognorm(s=sigma, scale=np.exp(mu)) + distribution_eval = distribution.pdf(evidence_eval) + ax.plot(evidence_eval, distribution_eval, label=label) + ax.set_xlim(0.9 * np.exp(min_extent), 1.1 * np.exp(max_extent)) + ax.legend() + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_title(title) + fig.tight_layout() + plt.show() + + +# SPM impedance results. +""" +plot_evidence_distribution( + [-4.86830e+00, -5.14230e+00, -5.45996e+00], + [-3.39258e+00, -9.72431e+00, -6.52256e+00], + ["17 μm", "42 μm", "80 μm"], + "SPM impedance\n89 % delithiation" +) +""" +plot_evidence_distribution( + [-7.06211e-01, -4.96983e00, -3.85247e00], + [-1.56320e00, -5.19572e00, -4.24408e00], + ["17 μm", "42 μm", "80 μm"], + "SPM impedance\n45 % lithiation", +) + +# DFN impedance results. +""" +plot_evidence_distribution( + [-8.30258e+00, -7.74164e+00, -6.14714e+00], + [-1.06599e+01, -1.76625e+01, -8.05868e+00], + ["17 μm", "42 μm", "80 μm"], + "DFN impedance\n89 % delithiation" +) +""" + +plot_evidence_distribution( + [-1.55120e00, -8.44433e00, -6.67521e00], + [-1.04638e01, -8.70146e00, -8.58370e00], + ["17 μm", "42 μm", "80 μm"], + "DFN impedance\n45 % lithiation", +) + +# SEI impedance results. +plot_evidence_distribution( + [-2.51683e00, -9.07830e00], + [-5.38383e00, -9.82745e00], + ["τ_SEI > τ_DL", "τ_SEI < τ_DL"], + "SEI model impedance", +) + +# Knee point model selection. +calculate_distribution_parameters(-3.00752, -15.8633) +calculate_distribution_parameters(-1.66418, -19.5028) diff --git a/examples/scripts/surrogate_modelling/electrolyte_inverse_model_training.py b/examples/scripts/surrogate_modelling/electrolyte_inverse_model_training.py new file mode 100644 index 000000000..ab47754d1 --- /dev/null +++ b/examples/scripts/surrogate_modelling/electrolyte_inverse_model_training.py @@ -0,0 +1,294 @@ +from copy import deepcopy +from multiprocessing import Pool + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pybamm +import sober +import torch +from ep_bolfi.models.solversetup import solver_setup, spectral_mesh_pts_and_method +from ep_bolfi.utility.preprocessing import calculate_desired_voltage +from scipy.stats import norm +from seaborn import kdeplot, pairplot +from sklearn.cluster import KMeans +from sklearn.metrics import silhouette_score +from sober import InverseModel + +import pybop + +# torch.set_default_dtype(torch.float64) +# sober.setting_parameters(device=torch.device('cpu')) + + +torch.set_default_dtype(torch.float32) +sober.setting_parameters(dtype=torch.float32) # , device=torch.device('cpu')) +torch.set_default_device(sober._settings._device) + +seed = 0 +model = pybamm.lithium_ion.DFN() +noise_generator = norm(0, 1e-4) + + +def simulator(parameters): + global model, noise_generator + + tdf_curvature = parameters[0] + D_e_scaling = parameters[1] + D_e_magnitude = parameters[2] + kappa_e_magnitude = parameters[3] + # kappa_e_peak = parameters[4] + # kappa_e_spread = parameters[5] + kappa_e_peak = 1.0 + kappa_e_spread = 1.0 + + def thermodynamic_factor(c_e, T): + return 1 + tdf_curvature * (c_e / 1000) ** 2 + + def electrolyte_diff(c_e, T): + return D_e_magnitude * pybamm.exp(-D_e_scaling * (c_e / 1000 - 1)) + + def electrolyte_cond(c_e, T): + return ( + kappa_e_magnitude + / (c_e / 1000) + * pybamm.exp( + -((pybamm.log(c_e / 1000 / kappa_e_peak) / kappa_e_spread) ** 2) + ) + ) + + model_parameters = pybamm.ParameterValues("Ramadass2004") + model_parameters["Current function [A]"] = 5.0 + model_parameters["Thermodynamic factor"] = thermodynamic_factor + model_parameters["Electrolyte diffusivity [m2.s-1]"] = electrolyte_diff + model_parameters["Electrolyte conductivity [S.m-1]"] = electrolyte_cond + discretization = { + "order_s_n": 10, + "order_s_p": 10, + "order_e": 10, + "volumes_e_n": 1, + "volumes_e_s": 1, + "volumes_e_p": 1, + "halfcell": False, + } + solver = solver_setup( + deepcopy(model), + model_parameters, + *spectral_mesh_pts_and_method(**discretization), + verbose=False, + ) + relaxation_t = np.logspace(-3, 2, 8) + solution = solver(t_eval=relaxation_t) + # relaxation_t = solution["Time [s]"].entries + relaxation_U = calculate_desired_voltage( + solution, relaxation_t, 1e-3, overpotential=True + ) + relaxation_U += noise_generator.rvs(size=len(relaxation_U)) + + return relaxation_t, relaxation_U, solution + + +def parallel_simulator(parameters): + return simulator(parameters)[1] + + +def training_simulator(parameters): + parameters = parameters.detach().cpu().numpy() + with Pool() as p: + results = p.map(parallel_simulator, parameters) + return torch.tensor(results).to(dtype=torch.float32) + + +def clustering(data, max_number_of_clusters=10): + clusterings = [] + labellings = [] + scores = [] + for i in range(2, max_number_of_clusters + 1): + clusterings.append(KMeans(n_clusters=i)) + labellings.append(clusterings[-1].fit_predict(data)) + scores.append(silhouette_score(data, labellings[-1])) + return clusterings, labellings, scores + + +if __name__ == "__main__": + # Ramadass2004 values: D_e = 7.5e-10, kappa_e = 1.0, TDF = 1.0. + # _, _, solution = simulator([1.0, 0.7, 7.5e-10, 1.0, 1.0, 1.0]) + bounds = torch.tensor( + [ + [0.1, 0.2, 3e-10, 0.5], # , 0.7, 0.7], + [4.0, 2, 1e-9, 2.0], # , 1.3, 1.5], + ] + ) + transforms = [(lambda x: torch.log(x), lambda x: torch.exp(x))] * 4 # 6 + names = [ + "TDF curv.", + "Dₑ scaling", + "Dₑ magn.", + "κₑ magn.", + ] # "κₑ peak", "κₑ spread", + + inverse_modelling = InverseModel( + training_simulator, + model_initial_samples=2**5, + bounds=bounds, + prior="Uniform", + transforms=transforms, + seed=seed, + disable_numpy_mode=True, + parallelization=False, + visualizations=False, + names=names, + ) + inverse_modelling.optimize_inverse_model_with_SOBER( + stopping_criterion_variance=1e-8, + maximum_number_of_batches=7, + sober_iterations_per_convergence_check=7, + sober_iterations_per_training_data_updates=1, + model_samples_per_iteration=2**5, + integration_nodes=2**5, + visualizations=False, + verbose=True, + ) + + relaxation_t, relaxation_U, solution = simulator( + [1.0, 0.7, 7.5e-10, 1.0] # , 1.0, 1.0] + ) + features = training_simulator( + torch.tensor([[1.0, 0.7, 7.5e-10, 1.0]]) # , 1.0, 1.0]]) + ) + + mean, _, (lower_bounds, upper_bounds) = inverse_modelling.evaluate( + features, one_dimensional_confidence=True + ) + + print("Prediction:", mean) + print("Lower bounds:", lower_bounds) + print("Upper bounds:", upper_bounds) + + """ + _, _, sol_lower = simulator(lower_bounds[0].numpy()) + _, _, sol_upper = simulator(upper_bounds[0].numpy()) + t_eval = torch.linspace( + 0, + min([ + sol_lower["Time [s]"].entries[-1], + sol_upper["Time [s]"].entries[-1] + ]), + 101 + ).numpy() + """ + + samples = [ + s[0] for s in inverse_modelling.sample(features, 64).detach().cpu().numpy() + ] + with Pool() as p: + simulations = p.map(simulator, samples) + + fig, ax = plt.subplots(figsize=(4 * 2**0.5, 4)) + ax.plot(relaxation_t / 3600, relaxation_U) + """ + ax.fill_between( + t_eval / 3600, + sol_lower["Voltage [V]"](t_eval), + sol_upper["Voltage [V]"](t_eval), + alpha=0.3, + color='grey' + ) + """ + for sim in simulations: + _, _, sol = sim + ax.plot( + sol["Time [h]"].entries, + calculate_desired_voltage( + sol, sol["Time [s]"].entries, 1e-3, overpotential=True + ), + alpha=0.5, + lw=0.5, + color="orange", + ) + ax.set_xlabel("Time / h") + ax.set_ylabel("Cell overpotential / mV") + ax.set_xscale("log") + + # Get the whole multivariate normal distribution for correlations. + prediction = inverse_modelling(features) + covariance = prediction.covariance_matrix + std = prediction.variance[0].sqrt() + correlation = (covariance / std[:, None]) / std[None, :] + fig_corr, ax_corr = plt.subplots(figsize=(3.75, 3)) + pybop.plot.correlation( + fig_corr, + ax_corr, + correlation.detach().cpu().numpy(), + names, + "Electrolyte characterization from pulse test", + entry_color="black", + ) + + np.set_printoptions(precision=3) + # Cluster the optimal evaulation points for analysis. + normalized_X = inverse_modelling.tm.numpy(inverse_modelling.X_all) + observations = inverse_modelling.tm.numpy(inverse_modelling.observations_all) + clusterings, labellings, scores = clustering(normalized_X, max_number_of_clusters=4) + for i, (clustering, labels, score) in enumerate( + zip(clusterings, labellings, scores, strict=False) + ): + print() + print("Silhouette score of #" + str(i + 2) + " clusters:", score) + labelled_X = ( + pd.DataFrame(normalized_X) + .set_axis(names, axis="columns") + .assign(labels=labels) + ) + correlation_plot = pairplot(labelled_X, height=1.25, hue="labels") + correlation_plot.map_lower(kdeplot, levels=4) + centers = clustering.cluster_centers_ + for j, center in enumerate(centers): + # Find nearest observations. + # index = argmin(sum((normalized_X - center)**2, axis=1)) + # X_index = normalized_X[index] + # observation = observations[index] + # Create new observations. + features = training_simulator(torch.tensor(center).unsqueeze(0)) + prediction = inverse_modelling(features) + model_center = ( + inverse_modelling.reverse_transform( + inverse_modelling.denormalize_input( + torch.tensor( + center, dtype=inverse_modelling.tm.dtype + ).unsqueeze(0) + ) + )[0] + .detach() + .cpu() + .numpy() + ) + model_prediction = ( + inverse_modelling.reverse_transform( + inverse_modelling.denormalize_input(prediction.mean) + )[0] + .detach() + .cpu() + .numpy() + ) + print( + "Parameters ", + model_center, + " - Prediction ", + model_prediction, + " = ", + model_center - model_prediction, + ) + covariance = prediction.covariance_matrix + std = prediction.variance[0].sqrt() + correlation = (covariance / std[:, None]) / std[None, :] + fig_corr, ax_corr = plt.subplots(figsize=(3.75, 3)) + pybop.plot.correlation( + fig_corr, + ax_corr, + correlation.detach().cpu().numpy(), + names, + str(i + 2) + " clusters, label " + str(j) + ", " + str(center), + entry_color="black", + ) + plt.show() diff --git a/examples/scripts/surrogate_modelling/global_inverse_modelling.py b/examples/scripts/surrogate_modelling/global_inverse_modelling.py new file mode 100644 index 000000000..7de412193 --- /dev/null +++ b/examples/scripts/surrogate_modelling/global_inverse_modelling.py @@ -0,0 +1,174 @@ +from copy import deepcopy +from itertools import cycle +from multiprocessing import Pool + +import matplotlib.pyplot as plt +import pybamm +import sober +import torch +from ep_bolfi.models.solversetup import simulation_setup, spectral_mesh_pts_and_method +from ep_bolfi.utility.fitting_functions import fit_sqrt +from scipy.stats import norm +from sober import InverseModel + +torch.set_default_dtype(torch.float64) +sober.setting_parameters(device=torch.device("cpu")) + +seed = 0 +model = pybamm.lithium_ion.DFN() +rest_duration = 300 +rest_fraction_used = 0.1 +period = 0.1 +noise_generator = norm(0, 1e-4) + + +def simulator(parameters): + global model, rest_duration, rest_fraction_used, period, noise_generator + + pulse_strength = parameters[0] + pulse_length = parameters[1] + model_parameters = model.default_parameter_values + + procedure = [ + "Discharge at " + + str(pulse_strength) + + " C for " + + str(pulse_length) + + " seconds (" + + str(period) + + " second period)", + "Rest for " + str(rest_duration) + " seconds (1 second period)", + ] + discretization = { + "order_s_n": 10, + "order_s_p": 10, + "order_e": 10, + "volumes_e_n": 1, + "volumes_e_s": 1, + "volumes_e_p": 1, + "halfcell": False, + } + solver, _ = simulation_setup( + deepcopy(model), + procedure, + model_parameters, + *spectral_mesh_pts_and_method(**discretization), + verbose=False, + ) + solution = solver(calc_esoh=False) + pulse_end = int(pulse_length / period) + 1 + relaxation_t = solution["Time [s]"].entries[ + pulse_end : pulse_end + int(rest_fraction_used * rest_duration) + ] + relaxation_U = solution["Voltage [V]"].entries[ + pulse_end : pulse_end + int(rest_fraction_used * rest_duration) + ] + relaxation_U += noise_generator.rvs(size=len(relaxation_U)) + + return relaxation_t, relaxation_U, solution + + +def training_simulator(parameters): + parameters = parameters.detach().cpu().numpy() + relaxation_t, relaxation_U, _ = simulator(parameters) + sqrt_features = fit_sqrt(relaxation_t, relaxation_U)[2] + sqrt_features = torch.tensor(sqrt_features) + sqrt_features[1] = torch.log(sqrt_features[1]) + return sqrt_features + + +if __name__ == "__main__": + bounds = torch.tensor( + [ + [0.02, 10.0], + [1.0, 600.0], + ] + ) + transforms = [ + (lambda x: torch.log(x), lambda x: torch.exp(x)), + (lambda x: torch.log(x), lambda x: torch.exp(x)), + ] + + inverse_modelling = InverseModel( + training_simulator, + model_initial_samples=128, + bounds=bounds, + prior="Uniform", + transforms=transforms, + seed=seed, + disable_numpy_mode=True, + visualizations=False, + names=["Pulse strength [C]", "Pulse length [s]"], + ) + inverse_modelling.optimize_inverse_model_with_SOBER( + stopping_criterion_variance=1e-12, + maximum_number_of_batches=3, + model_samples_per_iteration=128, + integration_nodes=100, + visualizations=False, + verbose=True, + ) + + relaxation_t, relaxation_U, solution = simulator([0.06, 80.0]) + features = training_simulator(torch.tensor([0.06, 80.0])) + + mean, _, (lower_bounds, upper_bounds) = inverse_modelling.evaluate( + features, one_dimensional_confidence=True + ) + + print("Prediction:", mean) + print("Lower bounds:", lower_bounds) + print("Upper bounds:", upper_bounds) + + """ + _, _, sol_lower = simulator(lower_bounds[0].numpy()) + _, _, sol_upper = simulator(upper_bounds[0].numpy()) + t_eval = torch.linspace( + 0, + min([ + sol_lower["Time [s]"].entries[-1], + sol_upper["Time [s]"].entries[-1] + ]), + 101 + ).numpy() + """ + + samples = [s[0] for s in inverse_modelling.sample(features, 32).numpy()] + with Pool() as p: + simulations = p.map(simulator, samples) + + fig, ax = plt.subplots(figsize=(3 * 2**0.5, 3)) + ax.plot( + (relaxation_t - relaxation_t[0]) / 3600, + relaxation_U, + label="observed portion of the data", + ) + """ + ax.fill_between( + t_eval / 3600, + sol_lower["Voltage [V]"](t_eval), + sol_upper["Voltage [V]"](t_eval), + alpha=0.3, + color='grey' + ) + """ + color_cycle = cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) + for simulation, sample, color in zip( + simulations, samples, color_cycle, strict=False + ): + _, _, solution = simulation + t = solution["Time [h]"].entries + U = solution["Voltage [V]"].entries + U += noise_generator.rvs(size=len(U)) + pulse_length = sample[1] + pulse_end = int(pulse_length / period) + 1 + t_pulse = t[:pulse_end] - t[pulse_end] + U_pulse = U[:pulse_end] + t_rest = t[pulse_end + int(rest_fraction_used * rest_duration) :] - t[pulse_end] + U_rest = U[pulse_end + int(rest_fraction_used * rest_duration) :] + ax.plot(t_pulse, U_pulse, alpha=0.5, lw=0.5, color=color) + ax.plot(t_rest, U_rest, alpha=0.5, lw=0.5, color=color) + ax.set_xlabel("Time / h") + ax.set_ylabel("Cell voltage / V") + ax.legend() + plt.show() diff --git a/pybop/__init__.py b/pybop/__init__.py index 5513ced29..8d2dc6215 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -30,7 +30,7 @@ # # Dataset class # -from ._dataset import Dataset, import_pyprobe_result +from ._dataset import Dataset, Datasets, import_pyprobe_result # # Transformation classes diff --git a/pybop/_dataset.py b/pybop/_dataset.py index 7cdf41953..66dd0d357 100644 --- a/pybop/_dataset.py +++ b/pybop/_dataset.py @@ -179,6 +179,44 @@ def get_subset(self, index: list | np.ndarray): return Dataset(data, domain=self.domain) +class Datasets: + """ + Represents a collection of experimental observations. + """ + + def __init__( + self, datasets: list, domain="Time [s]", control_variable="Current function [A]" + ): + self.datasets = [] + self.domain = domain + self.control_variable = control_variable + for dataset in datasets: + if not isinstance(dataset, Dataset): + dataset = Dataset(data_dictionary=dataset, domain=domain) + self.datasets.append(dataset) + + def get_subset(self, indices): + return Datasets( + [self.datasets[i] for i in indices], self.domain, self.control_variable + ) + + def __iter__(self): + self.count = -1 + return self + + def __next__(self): + self.count += 1 + if self.count >= len(self): + raise StopIteration + return self.datasets[self.count] + + def __len__(self): + return len(self.datasets) + + def __getitem__(self, i): + return self.datasets[i] + + def import_pyprobe_result( result: PyprobeResult, pybop_columns: list[str] | None = None, diff --git a/pybop/costs/endoflife_cost.py b/pybop/costs/endoflife_cost.py new file mode 100644 index 000000000..27091844d --- /dev/null +++ b/pybop/costs/endoflife_cost.py @@ -0,0 +1,60 @@ +from pybop._utils import FailedSolution +from pybop.costs.base_cost import BaseCost +from pybop.parameters.parameter import Inputs +from pybop.simulators.base_simulator import Solution + + +class EndOfLifeCost(BaseCost): + """ + Cost for optimising the lifetime of a battery modelled with PyBaMM. + + End-Of-Life is usually defined as the point when the battery has only + a certain fraction of its initial capacity left. + """ + + def __init__(self, relative_capacity_cutoff): + super().__init__() + self.minimising = False + self.domain = "Time [s]" + self.relative_capacity_cutoff = relative_capacity_cutoff + + def evaluate( + self, + sol: Solution | FailedSolution, + inputs: Inputs | None = None, + calculate_sensitivities: bool = False, + ) -> float: + if isinstance(sol, FailedSolution): + return self.failure(calculate_sensitivities) + eol_index = indices_of( + (sol["SEI thicknesses [m]"].entries - sol["SEI thicknesses [m]"].entries(0)) + * sol["Negative electrode surface area to volume ratio [m-1]"] + * F + / (3600 * sol["Nominal cell capacity [A.h]"]), + 1 - self.relative_capacity_cutoff, + )[0] + return sol["Time [s]"].entries[eol_index] + + +class RelativeEndOfLifeCost(EndOfLifeCost): + """ + Cost for optimising the lifetime gained from battery adjustments. + + Requires the battery to be modelled with PyBaMM and the extra adjustment + encoded as a pybamm.Parameter("Adjustment factor") relative to 1. + """ + + def __init__(self, relative_capacity_cutoff, eol_reference): + super().__init__(relative_capacity_cutoff) + self.eol_reference = eol_reference + + def evaluate( + self, + sol: Solution | FailedSolution, + inputs: Inputs | None = None, + calculate_sensitivities: bool = False, + ) -> float: + if isinstance(sol, FailedSolution): + return self.failure(calculate_sensitivities) + eol_time = super().evaluate(sol) + return (eol_time - eol_reference) / (sol["Adjustment factor"] - 1) diff --git a/pybop/models/_diffusive_decay.py b/pybop/models/_diffusive_decay.py new file mode 100644 index 000000000..4d4116137 --- /dev/null +++ b/pybop/models/_diffusive_decay.py @@ -0,0 +1,61 @@ +from numpy import pi, exp, sin, cos + + +class Diffusive_Relaxation: + """Solution to ∂ₜ u = D ∂ₓ² u with u(x, t=0) = f(x).""" + + def __init__(self, f, L, summands=10, radial=False): + self.f = f + self.L = L + self.summands = summands + self.radial = radial + if radial: + pass + else: + self.series = cos + self.coefficients = self.compute_zero_flow_coefficients() + + def compute_zero_flow_coefficients(self): + coefficients = tensor( + [ + 2.0 + / self.L + * quad(lambda x: self.f(x) * cos(n * pi * x / self.L), 0, self.L)[ + 0 + ] + for n in range(0, self.summands) + ] + ) + # In order to use a simple summation expression suitable for + # automatic differentiation, half the "zeroth" coefficient. + coefficients[0] = coefficients[0] / 2.0 + return coefficients + + def compute_zero_concentration_coefficients(self): + # In order to use the same summation expression for both cases, + # the "zeroth" coefficient is set to 0. + return tensor([0] + [ + 2.0 / self.L * quad( + lambda x: self.f(x) * sin(n * pi * x / self.L), 0, self.L)[0] + for n in range(1, self.summands) + ]) + + def concentration(self, x, t, D=1.0): + value = zeros_like(D * t) + for n in range(self.summands): + value += ( + self.coefficients[n] + * self.series(n * pi * x / self.L) + * exp(-(n**2) * pi**2 * D * t / self.L**2) + ) + return value + + def __call__(self, t, offset=0.0, timescale=1.0, magnitude=1.0): + D = self.L**2 / timescale + return offset + magnitude * ( + self.concentration(self.L, t, D) + - self.concentration(self.L, t, D) + - self.concentration(1.0, 0.0, D) + + self.concentration(0.0, 0.0, D) + ) + diff --git a/pybop/models/kneepoints.py b/pybop/models/kneepoints.py new file mode 100644 index 000000000..ddffcf7e9 --- /dev/null +++ b/pybop/models/kneepoints.py @@ -0,0 +1,63 @@ +import numpy as np +import pybop +from torch import tensor + +class KneepointModel(pybop.BaseSimulator): + def __init__(self, parameters, t, n_kneepoints=2): + super().__init__(parameters) + self.t = t + self.output_variables = ["Capacity fade"] + self.n_kneepoints = n_kneepoints + + def one_kneepoint_model(self, parameters): + first_slope = parameters[0].reshape(-1, 1) + kneepoint = parameters[1].reshape(-1, 1) + second_slope = parameters[2].reshape(-1, 1) + + return ( + (1.0 - first_slope * self.t) * (self.t < kneepoint) + + ((1.0 - first_slope * kneepoint) - second_slope * (self.t - kneepoint)) + * (self.t >= kneepoint) + ).T + + def two_kneepoints_model(self, parameters): + first_slope = parameters[0].reshape(-1, 1) + first_kneepoint = parameters[1].reshape(-1, 1) + second_slope = parameters[2].reshape(-1, 1) + second_kneepoint = parameters[3].reshape(-1, 1) + first_kneepoint + third_slope = parameters[4].reshape(-1, 1) + + return ( + (1.0 - first_slope * self.t) * (self.t < first_kneepoint) + + ( + (1.0 - first_slope * first_kneepoint) + - second_slope * (self.t - first_kneepoint) + ) + * (self.t >= first_kneepoint) + * (self.t < second_kneepoint) + + ( + (1.0 - first_slope * first_kneepoint) + - second_slope * (second_kneepoint - first_kneepoint) + - third_slope * (self.t - second_kneepoint) + ) + * (self.t >= second_kneepoint) + ).T + + def solve_batch(self, inputs, calculate_sensitivities=False): + inputs_array = tensor(np.asarray([entry for entry in inputs[0].values()])) + capacity_fade = self(inputs_array) + sols = [] + for entry, cf in zip(inputs, capacity_fade, strict=False): + sol = pybop.Solution(entry) + sol.set_solution_variable("Capacity fade", cf) + sols.append(sol) + return sols + + def __call__(self, parameters): + if self.n_kneepoints == 1: + return self.one_kneepoint_model(parameters) + elif self.n_kneepoints == 2: + return self.two_kneepoints_model(parameters) + else: + raise ValueError("Only one or two kneepoints are implemented.") + diff --git a/pybop/models/lithium_ion/silicon_relaxation.py b/pybop/models/lithium_ion/silicon_relaxation.py new file mode 100644 index 000000000..1d3e2bb8d --- /dev/null +++ b/pybop/models/lithium_ion/silicon_relaxation.py @@ -0,0 +1,144 @@ +from pybamm import citations, constants, Parameter, ParameterValues, exp, log, tanh, DummySolver +from pybamm import t as pybamm_t +from pybamm.models.base_model import BaseModel + +from pybop.models._diffusive_decay import Diffusive_Relaxation + + +class SiliconRelaxation(BaseModel): + """ + Calculates the mechanical voltage response from silicon at rest. + + This model from Koebbing et al. is a core-shell (silicon-SEI/oxide) model. + """ + + def __init__(self, name="Silicon voltage model after Koebbing2024", **kwargs): + super().__init__(name=name, **kwargs) + self._summary_variables = [] + + citations.register("""@article{ + Köbbing2024, + title={{Slow Voltage Relaxation of Silicon Nanoparticles with a Chemo-Mechanical Core-Shell Model}}, + author={Köbbing, L and Kuhn, Y and Horstmann, B}, + journal={ACS Applied Materials & Interfaces}, + volume={16}, + pages={67609-67619}, + year={2024} + }""") + + ############## + # Parameters # + ############## + + # Parameters that match onto existing PyBaMM notation. + R_core = Parameter("Negative particle radius [m]") + D = Parameter("Negative particle diffusivity [m]") + L_shell = Parameter("Initial SEI thickness [m]") + T = Parameter("Ambient temperature [K]") + F = constants.F + R = constants.R + # New parameters for the mechanical properties. + E_core = Parameter("Negative particle Young modulus [Pa]") + nu_core = Parameter("Negative particle Poisson ratio") + sigma_Y_core = Parameter("Negative particle yield stress [Pa]") + v_Li = Parameter("Negative particle partial molar volume [m3.mol-1]") + E_shell = Parameter("Negative shell Young modulus [Pa]") + nu_shell = Parameter("Negative shell Poisson ratio") + sigma_Y_shell = Parameter("Negative shell yield stress [Pa]") + eta_shell = Parameter("Negative shell Newtonian viscosity [Pa.s]") + + G_core = E_core / (2 * (1 + nu_core)) # Second Lame constant. + lambda_core = 2 * G_core * nu_core / (2 * (1 + nu_core)) # First Lame constant. + + alpha = 0.5 * (R_core / L_shell - 1) + + U_infty = Parameter("Terminal voltage of observed relaxation [V]") + slope = Parameter("Logarithmic slope of observed mechanical relaxation [V]") + timescale_exp = Parameter("Exponential timescale for decay of observed mechanical relaxation [s]") + + # Express the effective parameters from above by adjusting the parameters the model is written in. + lambda_ch = 1 # lumped parameter of the model, set to 1 without loss of generality + sigma_ref = slope * alpha * F * lambda_ch**3 / (2 * v_Li) + tau = lambda_ch * E_core * alpha * timescale_exp / sigma_ref + + # Optionally add a diffusive relaxation. + if kwargs.get("with diffusion"): + diff_portion = Parameter("Fraction of observed relaxation that is diffusive") + diff_timescale = Parameter("Timescale of observed diffusive relaxation [s]") + rel_magnitude = (1 - diff_portion) * U_infty + diff_magnitude = diff_portion * U_infty + else: + rel_magnitude = U_infty + diff_magnitude = 0 + + # Determine integration constant from t=0 with + # sigma_ev(t=0) = sigma_0 from sigma_ev = delta_U * F / v_Li. + integration_constant = tanh( + rel_magnitude * alpha * F * lambda_ch**3 / (2 * v_Li * sigma_ref) + ) + + # As PyBaMM has no arctanh, write it via arctanh = 0.5 * log((1 + x) / (1 - x)). + arctanh_argument = integration_constant * exp(-E_core * alpha * lambda_ch / (tau * sigma_ref) * pybamm_t) + delta_U = ( + 2 * v_Li * sigma_ref / (alpha * F * lambda_ch**3) + * 0.5 * log((1 + arctanh_argument) / (1 - arctanh_argument)) + ) + + if kwargs.get("with diffusion"): + delta_U += Diffusive_Relaxation(lambda x: x, L=1)( + pybamm_t, -diff_magnitude, diff_timescale, diff_magnitude + ) + + self.variables = { + "Voltage [V]": delta_U - U_infty, + "Time [s]": pybamm_t, + "Current [A]": 0, + "Current variable [A]": 0, + } + + @property + def default_geometry(self): + return {} + + @property + def default_parameter_values(self) -> ParameterValues: + return ParameterValues({ + "Negative particle radius [m]": 30e-9, + "Negative particle diffusivity [m]": 1e-17, + "Initial SEI thickness [m]": 20e-9, + "Ambient temperature [K]": 298, + "Negative particle Young modulus [Pa]": 200e9, + "Negative particle Poisson ratio": 0.22, + "Negative particle yield stress [Pa]": 3e9, + "Negative particle partial molar volume [m3.mol-1]": 9e-6, + "Negative shell Young modulus [Pa]": 100e9, + "Negative shell Poisson ratio": 0.3, + "Negative shell yield stress [Pa]": 2e9, + "Negative shell Newtonian viscosity [Pa.s]": 135e12, + "Terminal voltage of observed relaxation [V]": 0.12, + "Logarithmic slope of observed mechanical relaxation [V]": 0.025, + "Exponential timescale for decay of observed mechanical relaxation [s]": 2e4, + "Fraction of observed relaxation that is diffusive": 0.0, + "Timescale of observed diffusive relaxation [s]": 90.0, + }) + + @property + def default_quick_plot_variables(self): + return ["Voltage [V]",] + + @property + def default_submesh_types(self): + return {} + + @property + def default_var_pts(self): + return {} + + @property + def default_spatial_methods(self): + return {} + + @property + def default_solver(self): + return DummySolver() + diff --git a/pybop/optimisers/ep_bolfi_optimiser.py b/pybop/optimisers/ep_bolfi_optimiser.py index c428d5ee0..77a8e2f01 100644 --- a/pybop/optimisers/ep_bolfi_optimiser.py +++ b/pybop/optimisers/ep_bolfi_optimiser.py @@ -203,7 +203,8 @@ def __init__( # author={Minka, T}, # journal={Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence (UAI2001)}, # pages={362-369}, - # year={2013} + # year={2013}, + # doi={10.48550/arXiv.1301.2294} # }""") citations.register("""@article{ Barthelme2014, @@ -212,7 +213,8 @@ def __init__( journal={Journal of the American Statistical Association}, volume={109}, pages={315-333}, - year={2014} + year={2014}, + doi={10.1080/01621459.2013.864178} }""") citations.register("""@article{ Gutmann2016, @@ -221,7 +223,8 @@ def __init__( journal={Journal of Machine Learning Research}, volume={17}, pages={1-47}, - year={2016} + year={2016}, + doi={arXiv.1501.03291} }""") citations.register("""@article{ Kuhn2022, @@ -231,7 +234,8 @@ def __init__( volume={6}, pages={e202200374}, year={2023}, - publisher={Chemistry Europe} + publisher={Chemistry Europe}, + doi={10.1002/batt.202200374} }""") def _set_up_optimiser(self): @@ -371,7 +375,7 @@ def _run(self) -> "BayesianOptimisationResult": for entry in x_best_over_time ] self._logger.x_search_best = x_search_best_over_time[-1] - self._logger.cost_best = cost_best[0] + self._logger.cost_best = cost_best[-1] model_mean_dict = { key: value[0] for key, value in ep_bolfi_result["inferred parameters"].items() diff --git a/pybop/optimisers/sober_basq_optimiser.py b/pybop/optimisers/sober_basq_optimiser.py new file mode 100644 index 000000000..333eb3ead --- /dev/null +++ b/pybop/optimisers/sober_basq_optimiser.py @@ -0,0 +1,460 @@ +import copy +import time +from collections.abc import Callable +from contextlib import redirect_stderr, redirect_stdout +from dataclasses import dataclass +from os import cpu_count +from sys import stderr, stdout + +import numpy as np +import torch +from pybamm import citations +from torch import tensor + +import pybop +from pybop import BaseOptimiser +from pybop._logging import Logger +from pybop.optimisers.ep_bolfi_optimiser import BayesianOptimisationResult + + +@dataclass +class SOBER_BASQ_Options(pybop.OptimiserOptions): + """ + A class to hold SOBER and BASQ options for the optimisation as well + as the optional model selection process. + """ + + model_initial_samples: int = cpu_count() + maximise: bool = False + set_up_parabolic_hyperparameters: bool = False + weights: np.ndarray | None = None + custom_objective_and_loglikelihood: Callable | None = None + seed: int | None = None + batched_input: bool = False + + sober_iterations: int = 1 + model_samples_per_iteration: int = cpu_count() + + integration_nodes: int | None = cpu_count() + + normalise_evidence: bool = True + + verbose: bool = False + + def validate(self): + super().validate() + + +class SOBER_BASQ(BaseOptimiser): + """ + Wraps the Bayesian Optimization algorithm SOBER. SOBER includes the + Bayesian Model Selection algorithm BASQ. + + For implementation details and background information, consult the + relevant publications at https://doi.org/10.48550/arXiv.2404.12219 + and https://proceedings.neurips.cc/paper_files/paper/2022/file/697200c9d1710c2799720b660abd11bb-Paper-Conference.pdf + and visit https://github.com/ma921/SOBER/. + + Note that all properties may and should be given here as PyBOP + objects, but will be converted to sober.SoberWrapper instance + upon instantiation of this class. To change attributes, re-init. + + Only compatible with MulitvariateParameters and multivariate priors. + """ + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + citations.register("""@article{ + Adachi2024, + title={{A Quadrature Approach for General-Purpose Batch Bayesian Optimization via Probabilistic Lifting}}, + author={Adachi, M and Hawakawa, S and Jørgensen, M, and Hamid, S and Oberhauser, H and Osborne, M}, + journal={arXiv}, + year={2024}, + doi={10.48550/arXiv.2404.12219} + }""") + citations.register("""@article{ + Adachi2022, + title={{Fast Bayesian Inference with Batch Bayesian Quadrature via Kernel Recombination}}, + author={Adachi, M and Hayakawa, S and Jørgensen, M and Oberhauser, H and Osborne, M}, + journal={Advances in Neural Information Processing Systems}, + volume={35}, + pages={16533-16547}, + year={2022} + }""") + + def evaluate_problem(self, inputs_array): + sim = self.problem.simulator + inputs = { + k: i + for k, i in zip( + self.problem.parameters.keys(), inputs_array.T, strict=False + ) + } + return sim.solve(inputs)[sim.output_variables[0]].data + + def _set_up_optimiser(self, **kwargs): + import sober + + prior = self.problem.parameters.distribution + + self.mean = prior.properties.get("mean") + if self.mean is not None: + self.mean = torch.tensor(self.mean) + self.covariance = prior.properties.get("cov") + self.bounds = prior.properties.get("bounds") + if self.bounds is None and hasattr(prior, "bounds"): + self.bounds = prior.bounds + + if isinstance(prior, pybop.MultivariateUniform): + self.prior_name = "Uniform" + elif isinstance(prior, pybop.MultivariateGaussian): + self.prior_name = "TruncatedGaussian" + else: + raise ValueError( + "The provided prior must be a multivariate uniform or multivariate Gaussian one." + ) + + # ToDo: generalise to other transformations (has to be PyTorch, + # else the vmap-vectorisation for that within SoberWrapper fails). + self.transform_parameters = [ + (torch.log, torch.exp) + if isinstance(par.transformation, pybop.LogTransformation) + else (torch.nn.Identity(), torch.nn.Identity()) + for par in self.problem.parameters.values() + ] + + # ToDo: can only use one problem function for now, else multiprocessing breaks. + if isinstance(self.problem.target_data, dict): + target_data = tensor(np.asarray(list(self.problem.target_data.values())[0])) + else: + target_data = tensor(self.problem.target_data) + self.optimiser = sober.SoberWrapper( + self.evaluate_problem, + target_data, + self._options.model_initial_samples, + self.mean, + None if self.covariance is None else tensor(self.covariance), + None if self.bounds is None else tensor(self.bounds).T, + self.prior_name, + self._options.maximise, + self._options.set_up_parabolic_hyperparameters, + self._options.weights, + self._options.custom_objective_and_loglikelihood, + self.transform_parameters, + self._options.seed, + False, # disable_numpy_mode + not self._options.batched_input, # parallelization (i.e., False has to parallelise itself) + False, # visualizations, + None, # true_optimum, + True, # standalone + None, # names, + ) + self._logger = Logger( + minimising=not self._options.maximise, + verbose=self._options.verbose, + verbose_print_rate=self.verbose_print_rate, + ) + + def _run(self): + verbose_log_target = stdout if self._options.verbose else None + verbose_err_target = stderr if self._options.verbose else None + with redirect_stdout(verbose_log_target): + with redirect_stderr(verbose_err_target): + start = time.time() + self.optimiser.run_SOBER( + self._options.sober_iterations, + self._options.model_samples_per_iteration, + visualizations=False, + verbose=self._options.verbose, + ) + ( + raw_taken_samples, + MAP, + best_observed, + log_expected_marginal_likelihood, + log_approx_variance_marginal_likelihood, + ) = self.optimiser.run_BASQ( + self._options.integration_nodes, + return_raw_samples=True, + visualizations=False, + verbose=self._options.verbose, + ) + end = time.time() + + x_list = self.optimiser.X_all.numpy() + cost_list = self.optimiser.Y_all.numpy() + x_best = copy.deepcopy(x_list) + cost_best = copy.deepcopy(cost_list) + for i in range(1, len(cost_list)): + if cost_list[i] < cost_best[i - 1]: + x_best[i:, None] = x_list[i, None] + cost_best[i:] = cost_list[i] + self._logger.x_model = x_list.tolist() + self._logger.x_search = [ + [ + par.transformation.to_search(e)[0] + for e, par in zip(entry, self.problem.parameters.values(), strict=False) + ] + for entry in x_list + ] + self._logger.cost = cost_list + self._logger.iterations = [ + i // (self._options.sober_iterations) for i in range(len(cost_list)) + ] + self._logger.evaluations = [i + 1 for i in range(len(cost_list))] + self._logger.x_model_best = x_best[-1] + x_search_best_over_time = [ + [ + par.transformation.to_search(e)[0] + for e, par in zip(entry, self.problem.parameters.values(), strict=False) + ] + for entry in x_best + ] + self._logger.x_search_best = x_search_best_over_time[-1] + self._logger.cost_best = cost_best[-1] + posterior = pybop.MultivariateParameters( + self.problem.parameters, + distribution=pybop.MultivariateNonparametric( + self.optimiser.denormalize_input(raw_taken_samples).T + ), + ) + + self._logger.iteration = {"SOBER iterations": self._options.sober_iterations} + self._logger.evaluations = { + "model evaluations": self._options.sober_iterations + * self._options.model_samples_per_iteration + } + + if self._options.normalise_evidence: + log_expected_marginal_likelihood *= (2 * np.pi) ** ( + len(self.problem.parameters) / 2 + ) + + return BayesianOptimisationResult( + optim=self, + time=end - start, + method_name="SOBER + BASQ", + lower_bounds=None if self.bounds is None else self.bounds[:, 0], + upper_bounds=None if self.bounds is None else self.bounds[:, 1], + posterior=posterior, + maximum_a_posteriori=MAP, + log_evidence_mean=log_expected_marginal_likelihood, + log_evidence_variance=log_approx_variance_marginal_likelihood, + ) + + def name(self): + return "Solving Optimisation as Bayesian Estimation via Recombination" + + +@dataclass +class SOBER_BASQ_EPLFI_Options(SOBER_BASQ_Options): + """ + Extends SOBER and BASQ options with EP-LFI-specific options. + """ + + ep_iterations: int = 2 + ep_total_dampening: float = 0.5 + ep_integration_nodes: int | None = cpu_count() + + def validate(self): + super().validate() + + if self.maximise: + raise ValueError( + "EP-LFI only minimises; consider negating the cost function." + ) + if self.weights is not None: + raise ValueError( + "EP-LFI performs automatic importance weighting; you can't set the weights." + ) + if self.custom_objective_and_loglikelihood is not None: + raise ValueError( + "EP-LFI builds a custom objective from the features already; you can't set your own." + ) + + +class SOBER_BASQ_EPLFI(SOBER_BASQ): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + citations.register("""@article{ + Barthelme2014, + title={{Expectation propagation for likelihood-free inference}}, + author={Barthelmé, S and Chopin, N}, + journal={Journal of the American Statistical Association}, + volume={109}, + pages={315-333}, + year={2014}, + doi={10.1080/01621459.2013.864178} + }""") + + def _collect_functions(self, inputs_array): + inputs = { + k: i + for k, i in zip( + self.problem.parameters.keys(), inputs_array.T, strict=False + ) + } + return np.array( + [ + problem.simulator.solve(inputs)[ + problem.simulator.output_variables[0] + ].data + for problem in self.problem.problems + ] + ).T + + def _features(self, y): + # Since in PyBOP, costs are already part of the problems, just report y. + return y + + def _set_up_optimiser(self, **kwargs): + import sober + + prior = self.problem.parameters.distribution + + self.mean = prior.properties.get("mean") + self.covariance = prior.properties.get("cov") + self.bounds = prior.properties.get("bounds") + if self.bounds is None and hasattr(prior, "bounds"): + self.bounds = prior.bounds + + if isinstance(prior, pybop.MultivariateUniform): + self.prior_name = "Uniform" + elif isinstance(prior, pybop.MultivariateGaussian): + self.prior_name = "TruncatedGaussian" + else: + raise ValueError( + "The provided prior must be a multivariate uniform or multivariate Gaussian one." + ) + + # ToDo: generalise to other transformations (has to be PyTorch, + # else the vmap-vectorisation for that within SoberWrapper fails). + self.transform_parameters = [ + (torch.log, torch.exp) + if isinstance(par.transformation, pybop.LogTransformation) + else (torch.nn.Identity(), torch.nn.Identity()) + for par in self.problem.parameters.values() + ] + + # ToDo: can only use one Python problem function for now, else + # multiprocessing breaks. + self.optimiser = sober.ExpectationPropagationLFI( + self._collect_functions, + tensor( + np.asarray( + [ + list(problem.target_data.values())[0] + for problem in self.problem.problems + ] + ) + ), + self._features, + self._options.model_initial_samples, + self.mean, + self.covariance, + None if self.bounds is None else tensor(self.bounds).T, + self._options.set_up_parabolic_hyperparameters, + self.transform_parameters, + self._options.seed, + False, # disable_numpy_mode + not self._options.batched_input, # parallelization (i.e., False has to parallelise itself) + False, # visualizations, + None, # true_optimum, + ) + self._logger = Logger( + minimising=not self._options.maximise, + verbose=self._options.verbose, + verbose_print_rate=self.verbose_print_rate, + ) + + def _run(self): + verbose_log_target = stdout if self._options.verbose else None + verbose_err_target = stderr if self._options.verbose else None + with redirect_stdout(verbose_log_target): + with redirect_stderr(verbose_err_target): + start = time.time() + self.optimiser.run_Expectation_Propagation( + ep_iterations=self._options.ep_iterations, + final_dampening=self._options.ep_total_dampening, + sober_iterations=self._options.sober_iterations, + model_samples_per_iteration=self._options.model_samples_per_iteration, + integration_nodes=self._options.ep_integration_nodes, + visualizations=False, + verbose=False, + ) + ( + raw_taken_samples, + MAP, + best_observed, + log_expected_marginal_likelihood, + log_approx_variance_marginal_likelihood, + ) = self.optimiser.run_BASQ( + self._options.integration_nodes, + return_raw_samples=True, + visualizations=False, + verbose=False, + ) + end = time.time() + + x_list = self.optimiser.X_all.numpy() + cost_list = self.optimiser.Y_all.numpy() + x_best = copy.deepcopy(x_list) + cost_best = copy.deepcopy(cost_list) + for i in range(1, len(cost_list)): + if cost_list[i] < cost_best[i - 1]: + x_best[i:, None] = x_list[i, None] + cost_best[i:] = cost_list[i] + self._logger.x_model = x_list.tolist() + self._logger.x_search = [ + [ + par.transformation.to_search(e)[0] + for e, par in zip(entry, self.problem.parameters.values(), strict=False) + ] + for entry in x_list + ] + self._logger.cost = cost_list + self._logger.iterations = [ + i // (self._options.sober_iterations) for i in range(len(cost_list)) + ] + self._logger.evaluations = [i + 1 for i in range(len(cost_list))] + self._logger.x_model_best = x_best[-1] + x_search_best_over_time = [ + [ + par.transformation.to_search(e)[0] + for e, par in zip(entry, self.problem.parameters.values(), strict=False) + ] + for entry in x_best + ] + self._logger.x_search_best = x_search_best_over_time[-1] + self._logger.cost_best = cost_best + self._logger.cost_best = cost_best[-1] + posterior = pybop.MultivariateParameters( + self.problem.parameters, + distribution=pybop.MultivariateNonparametric( + self.optimiser.denormalize_input(raw_taken_samples).T + ), + ) + + self._logger.iteration = {"SOBER iterations": self._options.sober_iterations} + self._logger.evaluations = { + "model evaluations": self._options.sober_iterations + * self._options.model_samples_per_iteration + } + + return BayesianOptimisationResult( + optim=self, + time=end - start, + method_name="EP-LFI + SOBER + BASQ", + lower_bounds=None if self.bounds is None else self.bounds[:, 0], + upper_bounds=None if self.bounds is None else self.bounds[:, 1], + posterior=posterior, + maximum_a_posteriori=MAP, + log_evidence_mean=log_expected_marginal_likelihood, + log_evidence_variance=log_approx_variance_marginal_likelihood, + ) + + def name(self): + return ( + "Solving Optimisation as Bayesian Estimation via Recombination " + "within Expectation Propagation for Likelihood-Free Inference" + ) diff --git a/pybop/parameters/multivariate_distributions.py b/pybop/parameters/multivariate_distributions.py index 5161165f4..eef15ac7f 100644 --- a/pybop/parameters/multivariate_distributions.py +++ b/pybop/parameters/multivariate_distributions.py @@ -6,6 +6,25 @@ from pybop.parameters.distributions import Distribution +def insert_x_at_dim(other, x, dim): + """ + Inserts each element of x at dim in other. + + Parameters + ---------- + other : numpy.ndarray + The fixed other variables in all dimensions but dim. + x : numpy.ndarray + The point(s) at which to evaluate in dim. + dim : int + The dimension at which to insert x. + """ + x = np.atleast_1d(x) + return np.array( + [np.concatenate(other[: dim - 1], [entry], other[dim - 1 :]) for entry in x] + ) + + class BaseMultivariateDistribution(Distribution): """ A base class for defining multivariate parameter distributions. @@ -29,9 +48,52 @@ class BaseMultivariateDistribution(Distribution): def pdf(self, x): return self.distribution.pdf(x, **self.properties) + def pdf_marginal(self, x, dim): + """ + Integrates the probability density function (PDF) at x down to + one variable. + + Parameters + ---------- + x : numpy.ndarray + The point(s) at which to evaluate the pdf. + dim : int + The dimension to which to reduce the pdf to. + + Returns + ------- + float + The marginal probability density function value at x in dim. + """ + unnormalized_pdf = integrate.nquad( + lambda y: self.pdf(insert_x_at_dim(y, x, dim)), + [[-np.inf, np.inf] * (self._n_parameters - 1)], + )[0] + return unnormalized_pdf / np.sum(unnormalized_pdf) + def logpdf(self, x): return self.distribution.logpdf(x, **self.properties) + def logpdf_marginal(self, x, dim): + """ + Integrates the logarithm of the probability density function + (PDF) at x down to one variable. + + Parameters + ---------- + x : numpy.ndarray + The point(s) at which to evaluate the PDF. + dim : int + The dimension to which to reduce the PDF to. + + Returns + ------- + float + The log marginal probability density function value at x in + dim. + """ + return np.log(self.pdf_marginal(x, dim)) + icdf = None """Multivariate distributions have no invertible CDF.""" @@ -155,6 +217,7 @@ def __init__(self, mean=None, covariance=None, bounds=None, random_state=None): super().__init__(distribution=stats.multivariate_normal) self.name = "MultivariateGaussian" self.properties = {"mean": mean, "cov": covariance} + self.bounds = bounds self.sigma2 = covariance self._n_parameters = len(mean) diff --git a/pybop/parameters/multivariate_parameters.py b/pybop/parameters/multivariate_parameters.py index e89fc8cf1..0838a0648 100644 --- a/pybop/parameters/multivariate_parameters.py +++ b/pybop/parameters/multivariate_parameters.py @@ -44,8 +44,8 @@ def pdf(self, x: np.ndarray) -> np.ndarray: An array of the pdf values at each parameter vector. """ if len(x.shape) == 1: # one-dimensional - x = np.atleast_2d(x) - x = np.asarray([self.transformation.to_search(m) for m in x]) + x = np.atleast_2d(x).T + x = np.asarray([self.transformation.to_search(m) for m in x]).T return self.distribution.pdf(x) def sample_from_distribution( @@ -77,7 +77,7 @@ def rvs( samples = self.distribution.rvs(n_samples, random_state=random_state) if samples.ndim < 2: samples = np.atleast_2d(samples) - + samples = samples.T if apply_transform: samples = np.asarray([self.transformation.to_model(s) for s in samples]) diff --git a/pybop/plot/__init__.py b/pybop/plot/__init__.py index 47f5c8f06..08d84b79f 100644 --- a/pybop/plot/__init__.py +++ b/pybop/plot/__init__.py @@ -6,6 +6,7 @@ from .contour import contour from .dataset import dataset from .convergence import convergence +from .correlation import correlation from .parameters import parameters from .problem import problem from .nyquist import nyquist diff --git a/pybop/plot/correlation.py b/pybop/plot/correlation.py new file mode 100644 index 000000000..7a7483b16 --- /dev/null +++ b/pybop/plot/correlation.py @@ -0,0 +1,76 @@ +import matplotlib +import matplotlib.pyplot as plt +import numpy as np + + +def correlation( + fig, + ax, + correlation, + names=None, + title=None, + cmap=None, + entry_color="w", +): + """ + Produces a heatmap of a correlation matrix. + + :param fig: + The ``matplotlib.Figure`` object for plotting. + :param ax: + The ``matplotlib.Axes`` object for plotting. + :param correlation: + A two-dimensional (NumPy) array that is the correlation matrix. + :param names: + A list of strings that are names of the variables corresponding + to each row or column in the correlation matrix. + :param title: + The title of the heatmap. + :param cmap: + The matplotlib colormap for the heatmap. + :param entry_color: + The colour of the correlation matrix entries. + """ + if cmap is None: + cmap = plt.get_cmap("BrBG") + + # This one line produces the heatmap. + ax.imshow(correlation, cmap=cmap, norm=matplotlib.colors.Normalize(-1, 1)) + + # Define the coordinates of the ticks. + ax.set_xticks(np.arange(len(correlation))) + ax.set_yticks(np.arange(len(correlation))) + + # Display the names alongside the rows and columns. + if names is not None: + ax.set_xticklabels(names) + ax.set_yticklabels(names) + # Rotate the labels at the x-axis for better readability. + plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") + + # Plot the correlation matrix entries on the heatmap. + for i in range(len(correlation)): + for j in range(len(correlation)): + if i == j: + color = "w" + else: + color = entry_color + ax.text( + j, + i, + f"{correlation[i][j]:3.2f}", + ha="center", + va="center", + color=color, + in_layout=False, + ) + + ax.set_title(title or "Correlation matrix") + fig.colorbar( + matplotlib.cm.ScalarMappable( + norm=matplotlib.colors.Normalize(-1, 1), cmap=cmap + ), + ax=ax, + label="correlation", + ) + fig.tight_layout() diff --git a/pybop/plot/predictive.py b/pybop/plot/predictive.py new file mode 100644 index 000000000..cbee123a7 --- /dev/null +++ b/pybop/plot/predictive.py @@ -0,0 +1,73 @@ +import numpy as np +import plotly.express as px + +from pybop.plot.standard_plots import StandardPlot + + +def predictive( + result: "BayesianOptimisationResult", + simulator, + number_of_traces=8, + dataset_y="Voltage [V]", + data_legend_entry=None, + rvs_legend_entry=None, + pdf_plot=None, + pdf_label="PDF", + colour_scale='viridis', + show=True, + **layout_kwargs +): + """ + Plot the predictive posterior of a Bayesian parameterisation result. + """ + + posterior_resamples = result.posterior.rvs(number_of_traces, apply_transform=True) + posterior_resamples_pdf = result.posterior.pdf(posterior_resamples) + pdf_range = np.asarray([posterior_resamples_pdf.min(), posterior_resamples_pdf.max()]) + + simulations = simulator(posterior_resamples) + + plot_dict = StandardPlot( + x=result.problem.domain_data, + y=result.problem.cost._dataset[dataset_y], + layout_options=layout_kwargs, + trace_names=data_legend_entry, + ) + + for pdf, pred in zip(posterior_resamples_pdf, simulations): + plot_dict.add_traces( + x=result.problem.domain_data, + y=pred, + line={ + 'dash': 'dot', + 'color': px.colors.sample_colorscale(colour_scale, (pdf - pdf_range[0]) / (pdf_range[1] - pdf_range[0]))[0], + } + ) + + # Add the colourbar. + plot_dict.add_traces( + x=[None], + y=[None], + mode='markers', + marker={ + 'size': 0, + 'color': pdf_range, + 'colorscale': colour_scale, + 'showscale': True, + 'colorbar': {'title': {'text': "Posterior PDF", 'side': 'right'}}, + }, + ) + + if pdf_plot is not None: + plot_dict.add_traces( + x=pdf_plot[0], + y=pdf_plot[1], + trace_names=pdf_label, + ) + + fig = plot_dict(show=False) + fig.update_layout(**layout_kwargs) + if show: + fig.show() + return fig + diff --git a/pyproject.toml b/pyproject.toml index 485f1893f..4542dd67a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -62,15 +62,18 @@ scifem = [ bpx = [ "bpx>=0.5, <0.6", ] -ep-bolfi = [ - "GPy>=1.13.2; python_version < '3.13'", - "ep-bolfi>=3.0.2; python_version < '3.13'" +examples = [ # *** temporarily include all of these dependencies to test the examples *** + "GPy>=1.13.2; python_version < '3.13'", # for ep-bolfi + "ep-bolfi>=3.0.2; python_version < '3.13'", + "sober-bo@git+https://github.com/ma921/SOBER.git@3cd4449a249294e0a8b9261bb32afde23ed8341c", # 10/12/25 + "pybammeis", + "seaborn", ] pyprobe = [ "PyProBE-Data>=2.5.0;python_version >= '3.11' and python_version < '3.13'" ] -all = ["pybop[plot,scifem,bpx,pyprobe,ep-bolfi]"] +all = ["pybop[plot,scifem,bpx,pyprobe,examples]"] [tool.setuptools.packages.find] include = ["pybop", "pybop.*"]