Skip to content

Commit

Permalink
tests: adds integration tests for eis fitting, update problem.evaluat…
Browse files Browse the repository at this point in the history
…e eis args, correct remainging cost functions for complex numbers
  • Loading branch information
BradyPlanden committed Aug 16, 2024
1 parent e571b8d commit 563737e
Show file tree
Hide file tree
Showing 7 changed files with 258 additions and 50 deletions.
86 changes: 43 additions & 43 deletions examples/scripts/eis_fitting.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,72 +5,72 @@
# Define model
parameter_set = pybop.ParameterSet.pybamm("Chen2020")
parameter_set["Contact resistance [Ohm]"] = 0.0
initial_state = {"Initial SoC": 0.5}
n_frequency = 20
f_eval = np.logspace(-4, 5, n_frequency)
model = pybop.lithium_ion.SPM(
parameter_set=parameter_set,
eis=True,
options={"surface form": "differential", "contact resistance": "true"},
)

model.build(
inputs={
"Negative electrode active material volume fraction": 0.531,
"Positive electrode active material volume fraction": 0.732,
},
initial_state=initial_state,
)

sim = model.simulateEIS(
inputs={
"Negative electrode active material volume fraction": 0.531,
"Positive electrode active material volume fraction": 0.732,
},
f_eval=f_eval,
)

# Fitting parameters
parameters = pybop.Parameters(
pybop.Parameter(
"Positive electrode thickness [m]",
prior=pybop.Gaussian(60e-6, 1e-6),
bounds=[10e-6, 80e-6],
"Negative electrode active material volume fraction",
prior=pybop.Uniform(0.4, 0.75),
bounds=[0.375, 0.75],
),
pybop.Parameter(
"Negative electrode thickness [m]",
prior=pybop.Gaussian(40e-6, 1e-6),
bounds=[10e-6, 80e-6],
"Positive electrode active material volume fraction",
prior=pybop.Uniform(0.4, 0.75),
bounds=[0.375, 0.75],
),
)

sigma0 = 1e-4


def noise(sigma, values):
# Generate real part noise
real_noise = np.random.normal(0, sigma, values)

# Generate imaginary part noise
imag_noise = np.random.normal(0, sigma, values)

# Combine them into a complex noise
return real_noise + 1j * imag_noise


# Form dataset
dataset = pybop.Dataset(
{
"Frequency [Hz]": np.logspace(-4, 5, 30),
"Current function [A]": np.ones(30) * 0.0,
"Impedance": np.asarray( # [0.74, 0.42]
[
1.22932096e-01 - 1.61334852e-01j,
1.20183857e-01 - 8.62168750e-02j,
1.13263444e-01 - 5.18047726e-02j,
1.04052155e-01 - 3.35022824e-02j,
9.64515013e-02 - 2.23021233e-02j,
9.06786107e-02 - 1.52884662e-02j,
8.63208234e-02 - 1.06615416e-02j,
8.31090848e-02 - 7.35620839e-03j,
8.10200368e-02 - 4.85799774e-03j,
8.00077590e-02 - 3.25064001e-03j,
7.96461450e-02 - 2.82671946e-03j,
7.94561760e-02 - 3.80233438e-03j,
7.90211481e-02 - 6.75677894e-03j,
7.73396194e-02 - 1.30368488e-02j,
7.10702300e-02 - 2.41923459e-02j,
5.33354614e-02 - 3.65349286e-02j,
2.70601682e-02 - 3.59150273e-02j,
1.04187830e-02 - 2.35001901e-02j,
4.46828406e-03 - 1.31671362e-02j,
2.19296862e-03 - 7.49687043e-03j,
8.55961473e-04 - 4.24936512e-03j,
2.47416663e-04 - 2.22977873e-03j,
6.24540680e-05 - 1.11430731e-03j,
1.51553128e-05 - 5.48240587e-04j,
3.64126863e-06 - 2.68650763e-04j,
8.72757867e-07 - 1.31515901e-04j,
2.09065980e-07 - 6.43673754e-05j,
5.00740475e-08 - 3.15013181e-05j,
1.19929932e-08 - 1.54164989e-05j,
2.87236101e-09 - 7.54468952e-06j,
]
),
"Frequency [Hz]": f_eval,
"Current function [A]": np.ones(n_frequency) * 0.0,
"Impedance": sim["Impedance"] + noise(sigma0, len(sim["Impedance"])),
}
)

