Optimizing capacity fade and current using pymoo

[cyasing]

I am not sure if I should ask this here, but:

I am trying to conduct optimization of capacity fade while current gets be maximized (2 objectives). I implemented a simple setup such that the cell first gets discharged completely, or by half (depending on C rate and time), and then the pybamm setup is basically written inside pymoo setup such that the rest of charging can be observed by pymoo's algorithm.

Unfortunately, setting an arbitrary lower bound for current makes it so that it doesn't solve fully. It gives IDAS_TOO_MUCH_WORK, or IDAS_CONV_FAIL, or even pybamm.expression_tree.exceptions.SolverError: Events ['Minimum voltage'] are non-positive at initial conditions.

Unless I guess some right negative current (charging) values, it doesn't run. But then, how to determine the exact bounds for which the current would work?

I'm unable to understand how to set it up. If anyone here could tell how to proceed it will be best.

The bare pybamm setup looks like this:

import pybamm
import numpy as np

model_0 = pybamm.lithium_ion.DFN()
param = model_0.default_parameter_values
experiment = pybamm.Experiment(["Discharge at 1C for 0.5 hour"])
solver = pybamm.CasadiSolver(mode="fast")
sim = pybamm.Simulation(model_0, parameter_values=param, experiment=experiment, solver=solver)
solution_0 = sim.solve()
sim.plot()

model = model_0.set_initial_conditions_from(solution_0)
param["Current function [A]"] = "[input]"

solver = pybamm.CasadiSolver(mode="fast")
t_eval = np.linspace(0,1800)

sim = pybamm.Simulation(model, parameter_values=param, solver=solver)
solution = sim.solve(t_eval, inputs={"Current function [A]": -2 })
capacity_fade = solution["Total lithium capacity [A.h]"].entries # Rudimentary for now
current_val = solution["Current [A]"].entries
sim.plot()

With pymoo added, it looks like this:

import pybamm
import numpy as np
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize
from pymoo.operators.sampling.rnd import FloatRandomSampling
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM

class CellOptimization(Problem):
    def __init__(self, model_0, solution_0):
        super().__init__(n_var=1, n_obj=2, n_constr=0, xl=-1, xu=5, model_0=model_0, solution_0=solution_0)

    def _evaluate(self, x, out, *args, **kwargs):
        model = model_0.set_initial_conditions_from(solution_0)
        param = model.default_parameter_values
        param["Current function [A]"] = "[input]"

        solver = pybamm.CasadiSolver(mode="fast")
        t_eval = np.linspace(0,3600)

        sim = pybamm.Simulation(model, parameter_values=param, solver=solver)
        solution = sim.solve(t_eval, inputs={"Current function [A]": x[0]})
        
        capacity_fade = solution["Total lithium capacity [A.h]"].entries # Rudimentary for now
        current_val = solution["Current [A]"].entries
        out["F"] = [-current_val, -capacity_fade] # To maximize both

model_0 = pybamm.lithium_ion.DFN()
param = model_0.default_parameter_values
experiment = pybamm.Experiment(["Discharge at 1C for 1 hour"])
solver = pybamm.CasadiSolver(mode="fast")
sim = pybamm.Simulation(model_0, parameter_values=param, experiment=experiment, solver=solver)
solution_0 = sim.solve()
sim.plot()

problem = CellOptimization(model_0=model_0, solution_0=solution_0)

algorithm = NSGA2(
    pop_size=100,
    n_offsprings=100,
    sampling=FloatRandomSampling(),
    crossover=SBX(prob=0.9, eta=15),
    mutation=PM(eta=20),
    eliminate_duplicates=True
)

res = minimize(problem, algorithm, verbose=True)

print("Generations", res.algorithm.n_gen)
plot = Scatter(title="Problem")
plot.add(problem.pareto_front(use_cache=False, flatten=False), plot_type="line", color="black")
plot.add(res.F, facecolor="none", edgecolor="red", alpha=0.8, s=20)
plot.show()

Note, I did look at calculate_sensitivites feature, but I was not sure how to go about multi-objective optimization.

[cyasing]

Okay, if you control the C-rate and manage time of simulation depending on C-rate, it should work.