Batterijbeheer: State of Charge schatten

Een batterijmanagementsysteem weet twee dingen over een cel: hoeveel stroom er loopt, en welke spanning er op de klemmen staat. Geen van beide is de state of charge. De stroom geeft je de lading alleen als de stroomsensor de waarheid vertelt, en die vertelt je nooit of hij dat doet. De spanning zegt iets, maar alleen via een curve, en die curve is niet recht.

Deze casestudie bouwt de schatting die een echt BMS draait: een extended Kalman-filter dat naar beide signalen kijkt, de offset van het signaal dat het niet kan vertrouwen leert, en de lading meldt met zijn eigen onzekerheid.

De cel

Een 18650-lithium-ioncel als Thevenin-circuit: een ohmse weerstand, één polarisatietak, en een open-circuitspanning die van de lading afhangt.

De vier getallen hieronder zijn een uitgewerkt voorbeeld en geen datasheet. De eigen sheet van Panasonic voor de NCR18650B geeft 3,25 Ah minimum en begrenst de inwendige weerstand alleen vanaf boven met 100 mΩ; daar staat dus noch de 3,20 Ah noch de 0,05 Ω. Ze zijn gekozen voor de bemonstering in plaats van voor een cel: één pool op 0,05 rad/s, en twintig samples per tijdconstante.

V = V_{OC}(SOC) - V_1 - R_0 I \qquad \frac{dSOC}{dt} = -\frac{I}{3600\,Q_{cap}}, \qquad \frac{dV_1}{dt} = -\frac{V_1}{\tau} + \frac{I}{C_1}

Parameter	Betekenis	Waarde
$Q_{cap}$	Capaciteit	$3.20\ \text{Ah}$
$R_0$	Ohmse weerstand	$0.05\ \Omega$
$R_1$	Polarisatieweerstand	$0.02\ \Omega$
$C_1$	Polarisatiecapaciteit	$1000\ \text{F}$
$\tau = R_1C_1$	Tijdconstante van de polarisatie	$20\ \text{s}$

Twee toestanden. Eén ervan — de polarisatiespanning $V_1$ — is makkelijk, want dat is simpelweg het verschil tussen de klemspanning en wat de lading zegt dat de cel zou moeten aangeven. De andere is het probleem. De lading staat op geen meter.

Waarom tellen niet genoeg is

Integreer de stroom, deel door de capaciteit, trek af van de start. Geen model, geen curve, niets af te stellen, en een uur lang ziet het er perfect uit.

Het addertje is dat de meter niet de cel is. Een fout in de stroomsensor wordt lading die nooit terugkomt:

\frac{d\,SOC_{fout}}{dt} = \frac{\delta I}{3600\,Q_{cap}}

Belasting	Sensoroffset	$\delta I$	Drift
2 A	2 %	0,04 A	1.25 % per uur
10 A	2 %	0,20 A	6.25 % per uur
10 A	5 %	0,50 A	15.6 % per uur
2 A	10 %	0,20 A	6.25 % per uur
10 A	10 %	1,00 A	31.3 % per uur

Het tempo volgt de absolute offset $\delta I$ in ampere, en de capaciteit. Niet de belasting. De laatste twee regels zijn het punt, en de eerste twee zijn de valkuil: een percentage offset van een grotere belasting is een grotere absolute offset, dus de tabel leest alsof de drift met de belasting meegroeit. Houd de absolute offset op 0,20 A en het tempo is 6,25 % per uur, of de belasting nu 2 A of 10 A is. Een shunt wordt in procenten gespecificeerd, en daar zit de valkuil in.

Twee dingen volgen daaruit. De drift is een helling, en is dus onzichtbaar op een as van 0…100 % lading — een grafiek van de absolute lading ziet er de hele tijd perfect uit, en dat is precies de storing. En geen enkele belasting maakt het getal kleiner, dus het eerlijke antwoord is geen zachtere stroom. Het is een tweede mening over de lading.

Wat de spanning verklapt

Haal de stroom weg en de cel komt tot rust. De polarisatiespanning vervalt, en wat op de klemmen overblijft is de open-circuitspanning — een functie van de lading en van niets anders.

V = V_{OC}(SOC)

Dit is de ingang tot het probleem. Maar de curve is niet recht, en hoe recht hij is varieert enorm:

SOC	$\partial V_{OC}/\partial SOC$	tegenover het plateau
0–10 %	5.00 V per eenheid lading	8.3×
10–20 %	1.50	2.5×
20–90 %	0.60 – 0.90	het plateau
90–100 %	0.60	1.0×

Dezelfde 10 % lading is onderin de curve een halve volt waard en in het midden zes millivolt. Die factor 8.3 is geen curiositeit: hij bepaalt rechtstreeks de waarneembaarheid van het model. Met de meetrij $C = [\partial V_{OC}/\partial SOC,\; -1]$ en een diagonale toestandsmatrix is het determinant van de waarneembaarheidsmatrix