signal = ["Impedance"]
# Generate problem, cost function, and optimisation class
problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
cost = pybop.SumSquaredError(problem)
cost = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma0)
optim = pybop.CMAES(cost, max_iterations=100, sigma0=0.25, max_unchanged_iterations=30)

x, final_cost = optim.run()
Expand Down
14 changes: 12 additions & 2 deletions pybop/costs/_likelihoods.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,12 @@ def compute(self, inputs: Inputs) -> float:
np.sum(
self._offset
+ self._multip
* np.sum((self._target[signal] - self.y[signal]) ** 2.0)
* np.sum(
np.real(
(self._target[signal] - self.y[signal])
* np.conj(self._target[signal] - self.y[signal])
)
)
)
for signal in self.signal
]
Expand Down Expand Up @@ -205,7 +210,12 @@ def compute(self, inputs: Inputs) -> float:
np.sum(
self._logpi
- self.n_data * np.log(sigma)
- np.sum((self._target[signal] - self.y[signal]) ** 2.0)
- np.sum(
np.real(
(self._target[signal] - self.y[signal])
* np.conj(self._target[signal] - self.y[signal])
)
)
/ (2.0 * sigma**2.0)
)
for signal in self.signal
Expand Down
3 changes: 1 addition & 2 deletions pybop/costs/base_cost.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,6 @@ def __init__(self, problem: Optional[BaseProblem] = None):
self.signal = self.problem.signal
self.transformation = self.parameters.construct_transformation()
self._has_separable_problem = True
self.eis = self.problem.eis

@property
def n_parameters(self):
Expand Down Expand Up @@ -91,7 +90,7 @@ def evaluate(self, inputs: Union[Inputs, list]):
try:
if self._has_separable_problem:
self.y = self.problem.evaluate(
inputs, eis=self.eis, update_capacity=self.update_capacity
inputs, update_capacity=self.update_capacity
)

return self.compute(inputs)
Expand Down
1 change: 1 addition & 0 deletions pybop/models/base_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -807,6 +807,7 @@ def new_copy(self):
"var_pts": self.pybamm_model.default_var_pts.copy(),
"spatial_methods": self.pybamm_model.default_spatial_methods.copy(),
"solver": self.pybamm_model.default_solver.copy(),
"eis": copy.copy(self.eis),
}

return model_class(**model_args)
Expand Down
2 changes: 1 addition & 1 deletion pybop/plotting/nyquist.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def nyquist(problem, problem_inputs: Inputs = None, show=True, **layout_kwargs):
else:
problem_inputs = problem.parameters.verify(problem_inputs)

model_output = problem.evaluate(problem_inputs, eis=problem.eis)
model_output = problem.evaluate(problem_inputs)
domain_data = model_output["Impedance"].real
target_output = problem.get_target()

Expand Down
3 changes: 1 addition & 2 deletions pybop/problems/fitting_problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,6 @@ def set_initial_state(self, initial_state: Optional[dict] = None):
def evaluate(
self,
inputs: Inputs,
eis=False,
update_capacity=False,
) -> dict[str, np.ndarray[np.float64]]:
"""
Expand All @@ -121,7 +120,7 @@ def evaluate(
The model output y(t) simulated with given inputs.
"""
inputs = self.parameters.verify(inputs)
if eis:
if self.eis:
return self._evaluateEIS(inputs, update_capacity=update_capacity)
else:
try:
Expand Down
199 changes: 199 additions & 0 deletions tests/integration/test_eis_parameterisation.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,199 @@
import numpy as np
import pytest

import pybop


class TestEISParameterisation:
"""
A class to test the eis parameterisation methods.
"""

@pytest.fixture(autouse=True)
def setup(self):
self.sigma0 = 5e-4
self.ground_truth = np.asarray([0.55, 0.55]) + np.random.normal(
loc=0.0, scale=0.05, size=2
)

@pytest.fixture
def model(self):
parameter_set = pybop.ParameterSet.pybamm("Chen2020")
x = self.ground_truth
parameter_set.update(
{
"Negative electrode active material volume fraction": x[0],
"Positive electrode active material volume fraction": x[1],
}
)
return pybop.lithium_ion.SPM(
parameter_set=parameter_set,
eis=True,
options={"surface form": "differential"},
)

