MATLAB BattMo: Incorrectly assigned variable and multiplication error (separate issues) when trying to use Optimisation with Dual Layer Capacity enabled

[szopaw]

Hello,

I've come across two errors when enabling "Interface.useDoubleLayerCapacity" and assigning "Interface.doubleLayerCapacitance" a value.

First error which occurs is:

Error using *
Incorrect dimensions for matrix multiplication. Check
that the number of columns in the first matrix
matches the number of rows in the second matrix. To
operate on each element of the matrix individually,
use TIMES (.*) for elementwise multiplication.

Error in [ * ](matlab:matlab.lang.internal.introspective.errorDocCallback('ADI/mtimes', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\core\utils\ADI.m', 188)) ([line 188](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\core\utils\ADI.m',188,0)))
h.val = uh.val;
^
Error in [Interface/updateReactionCapacityRateEquation](matlab:matlab.lang.internal.introspective.errorDocCallback('Interface/updateReactionCapacityRateEquation', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Electrochemistry\Interface.m', 437)) ([line 437](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Electrochemistry\Interface.m',437,0)))
state.capacityRequation = jDL - (cDL/(Fdt))*((dphi - dphi0) + (R.T./F)./c.(c - c0));
^

In my own simulation, with my own json setup, I've managed to override this by changing the "" operator to ".". I can then run the simulation itself.

However, when trying to run the Examples with no other changes aside from:

json.(pe).(co).(am).(itf).useDoubleLayerCapacity= true;
jsonstruct.(pe).(co).(am).(itf).doubleLayerCapacitance = 10;
json.(ne).(co).(am).(itf).useDoubleLayerCapacity= true;
jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 10;

This issue is permanent and implementing the above change makes no difference.

A second issue occurs when trying to run the optimisation based on in "runParameterIdentification.m", please see attached MATLAB scripts. This works as intended without the Double Layer enabled, however when it is enabled (and the above fix to Interface is implemented), it runs into this error:

Solving timestep 1/5: -> 1 Second
Solving timestep 2/5: 1 Second -> 2 Seconds
Solving timestep 3/5: 2 Seconds -> 3 Seconds
Solving timestep 4/5: 3 Seconds -> 4 Seconds
Solving timestep 5/5: 4 Seconds -> 5 Seconds
*** Simulation complete. Solved 5 control steps in 185 Milliseconds ***
Solving timestep 1/5: -> 1 Second
Solving timestep 2/5: 1 Second -> 2 Seconds
Solving timestep 3/5: 2 Seconds -> 3 Seconds
Solving timestep 4/5: 3 Seconds -> 4 Seconds
Solving timestep 5/5: 4 Seconds -> 5 Seconds
*** Simulation complete. Solved 5 control steps in 138 Milliseconds ***
Solving reverse mode step 1 of 5
Solving reverse mode step 2 of 5
Unrecognized field name "cElectrolyte".