\det O = \left(\frac{\partial V_{OC}}{\partial SOC}\right)\left(1 - 0.9512\right)

Dus de lading is precies waarneembaar zolang de curve niet vlak is, en hoe goed waarneembaar ze is schaalt met die helling.

Een lineaire schatter, en waar hij faalt

De klassieke schatter neemt de meting, trekt af wat het model verwacht, en corrigeert de toestand met een versterking die uit de twee ruisniveaus volgt. Het is de juiste machinerie, en ze heeft één beslissing die van jou is: de relatie tussen lading en spanning is gebogen, dus de schatter lineariseert hem.

Eén keer. Het model wordt affien:

V_{OC}(SOC) \approx H_0\, SOC + \left[V_{OC}(s_{lin}) - H_0\, s_{lin}\right]

De helling $H_0$ is één keuze. De asafsnede tussen haakjes is een andere, en die bepaalt waar het filter vastthoudt — want een lineaire schatting convergeert naar het punt waar zijn eigen meetmodel de gemeten spanning reproduceert, en dat is niet waar de cel is.

Lineariseer op één werkpunt, zit in rust, en het resultaat is een systematische bias, geen vertraging:

Raaklijn op	Werkelijke lading	Blijvende bias
0.45	0.05	−53.3 %
0.45	0.85	+7.5 %
0.05	0.45	−28.8 %
0.05	0.05	−0.0 %

Met de raaklijn op 45 % en de werkelijke lading op 5 % komt de schatting tot rust op −0.48: een lading die geen cel kan hebben. Alleen een raaklijn dicht bij het werkelijke werkpunt geeft een kleine bias, dus één keuze van werkpunt is een beslissing waar je de hele run voor betaalt.

Er is een tweede, tegen-intuïtief detail. De schade komt van een raaklijn die te steil is, niet te vlak. Bij de verkeerde helling gelooft de schatting veel te veel in een ruisachtige meting, valt de winst in en kruipt de schatting — en waar ze stopt is ze nog 63 % van de waarheid ver, en dat is de laatste rij van de tabel hierboven. Een te vlakke raaklijn is vergeleken onschadelijk: de innovatie is dan groot genoeg om de schatting toch naar de waarheid te trekken.

De extended Kalman-filter

Aan de machinerie verandert niets. Eén regel wel — in plaats van de helling bij het begin vast te spijkeren, zoek je hem op waar de schatting nu staat:

H = \left[\left.\frac{\partial V_{OC}}{\partial SOC}\right|_{SOC_{hat}},\; -1\right]

De curve wordt niet een keer gelineariseerd en daarna vertrouwd. Hij wordt opnieuw gelineariseerd, op de plek waar de schatting is aangekomen, zodat de gevoeligheid die ze gebruikt de gevoeligheid is die er werkelijk is. De EKF is een Kalman-filter dat toegeeft dat het benadert.

Raaklijn op	Werkelijke lading	Lineaire bias	EKF-bias
0.45	0.05	−53.3 %	0.0 %
0.45	0.85	+7.5 %	−0.0 %
0.05	0.45	−28.8 %	−0.0 %
0.05	0.85	−63.1 %	+0.0 %

Eén set $Q$ en $R$ , beide gebieden, nergens bias. En de EKF is niet overal even snel: in het steile deel doet hij er 71 s over tegen de 2 s van de lineaire. De herlinearisatie is voorzichtig precies waar de meting werkelijk veel weet. Die afweging is het onderwerp van de afstelstap, geen gebrek dat verholpen moet worden.

Beginnen vanaf het verkeerde antwoord

Het interessante geval is niet de schatting in overleg met een geduldig experiment. Het is de schatting die de lading vindt uit het niets: schatting begint op 0.50, werkelijke lading 0.90, cel in rust, niets afgesteld om het makkelijk te maken.

Hij convergeert in ongeveer drie seconden. Niet "uiteindelijk" — de insnelcurve tegen de tijd is het hele punt van de stap, en de twee sporen liggen daarna op elkaar.

De EKF heeft één echte zwakheid, en die is het waard te benoemen. Het filter begint op het vlakke plateau, waar de curve weinig zegt, dus de eerste correctie wordt berekend met een kleine helling en schiet door tot in de steile staart. Het belandt waar de helling 8.3 keer zo groot is — maar de covariantie die daar al is ingestort, is voor de verkeerde helling berekend, en een ingestorte covariantie is een kleine versterking. Het herstelt, langzaam. Daarom is de insetteltijd ongeveer 3 s bij 0.45, 71 s bij 0.05, en 71 s bij 0.15.

Afstellen