@pytest.fixture
def parameters(self):
return pybop.Parameters(
pybop.Parameter(
"Negative electrode active material volume fraction",
prior=pybop.Uniform(0.4, 0.75),
bounds=[0.375, 0.75],
),
pybop.Parameter(
"Positive electrode active material volume fraction",
prior=pybop.Uniform(0.4, 0.75),
bounds=[0.375, 0.75],
),
)

@pytest.fixture(params=[0.5])
def init_soc(self, request):
return request.param

@pytest.fixture(
params=[
pybop.GaussianLogLikelihoodKnownSigma,
pybop.GaussianLogLikelihood,
pybop.RootMeanSquaredError,
pybop.SumSquaredError,
pybop.SumofPower,
pybop.Minkowski,
pybop.MAP,
]
)
def cost(self, request):
return request.param

def noise(self, sigma, values):
# Generate real part noise
real_noise = np.random.normal(0, sigma, values)

# Generate imaginary part noise
imag_noise = np.random.normal(0, sigma, values)

# Combine them into a complex noise
return real_noise + 1j * imag_noise

@pytest.fixture(
params=[
pybop.SciPyDifferentialEvolution,
pybop.CMAES,
pybop.CuckooSearch,
pybop.NelderMead,
pybop.SNES,
pybop.XNES,
]
)
def optimiser(self, request):
return request.param

@pytest.fixture
def optim(self, optimiser, model, parameters, cost, init_soc):
n_frequency = 12
# Set frequency set
f_eval = np.logspace(-4, 5, n_frequency)

# Form dataset
solution = self.get_data(model, init_soc, f_eval)
dataset = pybop.Dataset(
{
"Frequency [Hz]": f_eval,
"Current function [A]": np.ones(n_frequency) * 0.0,
"Impedance": solution["Impedance"]
+ self.noise(self.sigma0, len(solution["Impedance"])),
}
)

# Define the problem
signal = ["Impedance"]
problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)

# Construct the cost
if cost in [pybop.GaussianLogLikelihoodKnownSigma]:
cost = cost(problem, sigma0=self.sigma0)
elif cost in [pybop.GaussianLogLikelihood]:
cost = cost(problem, sigma0=self.sigma0 * 4) # Initial sigma0 guess
elif cost in [pybop.MAP]:
cost = cost(
problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0
)
elif cost in [pybop.SumofPower, pybop.Minkowski]:
cost = cost(problem, p=2)
else:
cost = cost(problem)

# Construct optimisation object
common_args = {
"cost": cost,
"max_iterations": 250,
"absolute_tolerance": 1e-6,
"max_unchanged_iterations": 25,
}

if isinstance(cost, pybop.MAP):
for i in cost.parameters.keys():
cost.parameters[i].prior = pybop.Uniform(
0.2, 2.0
) # Increase range to avoid prior == np.inf

# Set sigma0 and create optimiser
sigma0 = 0.05 if isinstance(cost, pybop.MAP) else None
optim = optimiser(sigma0=sigma0, **common_args)

return optim

@pytest.mark.integration
def test_eis_optimisers(self, optim):
x0 = optim.parameters.initial_value()

# Add sigma0 to ground truth for GaussianLogLikelihood
if isinstance(optim.cost, pybop.GaussianLogLikelihood):
self.ground_truth = np.concatenate(
(self.ground_truth, np.asarray([self.sigma0]))
)

initial_cost = optim.cost(x0)
x, final_cost = optim.run()

# Assertions
if np.allclose(x0, self.ground_truth, atol=1e-5):
raise AssertionError("Initial guess is too close to ground truth")

if isinstance(optim.cost, pybop.GaussianLogLikelihood):
np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2)
np.testing.assert_allclose(x[-1], self.sigma0, atol=5e-4)
else:
assert (
(initial_cost > final_cost)
if optim.minimising
else (initial_cost < final_cost)
)
np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2)

def get_data(self, model, init_soc, f_eval):
initial_state = {"Initial SoC": init_soc}
model.build(
inputs={
"Negative electrode active material volume fraction": self.ground_truth[
0
],
"Positive electrode active material volume fraction": self.ground_truth[
1
],
},
initial_state=initial_state,
)
sim = model.simulateEIS(
inputs={
"Negative electrode active material volume fraction": self.ground_truth[
0
],
"Positive electrode active material volume fraction": self.ground_truth[
1
],
},
f_eval=f_eval,
)

return sim

0 comments on commit 563737e

Please sign in to comment.