Error in [Interface/updateReactionCapacityRateEquation](matlab:matlab.lang.internal.introspective.errorDocCallback('Interface/updateReactionCapacityRateEquation', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Electrochemistry\Interface.m', 433)) ([line 433](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Electrochemistry\Interface.m',433,0)))
c0 = state0.cElectrolyte;
^^^^^^^^^^^^^^^^^^^
Error in [BaseModel/getEquations](matlab:matlab.lang.internal.introspective.errorDocCallback('BaseModel/getEquations', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Utilities\Various\BaseModel.m', 640)) ([line 640](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Utilities\Various\BaseModel.m',640,0)))
eval(funcCallList{ifunc});
^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [PhysicalModel/getAdjointEquations](matlab:matlab.lang.internal.introspective.errorDocCallback('PhysicalModel/getAdjointEquations', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\ad-core\models\PhysicalModel.m', 230)) ([line 230](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\ad-core\models\PhysicalModel.m',230,0)))
[problem, state] = model.getEquations(state0, state, dt, forces, varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [PhysicalModel/solveAdjoint](matlab:matlab.lang.internal.introspective.errorDocCallback('PhysicalModel/solveAdjoint', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\ad-core\models\PhysicalModel.m', 960)) ([line 960](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\ad-core\models\PhysicalModel.m',960,0)))
problem_p = model.getAdjointEquations(current, after, dt_next, forces_p,...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [computeSensitivitiesAdjointADBattmo](matlab:matlab.lang.internal.introspective.errorDocCallback('computeSensitivitiesAdjointADBattmo', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Utilities\Adjoint\computeSensitivitiesAdjointADBattmo.m', 96)) ([line 96](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Utilities\Adjoint\computeSensitivitiesAdjointADBattmo.m',96,0)))
[lambda, lambdaVec]= setup.model.solveAdjoint(linsolve, getState, getObjectiveState, setup.schedule, lambdaVec, step);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [evalObjectiveBattmoEXPAVG](matlab:matlab.lang.internal.introspective.errorDocCallback('evalObjectiveBattmoEXPAVG', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\evalObjectiveBattmoEXPAVG.m', 139)) ([line 139](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\evalObjectiveBattmoEXPAVG.m',139,0)))
sens = computeSensitivitiesAdjointADBattmo(setupNew, states, parameters, getObj, ...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [GITTautooptim>@(p)evalObjectiveBattmoEXPAVG(p,objective,simsetup,params,'objScaling',objScaling)](matlab:matlab.lang.internal.introspective.errorDocCallback('GITTautooptim>@(p)evalObjectiveBattmoEXPAVG(p,objective,simsetup,params,''objScaling'',objScaling)', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\GITTautooptim.m', 141)) ([line 141](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\GITTautooptim.m',141,0)))
objectiveGradient = @(p) evalObjectiveBattmoEXPAVG(p, objective, simsetup, params, 'objScaling', objScaling);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [unitBoxBFGS](matlab:matlab.lang.internal.introspective.errorDocCallback('unitBoxBFGS', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\optimization\optim\unitBoxBFGS.m', 134)) ([line 134](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\Externals\mrst\autodiff\optimization\optim\unitBoxBFGS.m',134,0)))
[v0,g0] = f(u0);
^^^^^
Error in [GITTautooptim](matlab:matlab.lang.internal.introspective.errorDocCallback('GITTautooptim', 'C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\GITTautooptim.m', 167)) ([line 167](matlab: opentoline('C:\Users\szopa\OneDrive - Swansea University\Daniel Reinowski PhD resources\SINTEF\BattMo\GITT Diff\GITTautooptim.m',167,0)))
[v, pOptTmp, history] = unitBoxBFGS(p0scaled , ....
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Upon inspection, state0 at that point contains solely the capacityR field, and nothing else. This is only occurring at the third pass through Interface/updateReactionCapacityRateEquation during that step, assigning properly at twice before.

I'm trying to add the .m scripts here, but this form isn't allowing to upload them. I am however using the optimisation routine from "runParameterIdentification.m". I'm also not sure how to check the version of BattMo I'm using, but I've cloned it roughly three weeks ago, at the end of July. I'm happy to share the script and other information if you tell me how, I'd be grateful.

Kind Regards,
Daniel Reinowski

[xavierr]

Thanks for reporting this error. Could you share an example for the first error? The best would be a single input json file but it can also be a combination of matlab code and json input.

[szopaw]

Hello,

I just realised that I made an input error in the generic example (leaving the "jsonstruct" instead of "json"), which now seems to not throw up the multiplication error. It now also proceeds to the "Unrecognized field name "cElectrolyte"." error.

I still get that multiplication error in my script, but I cannot see a similar mistake in the code, so please see the .txt file for the matlab code and the .json for the input. I'm not able upload here the .mat file with the experimental data that goes into the code, so I'm not sure what to do about that (or if it's necessary to see the multiplication error).

I haven't made any other changes to the actual BattMo scripts other that what is already mentioned here.

digC20_GITT_model.json

GITTautooptim.txt

[xavierr]

Since I do not have access to you OCP curves, I used the ones from the Chen dataset.

I run a simple discharge.

I did not find any issue related to the double layer capacitance.

Have a look at the script that follow. You can also see how I looked at the input.

Could you then try it, and start from there to replicate the issue you have had with the double layer capacitance?

Otherwise, I can see three potential issues with your data set, which prevented me to run a successfull simulation in the first place.

1. Order in the guestStoichiometry input values. we should have guestStoichiometry100 > guestStoichiometry0 for the negative electrode and the opposite for the positive electrode.
2. The very low diffusion coefficient in the positive active material
3. The large thickness of the negative electrode creates problem

clear all

mrstModule add ad-core optimization mpfa mrst-gui

close all

ne      = 'NegativeElectrode';
pe      = 'PositiveElectrode';
elyte   = 'Electrolyte';
thermal = 'ThermalModel';
am      = 'ActiveMaterial';
co      = 'Coating';
itf     = 'Interface';
opc     = 'openCircuitPotential';
sd      = 'SolidDiffusion';
ctrl    = 'Control';
sep     = 'Separator';

% READ DATA FILES
jsonstruct = parseBattmoJson('digC20_GITT_model.json');

jsonstruct.(pe).(co).(am).(itf).useDoubleLayerCapacity = true;
jsonstruct.(ne).(co).(am).electronicConductivity       = 0.1;
jsonstruct.(pe).(co).(am).(itf).reactionRateConstant   = 1e-10;
jsonstruct.(pe).(co).(am).(itf).doubleLayerCapacitance = 10;

jsonstruct = removeStructFields(jsonstruct                               , ...
                                {ne, co, am, itf, 'openCircuitPotential'}, ...
                                {pe, co, am, itf, 'openCircuitPotential'}, ...
                                {elyte, 'ionicConductivity'}             , ...
                                {elyte, 'diffusionCoefficient'}          , ...
                                {'TimeStepping', 'timeStepDuration'}     , ...
                                {'TimeStepping', 'totalTime'});

%%
% Remove current collector. In 1D model, they usually do no make really sense.
jsonstruct.include_current_collectors = false;

%%
% load chen data set. We need those to get the OCPs
%
jsonstruct_chen = parseBattmoJson('chen2020_lithium_ion_battery.json');
jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json'));
jsonstruct_chen = mergeStructs({jsonstruct_chen, jsonstruct_geometry});

%%
% load LCO ocp. The OCP for the negative electrode in the Chen data set creates problem because it does
% not contain the large increase at the end of the discharge (when the active material at this electrode
% gets fully saturated).
dochangeOCP = true;
if dochangeOCP
    jsonstruct_itf = parseBattmoJson(fullfile('ParameterData', 'MaterialProperties', 'LCO', 'lco_Ramadass2004.json'));
    jsonstruct_chen.(pe).(co).(am).(itf) = jsonstruct_itf.(itf);
end

%%
% We merge the input
%

jsonstruct = mergeStructs({jsonstruct, jsonstruct_chen});

%%
% We use some handy tools to inspect the differences

fjv = compareStructs(jsonstruct, jsonstruct_chen);

fjv.print('trimLongEntries', true);

fjv.print('filter', {'comparison', 'different'});


%%
% The guest stoichiometric and saturation concentrations are correlated with the OCP. We remove those values from the input
jsonstruct = removeStructFields(jsonstruct                                , ...
                                {ne, co, am, itf, 'guestStoichiometry0'}  , ...
                                {ne, co, am, itf, 'guestStoichiometry100'}, ...
                                {pe, co, am, itf, 'guestStoichiometry0'}  , ...
                                {pe, co, am, itf, 'guestStoichiometry100'});

jsonstruct = removeStructFields(jsonstruct                                , ...
                                {ne, co, am, itf, 'saturationConcentration'}, ...
                                {pe, co, am, itf, 'saturationConcentration'});

%%
% The diffusion coefficient of the positive electrode in the inital dataset is too low. The simulation
% just cannot start. We remove them from input (and values from Chen will be used).
jsonstruct = removeStructFields(jsonstruct                                , ...
                                {pe, co, am, sd, 'referenceDiffusionCoefficient'});

%%
% The large thickness of the electrode appears to create problem. The simulation stops earlier. We remove
% it.
jsonstruct = removeStructFields(jsonstruct, ...
                                {ne, co, 'thickness'});

%%
% We merge the initial input with Chen input. The initial input takes precedence if the value is given in
% the two inputs.
jsonstruct = mergeStructs({jsonstruct, jsonstruct_chen}, 'warn', false);

%%
% We can visualize and compare the values
%
fjv = compareStructs(jsonstruct, jsonstruct_chen);
fjv.print('filter', {'comparison', 'different'});

%%
% We run the simulation
output = runBattery(jsonstruct);

%%

figure
plot(output.time/hour, output.E, '*-');
grid on
xlabel 'time  / h';
ylabel 'potential  / V';

fjv.print('filter', {'comparison', 'different'});

return

figure
plot(output.time/hour, output.I, '*-');
grid on
xlabel 'time  / h';
ylabel 'Current  / A';

[szopaw]

For thoroughness sake, I've downloaded the current version of BattMo and completely deleted the one I've been running just before writing this. I think I've been running a previous version, as I've had to change runBattery() to runBatteryJson in my code, it was not recognising the function anymore.

I just want to emphasise, the issues described occur only during the "Solving reverse mode" of the optimisation. Even my own data runs just fine in a single simulation, with no alterations and DL enabled.

Because I've stupidly decided to report two issue in one post, I'll specify:

The multiplication issue

As mentioned, in the default example code it was an error on my part. But it persists in my code after the refresh and after I've changed the issues you've highlighted. I've swapped the guestStoichiometry in both electrodes and reduced the anode thickness by an order of magnitude. I've multiplied the diffusivity by 100 in the initial guess and params. I genuinely cannot see where the issue is.

I'm attaching the modified code, .json plus a .zip of the .mat file with the data. I hope that's not a problem.

GITTautooptim.txt

digC20_GITT_model.json

BattMoGITTdigC20.zip

Changing the "*" to ".*" in Interface/updateReactionCapacityRateEquation line 437 still seems to sidestep it, but runs into the next issue.

The variable assignment issue

After the refresh, I can no longer replicate the issue in the default "runParameterIdentification.m" code. It just works now. But it still persists in mine, exactly as before. I have even copied again the entire section beginning "%% Setup the optimization problem" from "runParameterIdentification.m".

The difference between the two I can for sure identify, is that my experimental data consists solely of Potential and Time, formatted to fit into the "leastSquaresEI" function.

Of course just to confuse matters, when I deleted the old installation of BattMo and unzipped the new one, something seemed to work, no errors occurred. But it also worked with just runBattery() instead of runBatteryJson(). But afterwards, the script refused to recognise the function and I needed to rename it.

I honestly don't even know how that was possible, I'm 90% sure I've closed MATLAB before opening and running the new startupBattMo. I'm just really confused by now.

I'm sorry for this headache, but I just really want to be able to run the DL optimisation against my data. If there is any simpler way to input my cathode parameters against a metal anode and run it against the experimental data, I would be very happy to apply it.

[xavierr]

For some reasons, I do not receive the notifications from github so that I did not see you answer. Sorry for that . I will have a look.