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..100be0b8c --- /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) + ], + 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_parameterisation_with_sober.py b/examples/scripts/battery_parameterisation/bayesian_parameterisation_with_sober.py new file mode 100644 index 000000000..d5c14ce52 --- /dev/null +++ b/examples/scripts/battery_parameterisation/bayesian_parameterisation_with_sober.py @@ -0,0 +1,407 @@ +import importlib.util +import sys + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from pybamm import citations, print_citations +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, +) + +set_default_dtype(float64) +setting_parameters(device=device("cpu"), dtype=float64) + +seed = 0 + + +def visualise_correlation( + fig, + ax, + correlation, + names=None, + title=None, + cmap=plt.get_cmap("BrBG"), + 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. + """ + + # 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() + + +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 = np.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) * np.cos(n * np.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 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 * np.pi * x / self.L) + * np.exp(-(n**2) * np.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) + ) + + +class SiliconVoltageRelaxation(pybop.BaseSimulator): + def __init__( + self, parameters, fixed_parameters, timepoints: np.ndarray | None = None + ): + 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} + }""") + super().__init__(parameters) + self.fixed_parameters = fixed_parameters + self.timepoints = np.asarray(timepoints) + self.output_variables = ["Voltage change [V]"] + + def voltage_relaxation(self, inputs_array): + # Parameters from Köbbing2024; notation is taken from there. + R = 50e-9 + D = 1e-17 + E_core = 200e9 + nu_core = 0.22 + lambda_core = 64e9 + G_core = 82e9 + sigma_Y_core = 3e9 + c_max = 311e3 + SOC_0 = 0.1 * c_max + SOC_100 = 0.9 * c_max + v_Li = 9e-6 + L_shell = 20e-9 + E_shell = 100e9 + nu_shell = 0.3 + sigma_Y_shell = 2e9 + eta_shell = 135e12 + # sigma_ref = 133e6 + # tau = 3e8 + T = 298 + F = 96485 + R_gas = 8.314 + + diffusion_model = Diffusive_Relaxation(lambda x: x, L=1) + diffusion = False + + R_core = R - L_shell + alpha = 0.5 * (R_core / L_shell - 1) + + U_infty = inputs_array[0].reshape(-1, 1) # torch: unsqueeze(1) + slope = inputs_array[1].reshape(-1, 1) + timescale_exp = inputs_array[2].reshape(-1, 1) + if diffusion: + diff_portion = inputs_array[3].reshape(-1, 1) + diff_timescale = inputs_array[4].reshape(-1, 1) + else: + diff_portion = 0 + rel_magnitude = (1 - diff_portion) * U_infty + diff_magnitude = diff_portion * U_infty + # Express effective parameters by adjusting model parameters. + lambda_ch = 1 # does not appear independently + sigma_ref = slope * alpha * F * lambda_ch**3 / (2 * v_Li) + tau = lambda_ch * E_core * alpha * timescale_exp / sigma_ref + # Determine integration constant from t=0 with + # sigma_ev(t=0) = sigma_0 from sigma_ev = delta_U * F / v_Li. + integration_constant = np.tanh( + rel_magnitude * alpha * F * lambda_ch**3 / (2 * v_Li * sigma_ref) + ) + # sigma_0 = np.arctanh(integration_constant) * ( + # 2 * sigma_ref / (alpha * lambda_ch**3) + # ) + delta_U = ( + 2 + * v_Li + * sigma_ref + / (alpha * F * lambda_ch**3) + * np.arctanh( + integration_constant + * np.exp( + -E_core * alpha * lambda_ch / (tau * sigma_ref) * self.timepoints + ) + ) + ) + if diffusion: + delta_U += diffusion_model( + self.timepoints, -diff_magnitude, diff_timescale, diff_magnitude + ) + + return (U_infty - delta_U).T + + def batch_solve(self, inputs, calculate_sensitivities=False): + inputs_array = tensor([entry for entry in inputs[0].values()]) + relaxations = self.voltage_relaxation(inputs_array) + sols = [] + for entry, rel in zip(inputs, relaxations): + sol = pybop.Solution(entry) + sol.set_solution_variable("Voltage change [V]", rel) + sols.append(sol) + return sols + + +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 change [V]": np.asarray( + relaxations[data_index]["Voltage change [V]"] + ), + } + ) + len_data = len(t) + + short_term = pybop.MeanSquaredError( + dataset, + "Voltage change [V]", + [1] * short_term_end + [0] * (len_data - short_term_end), + ) + mid_term = pybop.MeanSquaredError( + dataset, + "Voltage change [V]", + [0] * short_term_end + + [1] * (long_term_start - short_term_end) + + [0] * (len_data - long_term_start), + ) + long_term = pybop.MeanSquaredError( + dataset, + "Voltage change [V]", + [0] * long_term_start + [1] * (len_data - long_term_start), + ) + + unknowns = pybop.MultivariateParameters( + { + "U(t=∞) [V]": pybop.Parameter( + initial_value=0.1, + bounds=[0.01, 0.2], + transformation=pybop.LogTransformation(), + ), + "log-slope [V]": pybop.Parameter( + initial_value=0.01, + bounds=[0.001, 0.2], + transformation=pybop.LogTransformation(), + ), + "relaxation timescale [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]]) + ), + ) + simulator = SiliconVoltageRelaxation(unknowns, {}, timepoints=t) + # 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, + ) + pybop_wrapper = SOBER_BASQ_EPLFI(problem, options) + pybop_result = pybop_wrapper.run() + # Calculate the correlation matrix in addition to the full plot. + raw_taken_samples = pybop_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)) + visualise_correlation( + fig_corr, + ax_corr, + raw_corr, + ["U(t=∞) [V]", "log-slope [V]", "relaxation timescale [s]"], + "", + entry_color="white", + ) + # Re-sample the posterior for the predictive posterior. + posterior_resamples = pybop_result.posterior.rvs(64, apply_transform=True) + posterior_resamples_pdf = pybop_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..e4e7f0618 --- /dev/null +++ b/examples/scripts/design_optimisation/bayesian_optimisation.py @@ -0,0 +1,487 @@ +# Miscellaneous imports for plotting, arithmetic, and statistics. +from copy import deepcopy + +# In this example, we use parallel processing. +from multiprocessing import Pool + +import matplotlib.pyplot as plt +import numpy as np +import pybamm +from matplotlib.ticker import PercentFormatter +from pybamm import citations, 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 import BaseSimulator, Parameters, Solution +from pybop.costs.feature_distances import indices_of +from pybop.optimisers.sober_basq_optimiser import SOBER_BASQ, SOBER_BASQ_Options + +set_default_dtype(float64) +setting_parameters(dtype=float64, device=device("cpu")) +np.seterr(divide="ignore") + + +class SEIGrowthVonKolzenberg(BaseSimulator): + """ + Implements the Solid-Electrolyte Interphase (SEI) growth model from + von Kolzenberg et al. (2020). + Extra model assumption: since the positive electrode is oversized, + it can deliver extra lithium to the negative electrode. + Currents are assumed to be slightly higher to account for this. + """ + + def __init__( + self, + parameters, + fixed_parameters, + timepoints: np.ndarray | None = None, + currents: np.ndarray | None = None, + ): + citations.register("""@article{ + vonKolzenberg2020, + title={{Solid-Electrolyte Interphase During Battery Cycling: Theory of Growth Regimes}}, + author={von Kolzenberg, L and Latz, A and Horstmann, B}, + journal={ChemSusChem}, + volume={13}, + pages={3901}, + year={2020}, + doi={10.1002/cssc.202000867} + }""") + super().__init__(parameters) + self.fixed_parameters = fixed_parameters + self.timepoints = timepoints + self.currents = currents + self.delta_t = np.diff(timepoints, append=timepoints[-1] - timepoints[-2]) + self.output_variables = ["SEI thickness [m]"] + + def sei_growth(self, inputs, timepoints=None, currents=None): + timepoints = timepoints if timepoints is not None else self.timepoints + currents = currents if currents is not None else self.currents + delta_t = np.diff(timepoints, append=timepoints[-1] - timepoints[-2]) + + local_parameters = deepcopy(self.fixed_parameters) + local_parameters.update(inputs) + + negative_electrode_capacity = ( + (1 - local_parameters["Negative electrode porosity"]) + * local_parameters["Negative electrode thickness [m]"] + * local_parameters["Electrode width [m]"] + * local_parameters["Electrode height [m]"] + * local_parameters["Maximum concentration in negative electrode [mol.m-3]"] + * 96485.33212 + ) + SOC_start = ( + local_parameters["Initial concentration in negative electrode [mol.m-3]"] + / local_parameters["Maximum concentration in negative electrode [mol.m-3]"] + ) + socs = SOC_start + np.cumsum(currents) * delta_t / negative_electrode_capacity + ocps = local_parameters["Negative electrode OCP [V]"](socs) + + R = 8.314462618 # Ideal gas constant + T = local_parameters["Ambient temperature [K]"] + F = 96485.33212 # Faraday's constant + eff_surface_area = ( + 3 + * (1 - local_parameters["Negative electrode porosity"]) + / local_parameters["Negative particle radius [m]"] + ) + # Assumption: exchange-current density is given as the pre-factor without concentration dependencies. + exchange_current_density = local_parameters[ + "Negative electrode exchange-current density [A.m-2]" + ]( + local_parameters["Initial concentration in electrolyte [mol.m-3]"], + 0.5, + local_parameters["Maximum concentration in negative electrode [mol.m-3]"], + local_parameters["Ambient temperature [K]"], + ).value + # Approximation: electrolyte concentration remains constant. + exchange_current_densities = ( + exchange_current_density * (1 - socs) ** 0.5 * socs**0.5 + ) + + # Approximation: the current sinked in the SEI is much smaller than the intercalation current. + intercalation_currents = currents / ( + eff_surface_area * local_parameters["Negative electrode thickness [m]"] + ) + intercalation_overpotential = ( + R + * T + / F + * 2 + * np.arcsinh(0.5 * intercalation_currents / exchange_current_densities) + ) + + diffusion_critical_thicknesses = ( + local_parameters["Initial concentration in electrolyte [mol.m-3]"] + * local_parameters["SEI diffusivity [m2.s-1]"] + * 96485.33212 + / local_parameters["SEI formation rate constant [A.m-2]"] + * np.exp( + -(1 - local_parameters["SEI formation symmetry factor"]) + * F + / (R * T) + * ( + intercalation_overpotential + + ocps + + local_parameters["SEI lithium reference potential [J.mol-1]"] / F + ) + ) + ) + migration_critical_thicknesses = ( + 2 + * R + * T + * local_parameters["SEI ionic conductivity [S.m-1]"] + / (F * np.abs(intercalation_currents)) + ) + + sei_thicknesses = [local_parameters["Inital SEI thickness [m]"]] + [ + 0.0 for _ in range(len(socs)) + ] + for i, dt in enumerate(delta_t): + current_sei_thickness = sei_thicknesses[i] + # This is basically SEI thickness minus tunneling effects, + # but adjusted to allow for SEIs thinner than the tunneling length. + apparent_sei_thickness = ( + current_sei_thickness - local_parameters["SEI tunneling length [m]"] + ) / 2 + ( + ( + ( + current_sei_thickness + - local_parameters["SEI tunneling length [m]"] + ) + / 2 + ) + ** 10 + + local_parameters["Reference apparent SEI thickness [m]"] ** 10 + ) ** 0.1 + sei_current = ( + -local_parameters["SEI formation rate constant [A.m-2]"] + * np.exp( + -local_parameters["SEI formation symmetry factor"] + * F + / (R * T) + * ( + intercalation_overpotential[i] + + ocps[i] + + local_parameters["SEI lithium reference potential [J.mol-1]"] + / F + ) + ) + * (1 + apparent_sei_thickness / migration_critical_thicknesses[i]) + / ( + 1 + + apparent_sei_thickness / migration_critical_thicknesses[i] + + apparent_sei_thickness / diffusion_critical_thicknesses[i] + ) + ) + sei_thicknesses[i + 1] = ( + current_sei_thickness + - local_parameters["SEI molar volume [m3.mol-1]"] / F * sei_current * dt + ) + + return np.asarray(sei_thicknesses[:-1]) + + def batch_solve(self, inputs, calculate_sensitivities=False): + sols = [] + for entry in inputs: + sol = Solution(entry) + sol.set_solution_variable("SEI thickness [m]", self.sei_growth(entry)) + sols.append(sol) + return sols + + +class IdealisedSolarBatteryDegradation(BaseSimulator): + """ + 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. + """ + + def __init__(self, parameters, fixed_parameters, capacity_cutoff=0.4): + super().__init__(parameters) + self.fixed_parameters = fixed_parameters + self.capacity_cutoff = capacity_cutoff + # Calculate the reference case of a minimally sized battery. + timepoints, currents = self.day_night_cycle(1.0) + self.sei_growth_model = SEIGrowthVonKolzenberg( + Parameters(), self.fixed_parameters, timepoints, currents + ) + self.eol_reference = self.eol(1.0) + self.output_variables = ["EOL [d]", "SEI thickness [m]"] + + def day_night_cycle(self, oversize_factor): + """ + :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 sei_growth(self, oversize_factor, return_plot_data=False): + timepoints, currents = self.day_night_cycle(oversize_factor) + sei_thicknesses = self.sei_growth_model.sei_growth({}, timepoints, currents) + if return_plot_data: + return sei_thicknesses, timepoints + return sei_thicknesses + + def eol_formula(self, sei_thicknesses): + relative_capacity_loss = ( + (sei_thicknesses - sei_thicknesses[0]) + * ( + 3 + * (1 - self.fixed_parameters["Negative electrode porosity"]) + / self.fixed_parameters["Negative particle radius [m]"] + ) + * 96485.33212 + / (3600 * self.fixed_parameters["Nominal cell capacity [A.h]"]) + ) + eol_day = indices_of(relative_capacity_loss, self.capacity_cutoff)[0] / 24 + return eol_day + + def eol(self, oversize_factor): + sei_thicknesses = self.sei_growth(oversize_factor) + return self.eol_formula(sei_thicknesses) + + def relative_eol_gain_formula(self, sei_thicknesses, oversize_factor): + gain_vs_offset = (self.eol_formula(sei_thicknesses) - self.eol_reference) / ( + oversize_factor - 1 + ) + return gain_vs_offset + + def relative_eol_gain(self, oversize_factor): + sei_thicknesses = self.sei_growth(oversize_factor) + return self.relative_eol_gain_formula(sei_thicknesses, oversize_factor) + + def batch_solve(self, inputs, calculate_sensitivities=False): + sols = [] + for entry in inputs: + sol = Solution(entry) + sei_thicknesses = self.sei_growth(entry["Oversize factor"]) + sol.set_solution_variable("SEI thickness [m]", sei_thicknesses) + sol.set_solution_variable( + "EOL [d]", + [ + self.relative_eol_gain_formula( + sei_thicknesses, entry["Oversize factor"] + ) + ], + ) + sols.append(sol) + return sols + + +if __name__ == "__main__": + battery_parameters = pybamm.ParameterValues("Marquis2019") + + # We append parameters from a cell we know the SEI parameters for. + battery_parameters.update( + { + "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 formation rate constant [A.m-2]": 1e-5, + "SEI ionic 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 formation symmetry factor": 0.22, + "Inital SEI thickness [m]": 2e-9, + "Reference apparent SEI thickness [m]": 0.05e-9, + "SEI tunneling length [m]": 2.05e-9, + "SEI molar volume [m3.mol-1]": 1.078e-5, + "SEI diffusivity [m2.s-1]": 1e-15, + "SEI thickness [m]": 67e-9, + "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.1, + "SEI Bruggeman coefficient": 4.54, + }, + check_already_exists=False, + ) + + pybop_prior = pybop.MultivariateParameters( + {"Oversize factor": pybop.Parameter(initial_value=1.1, bounds=[1.0, 1.2])}, + distribution=pybop.MultivariateUniform(np.asarray([[1.0, 1.2]])), + ) + # 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, + # such a plot can not be reasonably produced, but the optimiser + # still works just as well. + solar_battery_model = IdealisedSolarBatteryDegradation( + pybop_prior, battery_parameters + ) + oversize_factors = np.asarray([1 + 0.002 * i for i in range(101)]) + with Pool() as p: + eol_days = p.map(solar_battery_model.eol, oversize_factors) + oversize_factors = oversize_factors.flatten() + # End-Of-Life is chosen to be at 40% of capacity lost. + 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..71f123165 --- /dev/null +++ b/examples/scripts/model_selection/counting_kneepoints.py @@ -0,0 +1,269 @@ +import importlib.util +import sys + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import torch +from numpy import array, exp, linspace, log, ndarray, sum +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 + + +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 batch_solve(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): + 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.") + + +def marginalize_pdf(raw_taken_samples, sober_basq, x_eval=None): + """Calculate the marginal PDF of the kneepoint(s).""" + kde = gaussian_kde(raw_taken_samples) + dim = len(raw_taken_samples) + kneepoint_pdf_x = [] + kneepoint_pdf_y = [] + for i in range(1, dim, 2): + # Note: the rough order of parameters may be switched around. + # sober_basq.diag_order keeps track of that order. + raw_i = sober_basq.diag_order[i] + kneepoint_pdf_part_x = [] + kneepoint_pdf_part_y = [] + if x_eval is None: + k_p_p_min = kde.dataset[raw_i].min() + k_p_p_max = kde.dataset[raw_i].max() + raw_kneepoint_edges = linspace(k_p_p_min, k_p_p_max, 101) + else: + raw_kneepoint_edges = np.array( + [ + sober_basq.apply_transform_and_normalize_one_variable(x, i) + for x in x_eval + ] + ) + k_p_p_min = raw_kneepoint_edges.min() + k_p_p_max = raw_kneepoint_edges.max() + for k0, k1 in zip(raw_kneepoint_edges[:-1], raw_kneepoint_edges[1:]): + lower_bound = [kde.dataset[j].min() for j in range(dim)] + lower_bound[raw_i] = k0 + upper_bound = [kde.dataset[j].max() for j in range(dim)] + upper_bound[raw_i] = k1 + kneepoint_pdf_part_x.append( + sober_basq.denormalize_and_reverse_transform_one_variable( + 0.5 * (k0 + k1), + i, # this method uses diag_order + ) + ) + kneepoint_pdf_part_y.append(kde.integrate_box(lower_bound, upper_bound)) + k_norm = sum(kneepoint_pdf_part_y) * (k_p_p_max - k_p_p_min) + kneepoint_pdf_part_y = array(kneepoint_pdf_part_y) / k_norm + kneepoint_pdf_x.append(kneepoint_pdf_part_x) + kneepoint_pdf_y.append(kneepoint_pdf_part_y) + return kneepoint_pdf_x, kneepoint_pdf_y, kde + + +""" +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 = ndarray.astype(measurements[data_index]["Time [s]"].T[0], np.float64)[1:] + t[5] = t[4] + 1 + data = 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), + (array([0.0002, 1500, 0.0008]), array([0.0002, 1500, 0.0008, 2000, 0.001])), + ( + array([[0.00001, 0.001], [200, 2000], [0.00001, 0.01]]), + 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]", + ], + ), + ): + initial_values = exp(0.5 * (log(bounds.T[0]) + 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) + }, + 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, kneepoint_pdf_y, kde = marginalize_pdf( + pybop_result.posterior.distribution.distribution.dataset, + sober_basq_wrapper.optimiser, + x_eval=np.linspace(t[1], t[-1], 201), # zeroth is 0 + ) + # Sample the predictive posterior. + posterior_resamples = pybop_result.posterior.rvs(64, apply_transform=True) + posterior_resamples_pdf = pybop_result.posterior.pdf(posterior_resamples) + simulations = simulator(tensor(posterior_resamples)) + fig, ax = 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, sim in zip( + posterior_resamples, posterior_resamples_pdf, simulations + ): + ax.plot( + t, + simulator(torch.atleast_2d(tensor(pr)).T), + ls=":", + color=cmap(norm(pr_pdf)), + ) + degr_plot = ax.plot(t, data, color="black", lw=2, label="degradation data")[0] + ax_pdf = ax.twinx() + ax_pdf.plot(kneepoint_pdf_x[0], np.sum(kneepoint_pdf_y, axis=0)) + fig.colorbar( + matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap), + ax=ax_pdf, + label="Posterior PDF from KDE approximation", + ) + ax.set_xlabel("Cycles") + ax.set_ylabel("Capacity") + ax_pdf.set_ylabel("Kneepoints PDF") + ax.set_ylim((-0.1, 1.1)) + fig.legend( + [degr_plot], + [degr_plot.get_label()], + loc="outside lower center", + ) + + print_citations() + + plt.show() 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..cd06e746a --- /dev/null +++ b/examples/scripts/model_selection/impedance_model_selection.py @@ -0,0 +1,803 @@ +import importlib.util +import json +import sys +from contextlib import redirect_stdout +from multiprocessing import Pool + +import matplotlib.pyplot as plt +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 numpy import sqrt, tan +from scipy import constants +from scipy.stats import gaussian_kde +from torch import 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) + +import numpy as np +from sober import SoberWrapper +from torch import exp, log + +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), + }, +} + +parameter_transforms = { + "pouch_SPM": { + "Positive particle diffusivity [m2.s-1]": (lambda x: log(x), lambda x: exp(x)), + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + }, + "pouch_DFN": { + "Positive electrode Bruggeman coefficient": (lambda x: x, lambda x: x), + "Positive electrode conductivity [S.m-1]": (lambda x: log(x), lambda x: exp(x)), + "Electrolyte diffusivity [m2.s-1]": (lambda x: log(x), lambda x: exp(x)), + "Positive particle diffusivity [m2.s-1]": (lambda x: log(x), lambda x: exp(x)), + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + }, + "LG_MJ1_variant_1": { + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "SEI relative permittivity": (lambda x: log(x), lambda x: exp(x)), + "SEI Bruggeman coefficient": (lambda x: x, lambda x: x), + }, + "LG_MJ1_variant_2": { + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: log(x), + lambda x: exp(x), + ), + "SEI relative permittivity": (lambda x: log(x), lambda x: exp(x)), + "SEI Bruggeman coefficient": (lambda x: x, lambda x: x), + }, +} +parameter_transforms_numpy = { + "pouch_SPM": { + "Positive particle diffusivity [m2.s-1]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + }, + "pouch_DFN": { + "Positive electrode Bruggeman coefficient": (lambda x: x, lambda x: x), + "Positive electrode conductivity [S.m-1]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Electrolyte diffusivity [m2.s-1]": (lambda x: np.log(x), lambda x: np.exp(x)), + "Positive particle diffusivity [m2.s-1]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + }, + "LG_MJ1_variant_1": { + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "SEI relative permittivity": (lambda x: np.log(x), lambda x: np.exp(x)), + "SEI Bruggeman coefficient": (lambda x: x, lambda x: x), + }, + "LG_MJ1_variant_2": { + "Positive electrode double-layer capacity [F.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "Positive electrode exchange-current density [A.m-2]": ( + lambda x: np.log(x), + lambda x: np.exp(x), + ), + "SEI relative permittivity": (lambda x: np.log(x), lambda x: np.exp(x)), + "SEI Bruggeman coefficient": (lambda x: x, lambda x: x), + }, +} + + +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) * sqrt(εₙ ** (-βₙ) * 2 * pi * f / (2 * Dₑ)) + k_SEI = (1 - 1j) * 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) * tan( + k_electrolyte * L_electrolyte_for_SEI_model + ) * tan(k_SEI * L_SEI) + Z_D_SEI = ( + L_SEI + * Theta + / (Dₑ * ε_SEI**β_SEI) + * (t_SEI_minus - ρₑ_plus / ρₑ) ** 2 + * 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=[], +): + 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 + + 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=[], +): + 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)] + ).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))).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): + 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..df559ec52 --- /dev/null +++ b/examples/scripts/model_selection/visualize_large_variance_evidences.py @@ -0,0 +1,90 @@ +import matplotlib.pyplot as plt +from numpy import array, exp, log, logspace +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 = log(1 + exp(log_variance) / exp(log_evidence) ** 2) ** 0.5 + mu = log_evidence - sigma**2 / 2 + mean = exp(log_evidence) + lower = mean / exp(sigma) + upper = mean * 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 = array(log_evidences) + log_variances = array(log_variances) + min_extent = float("inf") + max_extent = -float("inf") + for le, lv, label in zip(log_evidences, log_variances, labels): + mu, sigma = calculate_distribution_parameters(le, lv) + lower = mu - 2 * sigma + upper = mu + 2 * sigma + evidence_eval = logspace(lower, upper, 200, base=exp(1)) + min_extent = min(lower, min_extent) + max_extent = max(upper, max_extent) + distribution = lognorm(s=sigma, scale=exp(mu)) + distribution_eval = distribution.pdf(evidence_eval) + ax.plot(evidence_eval, distribution_eval, label=label) + ax.set_xlim(0.9 * exp(min_extent), 1.1 * 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..6469d029a --- /dev/null +++ b/examples/scripts/surrogate_modelling/electrolyte_inverse_model_training.py @@ -0,0 +1,362 @@ +from copy import deepcopy +from multiprocessing import Pool + +import matplotlib +import matplotlib.pyplot as plt +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 numpy import arange, logspace, set_printoptions +from pandas import DataFrame +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 + +# torch.set_default_dtype(torch.float64) +# sober.setting_parameters(device=torch.device('cpu')) + + +def visualise_correlation( + fig, + ax, + correlation, + names=None, + title=None, + cmap=plt.get_cmap("BrBG"), + 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. + """ + + # 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(arange(len(correlation))) + ax.set_yticks(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() + + +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 = 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)) + visualise_correlation( + fig_corr, + ax_corr, + correlation.detach().cpu().numpy(), + names, + "Electrolyte characterization from pulse test", + entry_color="black", + ) + + 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) + ): + print() + print("Silhouette score of #" + str(i + 2) + " clusters:", score) + labelled_X = ( + 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)) + visualise_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..cdc6e9a46 --- /dev/null +++ b/examples/scripts/surrogate_modelling/global_inverse_modelling.py @@ -0,0 +1,172 @@ +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): + _, _, 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 68dab1a03..9c12e352c 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/optimisers/ep_bolfi_optimiser.py b/pybop/optimisers/ep_bolfi_optimiser.py index 4de567d83..9e89e1508 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): 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..92a424dbe --- /dev/null +++ b/pybop/optimisers/sober_basq_optimiser.py @@ -0,0 +1,452 @@ +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._result 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)} + 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()) + ] + 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()) + ] + 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, + logger=self._logger, + time=end - start, + optim_name="SOBER + BASQ", + posterior=posterior, + lower_bounds=None if self.bounds is None else self.bounds[:, 0], + upper_bounds=None if self.bounds is None else self.bounds[:, 1], + 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)} + 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()) + ] + 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()) + ] + 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, + logger=self._logger, + time=end - start, + optim_name="EP-LFI + SOBER + BASQ", + posterior=posterior, + lower_bounds=None if self.bounds is None else self.bounds[:, 0], + upper_bounds=None if self.bounds is None else self.bounds[:, 1], + 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])