$Q$ is wat het model toegeeft dat het niet weet. $R$ is wat de sensor toegeeft dat hij niet weet. De versterking is het dispuut tussen die twee, en er is geen instelling die goed is — alleen een die goed is voor de situatie waarin de cel verkeert.

Verhoog $R$ en de schatting leunt op haar model: glad, en traag met het opmerken van veranderingen. Verlaag het en ze leunt op de spanning: alert, en ze gelooft de ruis. $Q$ beweegt de andere kant van dezelfde balans, en dat is subtieler. De ladingstoestand staat dicht bij een zuivere integrator, dus zijn voornaamste fout is de gain van de stroomsensor — constant en multiplicatief in plaats van wit — en een additieve $q_{SOC}$ is helemaal geen sensorruis, maar toegegeven modelfout. Het is een dichtheid per seconde, dus het groeit mee:

$q_{SOC}$	σ na een uur	Wat het voorstelt
1e-6	6.0 % lading	meer spookonzekerheid dan de grootheid die je meet
1e-8	0.6 % lading	ongemodelleerde drift: temperatuur, veroudering

De vraag over de sampletijd

Beide filters draaien op 1 Hz, waar een BMS thuis is. Dat zit comfortabel binnen de eis in plaats van er net op: de snelste pool in de plant is $1/\tau = 0.05\ \text{rad/s} = 7.96\ \text{mHz}$ , wat 20 samples per tijdconstante geeft en 63× marge tegen de Nyquist-grens. De enige pool in het gediscretiseerde model zit op 0.9512, dus het filter is onvoorwaardelijk stabiel en er geldt geen stapgroottelimiet.

Plant en filter worden apart en exact gediscretiseerd, over hetzelfde interval:

A_d = \text{diag}\left(1,\; e^{-T_s/\tau}\right) = \text{diag}(1,\; 0.9512) \qquad B_d = \begin{bmatrix} -T_s/(3600\,Q_{cap}) \\ R_1\,(1 - e^{-T_s/\tau}) \end{bmatrix} = \begin{bmatrix} -8.681 \times 10^{-5} \\ 9.754 \times 10^{-4} \end{bmatrix}

Een exacte $A_d$ mengen met een vooruit-Euler- $B_d$ is een echte fout, en geen slordigheid: de tweede term ligt 2,52 % boven de exacte waarde (de exacte waarde ligt 2,46 % eronder), en dat verschil duikt op in de innovatie van de schatting op, waar het filter het niet kan onderscheiden van een fout in het model. Het blok herberekent de discretisatie — de ruis inbegrepen, via Van Loan — uit de parameters, zodat het veranderen van de capaciteit van de cel of van de sampletijd ook het eigen model van de schatting meeneemt.

Verder lezen

Het model waar deze gids op gebouwd is staat niet op deze pagina uit het niets. De thevenin-package van het Amerikaanse National Renewable Energy Laboratory geeft precies de vergelijkingen die deze gids geeft —

$\frac{d\,SOC}{dt} = \frac{-I\,\eta}{3600\,Q_{\max}} \qquad \frac{dV_j}{dt} = -\frac{V_j}{R_jC_j} + \frac{I}{C_j} \qquad V_{\text{cel}} = V_{\text{OCV}}(SOC) - \sum_j V_j - I R_0$

— met de hysteresehaak op nul, en dat is waar de twee-taksvorm hieronder voor is.

Randall, C. R. thevenin: Equivalent circuit models in Python (NREL, SWR-24-132).

BSD-3-Clause.

— de modelbeschrijving, met de vergelijkingen hierboven. De vier getallen uit deze gids komen daar niet vandaan; het zijn een werkeexample.

Thévenin-equivalentcircuitmodellering van lithiumcellen, en de uitbreiding met twee takken

die de tweede RC-tijdconstante toevoegt die een een-taksmodel laat vallen

Extended en iterated extended Kalman filtering voor batterijbeheer, en de rol van de

helling van de OCV-curve in de haalbare waarneembaarheid

De genormaliseerde innovatiekwadratuur, het getal dat vertelt of de aangenomen $Q$ en

$R$ ook de zijn die de cel heeft

Zelf proberen

De laatste stap van de interactieve casestudy is een speelveld, en de twee dingen die je wilt verslepen zijn de twee waar deze gids het over heeft:

Meetruisvariantie $R$ en procesruis $q_{SOC}$ , de twee grootheden uit de tabel hierboven. Beide staan op decadeschuifregelaars, want vier decades $R$ op een lineaire track zet elke bruikbare instelling in het eerste twintigste van de slag.
De werkelijke lading van de cel. Zet hem op 15 % en de schatting doet er 71 s over voor dat ene punt dat 90 % in 28 s haalt, om de reden in de alinea over het plateau.

De belasting, de offset van de sensor en de ruis van de sensor zelf staan er ook bij, en dat zijn de drie waar stap 2 van de les over gaat.