An AEKF Battery State of Charge Estimation Method Combining MAP and Fuzzy Control

Abstract

The present invention discloses an AEKF battery state of charge estimation method combining MAP and fuzzy control. First, the least squares method is used to fit the distribution of the variance of the error innovation sequence to obtain the probability density function of the prior variance within the recognition window. The distribution change of the error innovation is estimated by the MAP method. According to the MAP estimation result, combined with a fuzzy controller, the input and output are fuzzified, and corresponding fuzzy control rules are formulated for fuzzy decision-making to adaptively update the window size in the AEKF algorithm. The present invention realizes the recognition of the moment when the error innovation distribution changes, and at the same time changes the window size based on the maximum posterior probability function difference and its change rate to adapt to the change of the error innovation distribution, improving the accuracy, robustness of the lithium-ion battery state estimation and the noise processing ability of the algorithm.

Description

Technical Field

[0001] The present invention belongs to the technical field of new energy vehicle and energy storage batteries, and particularly relates to a method for estimating the state of charge (SOC) of a battery by combining MAP and fuzzy control with an AEKF. Background Art

[0002] The state of charge (SOC) of a power battery refers to the relative proportion of the electric energy currently stored in the battery, usually expressed as a percentage of the remaining battery charge, ranging from 0% (fully discharged) to 100% (fully charged). Accurately estimating the SOC is crucial for electric vehicles (EVs) because it directly affects the driving range, battery management, and energy optimization. Precise SOC estimation can help optimize the charging and discharging strategies of the battery, extend the battery life, avoid performance degradation caused by overcharging and over-discharging, and at the same time improve the efficiency of the energy recovery system, enhancing the overall economy and sustainability.

[0003] The estimation of the state of charge (SOC) of a power battery faces multiple problems and challenges in practical applications. First, the SOC estimation of the battery is affected by various factors, such as battery aging, temperature changes, charging and discharging rates, and measurement noise, which make the accuracy of SOC estimation complex and unstable. Second, existing estimation methods, such as the open-circuit voltage method, the ampere-hour integration method, and the Kalman filtering method, all have certain limitations. For example, although the open-circuit voltage method is simple, it is greatly affected by temperature changes and requires the battery to be in a stationary state for accurate measurement; although the ampere-hour integration method has strong real-time performance, it is prone to error accumulation. Especially during long-term use, errors in the battery current sensor can lead to inaccurate SOC estimation; although the currently widely used Kalman filtering method can combine multiple sensor data to improve accuracy, the processing of process noise and measurement noise is still not ideal, and there may still be errors in complex environments. At the same time, improved filtering algorithms often encounter problems with poor robustness.

[0004] To address these challenges, data-driven methods (such as machine learning and deep learning) have been proposed as a new SOC estimation solution. By training models with a large amount of historical data, data-driven methods can improve the accuracy and adaptability of SOC estimation, but still need to solve problems such as the quality, real-time performance, and computational complexity of the model training data.

[0005] As described above, for the method for estimating the state of charge (SOC) of a power battery, there is currently a lack of a method for estimating the state of charge (SOC) of a power battery that is accurate, efficient, has strong robustness, and has strong noise processing ability. Summary of the Invention

[0006] To solve the above problems existing in the prior art, the present invention provides an AEKF state of charge estimation method combining MAP and fuzzy control (abbreviated as MFAEKF algorithm). It mainly estimates the distribution change of the error innovation through the MAP method. According to the MAP estimation result, the window size in the AEKF algorithm is adaptively updated by combining a fuzzy controller, realizing the identification of the moment when the error innovation distribution changes. At the same time, based on the difference of the maximum a posteriori probability function, the window size is changed through a fuzzy control strategy to adapt to the change of the error innovation distribution, improving the accuracy, robustness of the lithium-ion battery state estimation and the noise processing ability of the algorithm.

[0007] The object of the present invention is achieved by the following technical solutions:

[0008] An AEKF state of charge estimation method combining MAP and fuzzy control includes the following steps:

[0009] S1. Parameter initialization: Initialize the error innovation sequence array Ek1, the error innovation sequence array Ek2, and the variance value array respectively ;

[0010] S2. Fitting the variance distribution of the error innovation sequence: Add the error innovation data to the error innovation sequence array Ek1 and the error innovation sequence array Ek2 respectively, and use the least squares algorithm to fit the variance value of the error innovation sequence array Ek1 , and store the fitted variance estimation value in the variance value array ;

[0011] S3. Using the maximum a posteriori probability estimation method to determine the change of the error innovation sequence distribution: Take the middle moment of the array Ek2 as the test point for the change of the error innovation sequence distribution, and calculate the maximum a posteriori probability function value when the error innovation sequence distribution stored in the error innovation sequence array Ek2 does not change and the maximum a posteriori probability function value when it changes. Use and The difference and the relative change rate of the difference as the determination index for the change of the error innovation sequence distribution;

[0012] S4. Update the size of the recognition window L by combining a fuzzy control strategy: Determine the control quantity as the starting position strat_index of the updated recognition window L, and the observed quantity as the difference and the relative change rate ​; Fuzzify the input quantity and output quantity; formulate fuzzy control rules; calculate the fuzzy relation of each rule using the Mamdani algorithm and output and obtain the total output; defuzzify the control quantity according to the principle of the maximum membership degree.

[0013] Furthermore, the step S2 includes:

[0014] S21. Obtain the error innovation data , and add the error innovation data to the error innovation sequence array Ek1 and the error innovation sequence array Ek2 respectively, which are used for error innovation variance fitting and determination of the change in the error innovation sequence distribution; when the data length of the error innovation sequence array Ek1 reaches the set maximum length , use the least squares algorithm to fit the variance value of the error innovation sequence array Ek1 , and store the fitted variance estimate value in the variance value array ;

[0015] S22. When the data length of the variance value array reaches , use the least squares algorithm again to fit the variance distribution of the error innovation, and use this as the prior probability distribution of the error innovation variance value in the current calculation cycle ,

[0016]

[0017] wherein, is the error innovation variance distribution array, is its variance value;

[0018] S23. The probability density function of the error innovation sequence in the error innovation sequence array Ek2 is expressed as:

[0019]

[0020] wherein, is the error innovation sequence data, is the variance value of the error innovation sequence array.

[0021] Furthermore, the step S3 includes:

[0022] S31. Assume that the error innovation sequence distribution stored in the error innovation sequence array Ek2 has not changed, then the error innovation sequence it stores is subject to independent and identical distribution, and calculate the maximum a posteriori probability function value of the error innovation sequence distribution in the error innovation sequence array Ek2:

[0023]

[0024] Wherein, is the maximum a posteriori probability function of the error innovation sequence of the error innovation sequence array Ek2; is the error innovation sequence data, is the error innovation variance value, is the prior probability distribution of the error innovation variance value fitted in step S2; The preset maximum data length of the error innovation sequence array Ek2;

[0025] S32. Assume that the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed. Each time the recognition window length is updated, the middle moment in the error innovation sequence array Ek2 must be tested. If the tested time point is exactly the critical point of the change in the error innovation sequence distribution, then the maximum a posteriori probability function of the error innovation sequence in the array Ek2 is equal to the product of the probability density functions of the random variable sequences on both sides of the critical point (middle moment), and calculate the maximum a posteriori probability function value after the change in the error innovation sequence distribution :

[0026]

[0027] Wherein, is the maximum a posteriori probability function of the error innovation sequence of the error innovation sequence array Ek2; and are respectively the sequences before the critical point and the sequence after the critical point in the array Ek2; The preset maximum data length of the error innovation sequence array Ek2; is the prior probability distribution of the error innovation variance value fitted in step S2;

[0028] S33. Use and The difference and the relative change rate of the difference as the determination index for the change in the error innovation sequence distribution:

[0029]

[0030] If , then the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed;

[0031] Define the difference :

[0032]

[0033] In the formula, represents the relative change rate of the value, represents the valid value in the previous calculation period.

[0034] Furthermore, the step S4 includes:

[0035] S41. Set an identification window L with an initial length of N init , where the initial length N init is the same as the length of the error innovation sequence array Ek2. The identification window L is used to store the error innovation . The addition of the error innovation within the identification window L is synchronized with the arrays Ek2 and Ek1; Initialize the starting position start_index = 0 of the identification window L that the fuzzy controller needs to return;

[0036] S42. Determine the observed quantity and the control quantity: The control quantity is the starting position strat_index of the updated identification window L, and the observed quantity is and ;

[0037] S43. Fuzzify the input quantity and the output quantity:

[0038] Divide the difference into three fuzzy sets: zero, positive small, and positive large; Divide the change range of the difference into five levels: , , , , ; Obtain the fuzzy relationship table of the change of the difference ;

[0039] Divide the relative change rate into five fuzzy sets: negative large, negative small, zero, positive small, and positive large; Divide the change range of the relative change rate into seven levels: , , , , , , ; Obtain the fuzzy relationship table of the change of the relative change rate ;

[0040] Divide the starting position strat_index of the updated identification window into three fuzzy sets: zero, positive small, and positive large; Divide the change range of the control quantity into four levels: 0, 1 / 8Ninit , 1 / 4N init , 1 / 2N init ; Obtain the fuzzy table of the control quantity change;

[0041] S44. Formulate the fuzzy control rules to obtain the fuzzy control table;

[0042] S45. Use the Mamdani algorithm to calculate the fuzzy relationship and the output , and obtain the total output;

[0043] S46. Defuzzify the control quantity according to the principle of the largest membership degree, and convert the fuzzy control action into the final control output quantity (start_index) and input it into the AEKF state of charge estimator for SOC calculation.

[0044] Furthermore, the fuzzy rules formulated in step S44 are:

[0045] Rule (1): If the difference is 0 and the relative change rate is 0, then the control quantity is 0;

[0046] Rule (2): If the difference is positive small and the relative change rate is positive small, then the control quantity is positive small;

[0047] Rule (3): If the difference is positive large and the relative change rate is positive large, then the control quantity is positive large;

[0048] Rule (4): If the difference is positive large and the relative change rate is negative large, then the control quantity is positive small;

[0049] Rule (5): If the difference is positive large and the relative change rate is negative small, then the control quantity is positive large;

[0050] Rule (6): If the difference is positive small and the relative change rate is negative small, then the control quantity is positive small;

[0051] Rule (7); If the difference is positive small and the relative change rate is negative large, then the control quantity is positive small;

[0052] Rule (8): If the difference is positive large and the relative change rate is positive small, then the control quantity is positive small;

[0053] Rule (9): If the difference is small and positive, and the relative change rate is large and positive, then the control variable is small and positive.

[0054] Furthermore, the step S45 includes:

[0055] Calculating the fuzzy relation and output of each rule by using the Mamdani algorithm

[0056]

[0057] to obtain the total output:

[0058]

[0059] Compared with the prior art, the present invention has the following advantages:

[0060] The present invention introduces the MAP (Maximum A Posteriori Probability Estimation) method to determine the change of the noise distribution sequence. First, the distribution of the variance of the error innovation sequence is fitted by using the least squares method to obtain the probability density function of the prior variance within the recognition window. Using the difference between the maximum a posteriori probability function PP2 determined by both sides of the decision point (the middle moment of the error innovation sequence array Ek2) and the maximum a posteriori probability function PP1 within the overall recognition window as the feature of the change of the error innovation sequence distribution, and thus making the decision of updating the recognition window.

[0061] The present invention designs a fuzzy control strategy to update the recognition window L. By comparing the magnitude of the difference between PP2 and PP1 and the magnitude of its change rate, the input and output are fuzzified, and at the same time, corresponding fuzzy control rules are formulated for fuzzy decision-making, that is, the adaptive update of the length of the recognition window L is realized. Using the fuzzy rules as the control principle enables the update of the recognition window L to accurately track the change of the error innovation distribution within the recognition window. The ability of the algorithm to process noise is improved.

[0062] The fuzzy rules proposed by the present invention classify the change amplitude and relative change rate of the error innovation by using a hierarchical membership function, and establish 7 empirical fuzzy rules to ensure that the window size can accurately follow the change of the noise distribution. Compared with the fixed parameter adjustment method, this fuzzy control strategy avoids the limitations of manual parameter adjustment and enhances the adaptive ability.

[0063] The present invention adopts an adaptive extended Kalman filtering (MFAEKF) framework combining MAP and fuzzy control. The MPA method is used to identify the moment of distribution change. When the difference (or the change rate of the difference) between PP2 and PP1 is large, it indicates that the distribution of the error innovation within the identification window has changed significantly at this time, and more historical data (up to N init / 2) should be discarded to adapt to the change in its distribution; conversely, appropriate historical data should be retained to allow a transitional process for the change in the distribution of the error innovation. This can combine the probability modeling ability of MAP and the non-linear decision-making ability of fuzzy control to construct a hierarchical adaptive extended Kalman filtering framework, and the two cooperate to optimize the filtering process. Brief Description of the Drawings

[0064] Figure 1 It is the overall flowchart of the AEKF battery state of charge estimation method combining MAP and fuzzy control according to the present invention.

[0065] Figure 2 It is the battery terminal voltage change curve obtained from the experiment in the embodiment of the present invention.

[0066] Figure 3 It is the HPPC working condition excitation current curve obtained from the experiment in the embodiment of the present invention.

[0067] Figure 4 It is the Gaussian white noise curve added in the embodiment of the present invention.

[0068] Figure 5 It is a comparison chart of the SOC results estimated by the AEKF battery state of charge estimation method combining MAP and fuzzy control described in the embodiment of the present invention and the SOC results estimated by the traditional SEKF algorithm.

[0069] Figure 6 It is Figure 5 The partial enlarged view at J in

[0070] Figure 7 It is the schematic diagram of the adaptive adjustment of the error innovation window size in the embodiment of the present invention.

[0071] Figure 8 It is Figure 7 The partial enlarged view of the blue rectangular area in Detailed Embodiment

[0072] The present invention will be further explained below in combination with specific examples and drawings.

[0073] As Figure 1 shown, the present invention provides an AEKF battery state of charge estimation method combining MAP and fuzzy control, including the following steps:

[0074] S1. Initialization of MFAEKF algorithm parameters:

[0075] S11. Initialize the error innovation sequence array Ek1 and the error innovation sequence array Ek2 respectively. Among them, the error innovation sequence stored in the error innovation sequence array Ek1 is used for error information variance fitting, and the maximum data length of the array Ek1 is set to , and the error innovation sequence stored in the error innovation sequence array Ek2 is used for determining the change of sequence distribution. The maximum data length of the array Ek2 is set to ;

[0076] S12. Initialize the variance value array , which is used to store the variance values fitted from the error innovation sequence array Ek1. Set the maximum data length of the variance value array to be ;

[0077] S2. Fitting of the variance distribution of the error innovation sequence:

[0078] S21. Obtain the error innovation data , and the error innovation data is the difference between the measured terminal voltage value at the current moment and the estimated value of the MFAEKF algorithm; add the error innovation data to the error innovation sequence array Ek1 and the error innovation sequence Ek2 respectively, and use them for error innovation variance fitting and determination of the change of the error innovation sequence distribution; the data addition method for each array is to add one by one. If the number of elements in the array Ek1 or the array Ek2 exceeds the maximum length value of the array, the first element of the array should be removed and the new error innovation should be added to the end of the array; when the data length of the error innovation sequence array Ek1 reaches the set maximum length , that is, Ek1 = , use the least squares algorithm to fit the variance value of the error innovation sequence array Ek1 , and store the fitted variance estimate value in the variance value array ;

[0079] S22. When the data length of the variance value array reaches , that is , use the least squares algorithm again to fit the variance distribution of the error innovation (assuming that the error innovation and its variance distribution both satisfy a normal distribution with a mean of 0 and variances of and respectively), and use this as the prior probability distribution of the error innovation variance value in the current calculation loop ,

[0080]

[0081] In the formula, is the error innovation variance distribution array, is its variance value;

[0082] After fitting, clear the array and the array Ek1, and prepare for the operation of the next cycle;

[0083] S23. Assume that the error innovation sequence in the error innovation sequence array Ek2 follows a Gaussian white noise distribution, its mean is zero, and the variance is , then the probability density function of the error innovation sequence in the error innovation sequence array Ek2 can be expressed as:

[0084]

[0085] In the formula, is the error innovation series array, is its variance value.

[0086] S3. Use the MAP (Maximum A Posteriori Estimation) method to determine the change in the error innovation sequence distribution:

[0087] S31. Assume that the distribution of the error innovation sequence stored in the error innovation sequence array Ek2 has not changed, and calculate the maximum a posteriori probability function value of the error innovation sequence distribution in the error innovation sequence array Ek2 :

[0088] Assume that the distribution of the error innovation sequence stored in the error innovation sequence array Ek2 has not changed, then the error innovation sequences within the recognition window are independent and identically distributed, and the maximum a posteriori probability of the random variable sequence is:

[0089]

[0090] In the formula, is the maximum a posteriori probability function of the error innovation sequence of the array Ek2 (assuming that the distribution of the error innovation sequence of the array Ek2 has not changed); is the error innovation sequence data, is its variance value, is the prior probability distribution of the error innovation sequence variance value fitted in step S2.

[0091] Take the logarithm of both sides of the equation to obtain the maximum a posteriori probability function:

[0092]

[0093] Let , and solve to obtain:

[0094]

[0095] Substitute into the above formula, and the value of the logarithmic maximum a posteriori probability function can be obtained :

[0096]

[0097] S32. Assume that the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed, and calculate the value of the maximum a posteriori probability function after the change in the distribution of the error innovation sequence :

[0098] Assume that the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed. Each time the recognition window length is updated, the middle moment in the array Ek2 must be checked. If the time point being checked is exactly the critical point of the change in the distribution of the error innovation sequence, then the maximum a posteriori probability function of the error innovation sequence in the array Ek2 is equal to the product of the probability density functions of the random variable sequences on both sides of the critical point (middle moment):

[0099]

[0100] In the formula, is the maximum a posteriori probability function of the error innovation sequence of the array Ek2 (assuming that the distribution of the error innovation sequence of the array Ek2 has changed); the sequences before and after the critical point in the array Ek2 follow two different Gaussian distributions, with variances and respectively. Take the logarithm of the above formula to obtain its maximum a posteriori probability function:

[0101]

[0102] In the formula, is the error innovation sequence before the critical point, is the error innovation sequence after the critical point;

[0103] Let , , then there is:

[0104]

[0105]

[0106] Substitute the above two formulas into the formula, and the value of the maximum a posteriori probability function after the change in the distribution of the error innovation sequence is obtained :

[0107]

[0108] S33. Use and the difference and the relative change rate of the difference as the determination index for the change in the distribution of the error innovation sequence: As the determination index for the change in the distribution of the error innovation sequence:

[0109] By comparing and to determine whether the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed. If the distribution of the error innovation sequence in the array Ek2 has changed, then the value of the maximum a posteriori probability function will be greater than the value of the maximum a posteriori probability function ; conversely, the value will be greater than the value.

[0110] Therefore, the determination of the change in the distribution of the error innovation sequence can be expressed by and the difference as follows:

[0111]

[0112] For and make the following simplification:

[0113]

[0114]

[0115] Then:

[0116]

[0117] In the formula, , , and the values have all been obtained in the previous steps;

[0118] If , then the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed.

[0119] Define the difference :

[0120]

[0121] In the formula, represents the relative change rate of the value, represents the valid value.

[0122] by and As a basis for determining the distribution change of the error innovation sequence in the error innovation sequence array Ek2, the establishment of the determination index is completed.

[0123] S4. Update the recognition window size by combining fuzzy control strategy:

[0124] S41. Set an initial length to N init The recognition window L, initial length N init The same length as the error innovation sequence array Ek2, the identification window L is used to store the error innovation , identify the error information within the window L The addition is synchronized with the arrays Ek2 and Ek1; the starting position of the recognition window L that needs to be returned to initialize the fuzzy controller is start_index=0;

[0125] S42. Determine the observation and control amount: the control amount is the starting position strat_index of the updated identification window L (determines the identification window length used to update the error window function ), the observed quantity is and ;

[0126] S43. Fuzzify input and output:

[0127] The difference can be Divided into three fuzzy sets: zero (ZO), positive small (PS), positive large (PB), according to the difference The range of changes is divided into five levels: , , , , . Get the difference The fuzzy relationship table of the changes is shown in Table 1. When it is a negative value, it means that the distribution of the new information sequence within the recognition window has not changed. At this time, it is only necessary to keep the window length calculated at the last moment, so it is not included in the fuzzy control range.

[0128] Table 1 Difference Changing fuzzy relationship table

[0129]

[0130] In the same way, the relative rate of change can be It is divided into five fuzzy sets: Negative Big (NB), Negative Small (NS), Zero (ZO), Positive Small (PS), Positive Big (PB). According to the relative change rate it is divided into seven levels according to the change range: , , , , , , to obtain the fuzzy relation table of the relative change rate change, as shown in Table 2.

[0131] Table 2 Fuzzy relation table of relative change rate change

[0132]

[0133] The control quantity is the starting position strat_index of the updated recognition window. It is divided into three fuzzy sets: Zero (ZO), Positive Small (PS), Positive Big (PB). It is divided into four levels according to the change range of the control quantity: 0, 1 / 8N init , 1 / 4N init , 1 / 2N init to obtain the fuzzy table of the control quantity change, as shown in Table 3.

[0134] Table 3 Fuzzy table of control quantity change

[0135]

[0136] S44. Formulate fuzzy control rules:

[0137] Rule (1): If the difference is 0 and the relative change rate is 0, then the control quantity is 0;

[0138] Rule (2): If the difference is positive small and the relative change rate is positive small, then the control quantity is positive small;

[0139] Rule (3): If the difference is positive big and the relative change rate is positive big, then the control quantity is positive big;

[0140] Rule (4): If the difference is positive big and the relative change rate is negative big, then the control quantity is positive small;

[0141] Rule (5): If the difference is positive big and the relative change rate is negative small, then the control quantity is positive big;

[0142] Rule (6): If the difference is small and positive, and the relative change rate is small and negative, then the control variable is small and positive;

[0143] Rule (7); If the difference is small and positive, and the relative change rate is large and negative, then the control variable is small and positive;

[0144] Rule (8): If the difference is large and positive, and the relative change rate is small and positive, then the control variable is small and positive;

[0145] Rule (9): If the difference is small and positive, and the relative change rate is large and positive, then the control variable is small and positive.

[0146] According to the above fuzzy rules, a fuzzy control table can be obtained as shown in Table 4; in Table 4, is the difference between PP2 and PP1, is the relative change rate, and u is the starting position start_index that the updated recognition window L needs to intercept;

[0147] Table 4 Fuzzy control table

[0148]

[0149] S45. Calculate the fuzzy relation and output of each rule using the Mamdani algorithm, and obtain the total output:

[0150] The fuzzy relation of each rule is determined by the membership degree values of the input variables. For "Rule (2)":

[0151]

[0152] The fuzzy relations of the remaining rules are calculated in the same way as above.

[0153] According to the Mamdani algorithm:

[0154]

[0155] Taking "Rule (2)" as an example

[0156]

[0157] Calculate the output of each rule in the same way, and finally obtain the total output:

[0158]

[0159] S46. Defuzzify the control quantity according to the principle of the largest membership degree, and convert the fuzzy control action into the final control output quantity. (start_index) Input it into the AEKF state of charge estimator for SOC calculation.

[0160] Embodiment:

[0161] This embodiment is an AEKF battery state of charge estimation method combining MAP and fuzzy control, including the following steps:

[0162] Step S1. Initialize the parameters of the MFAEKF algorithm:

[0163] Initialize arrays Ek1 and Ek2 (with maximum lengths of 1000 and 200 respectively). Ek1 stores the error innovation sequence for variance fitting, and the error innovation sequence stored in Ek2 is used to determine its distribution change.

[0164] Initialize the array (with a maximum length of ) to store the fitted variance values.

[0165] Step S2. Fit the variance distribution of the error innovation sequence:

[0166] S21. Add the error innovation data (the error innovation data is the difference between the terminal voltage measurement value at the current moment and the estimated value of the MFAEKF algorithm) to arrays Ek1 and EK2 respectively (if the number of elements in Ek1 and Ek2 exceeds the maximum length value of the array, the first element of the array should be removed and the new error innovation should be added to the end of the array), and use them for the determination of the variance distribution change of the error innovation and the variance fitting of the error innovation respectively. When the array Ek1 reaches the preset maximum length of 1000 ， i.e., Ek1 = , use the least squares algorithm to fit the variance value of the Ek1 array ; use the array to store the variance estimation value fitted in ① ;

[0167] S22. When the array data exceeds 1000, i.e., , use the least squares algorithm to fit the variance distribution of the error innovation again (assuming that the error innovation and its variance distribution both satisfy the normal distribution with a mean of 0 and variances of and respectively), and use this as the prior probability distribution of the error innovation variance value in the current calculation loop :

[0168]

[0169] Wherein is the error innovation variance distribution array, is its variance value.

[0170] After fitting, clear the array and the array Ek1, and prepare for the operation of the next loop.

[0171] S23. Assume that the error innovation sequence in the Ek2 array follows a Gaussian white noise distribution, with a mean of zero and a variance of , so the probability density function of the error innovation sequence in the error innovation sequence array Ek2 can be expressed as:

[0172]

[0173] Wherein is the error innovation series array, is its variance value.

[0174] Step S3. Use the MAP (Maximum A Posteriori Probability Estimation) method to determine the change in the error innovation sequence distribution:

[0175] S31. Assume that the distribution of the error innovation sequence stored in the array Ek2 has not changed, then the error innovation sequence within the window follows independent and identically distributed, and the maximum a posteriori probability of the random variable sequence is:

[0176]

[0177] Wherein is the maximum a posteriori probability function of the error innovation sequence of the array Ek2 (assuming that the distribution of the error innovation sequence of the array Ek2 has not changed); is the error innovation sequence array, is its variance value, is the prior probability distribution of the variance value of the error innovation sequence fitted in step (2).

[0178] Take the logarithm of both sides of the equation to obtain the maximum a posteriori probability function:

[0179]

[0180] Let , and solve to obtain

[0181]

[0182] Substitute into the above formula, and the logarithmic maximum a posteriori probability function value can be obtained:

[0183]

[0184] S32. Assume that the distribution of the error innovation sequence in the array Ek2 has changed. Each time the recognition window length is updated, it is necessary to check the intermediate time in the array Ek2. If the time point being checked happens to be the critical point of the change in the distribution of the error innovation sequence, then the maximum a posteriori probability function of the error innovation sequence in the array Ek2 is equal to the product of the probability density functions of the random variable sequences on both sides of the critical point (intermediate time):

[0185]

[0186] In the formula, is the maximum a posteriori probability function of the error innovation sequence of the array Ek2 (assuming that the distribution of the error innovation sequence of the array Ek2 has changed); the sequences before and after the critical point in the array Ek2 follow two different Gaussian distributions, with variances of and respectively. Taking the logarithm of the above formula gives its maximum a posteriori probability function:[[]]END]]

[0187]

[0188] In the formula is the error innovation sequence before the critical point, is the error innovation sequence after the critical point.

[0189] Let , , then there is

[0190]

[0191]

[0192] Substitute the above two formulas The value of the maximum a posteriori probability function after the change in the distribution of the error innovation sequence can be obtained :[[]]END]]

[0193]

[0194] S33. By comparing and the size of, it can be determined whether the distribution of the error innovation sequence in the array Ek2 has changed. If the distribution of the error innovation sequence in the array Ek2 has changed, then the value of the maximum a posteriori probability function will be greater than the value of the maximum a posteriori probability function ; conversely, the value of will be greater than value. Therefore, the determination of the change in the distribution of the error innovation sequence can be expressed by the and difference as follows:

[0195]

[0196] For and make the following simplification:

[0197]

[0198]

[0199] Then:

[0200]

[0201] In the formula , , and values have all been obtained in the previous steps.

[0202] Define:

[0203]

[0204] In the formula represents the relative change rate of the value, represents the valid value in the previous calculation period.

[0205] Taking and as the judgment basis for the change in the distribution of the error innovation sequence in the error innovation sequence array Ek2, the establishment of the judgment index is completed here.

[0206] Step S4. Update the recognition window size in combination with the fuzzy control strategy:

[0207] S41. Set a recognition window L with an initial length of 200. The recognition window L is used to store the error innovation , and the addition of the error innovation in the recognition window L is synchronized with the arrays Ek2 and Ek1; Initialize the starting position start_index = 0 of the recognition window L that the fuzzy controller needs to return,

[0208] S42. Determine the observed quantity and the control quantity. The control quantity is the starting position strat_index of the updated recognition window (which determines the length of the recognition window used as the function for updating the error window ), and the observed quantity is and ;

[0209] S43. Fuzzify the input quantity and the output quantity.

[0210] The difference can be divided into three fuzzy sets: zero (ZO), positive small (PS), and positive big (PB). (Since when the difference is negative, it means that the distribution of the innovation sequence within the recognition window has not changed. At this time, only the window length calculated at the previous moment needs to be maintained, so it is not included in the scope of fuzzy control.) According to the change range of the difference , it is divided into five levels: 0, 2.5e4, 5e4, 7.5e4, 1e5. Obtain the fuzzy relation table of the change of the difference :

[0211] Table 5 Fuzzy relation table of the change of the difference in this embodiment

[0212]

[0213] In the same way, the relative change rate can be divided into five fuzzy sets: negative big (NB), negative small (NS), zero (ZO), positive small (PS), and positive big (PB). According to the change range of the relative change rate , it is divided into seven levels: -0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15. Obtain the fuzzy relation table of the change of the relative change rate :

[0214] Table 6 Fuzzy relation table of the change of the relative change rate in this embodiment

[0215]

[0216] The control quantity is the starting position start_index of the updated recognition window. It is divided into three fuzzy sets: zero (ZO), positive small (PS), and positive big (PB). According to the change range of the control quantity, it is divided into four levels: 0, 25, 50, 100. Obtain the fuzzy table of the change of the control quantity:

[0217] Table 7 Fuzzy table of the change of the control quantity in this embodiment

[0218]

[0219] S44. Fuzzy rule description:

[0220] Rule (1): If the difference is 0 and the relative change rate is 0, then the control quantity is 0; ​​

[0221] Rule (2): If the difference is small and positive, and the relative change rate is small and positive, then the control quantity is small and positive;

[0222] Rule (3): If the difference is large and positive, and the relative change rate is large and positive, then the control quantity is large and positive;

[0223] Rule (4): If the difference is large and positive, and the relative change rate is large and negative, then the control quantity is small and positive;

[0224] Rule (5): If the difference is large and positive, and the relative change rate is small and negative, then the control quantity is large and positive;

[0225] Rule (6): If the difference is small and positive, and the relative change rate is small and negative, then the control quantity is small and positive;

[0226] Rule (7); If the difference is small and positive, and the relative change rate is large and negative, then the control quantity is small and positive;

[0227] Rule (8): If the difference is large and positive, and the relative change rate is small and positive, then the control quantity is small and positive;

[0228] Rule (9): If the difference is small and positive, and the relative change rate is large and positive, then the control quantity is small and positive.

[0229] According to the above empirical rules, a fuzzy control table ( is the difference between PP2 and PP1, is the relative change rate, and u is the starting position start_index where the updated recognition window L needs to be intercepted):

[0230] Table 8 Fuzzy control table of this embodiment

[0231]

[0232] S45. Calculate the fuzzy relationship and output

[0233] The fuzzy relationship of each rule is determined by the membership degree values of the input variables. For rule 2

[0234]

[0235] The fuzzy relationship calculations for each of the remaining rules are the same as above.

[0236] According to the Mamdani algorithm:

[0237]

[0238] Taking Rule 2 as an example:

[0239]

[0240] Calculate the output of each rule in the same way, and finally obtain the total output

[0241]

[0242] S46. Defuzzify the control quantity according to the principle of the largest membership degree, and convert the fuzzy control action into the final control output quantity (start_index) Input it into the AEKF state of charge estimator for SOC calculation.

[0243] Apply the HPPC test data of 18650 batteries with a rated capacity of 2.2 Ah for algorithm verification. The model parameter identification adopts the FFRLS algorithm. Figure 2 and Figure 3 are the battery terminal voltage change and the HPPC working condition current excitation curve obtained from the experiment respectively. At the same time, add Gaussian white noise with a covariance of [0.0000001; 0.0000001; 0.000000001] and a mean of zero to the state variables. The initial SOC values are all set to 1. As Figure 4 shown is the added Gaussian white noise curve.

[0244] As Figure 5 shown is the comparison chart of the SOC results estimated by this method and the results estimated by the traditional AEKF algorithm. Due to the addition of Gaussian white noise to the state variables, the state estimation error of the AEKF algorithm is too large at 1471 s, the calculation becomes unstable, and finally the filtering diverges; while this method has a strong adaptive adjustment ability for noise and can keep the estimation result near the true value through the feedback adjustment of the recognition window. By comparing the errors with the reference SOC, the maximum error of this method is 0.0065, and the average error is 0.0021 (the estimation error after the algorithm is stable); the maximum error of the AEKF is 0.0125, and the average error is 0.0071. It can be seen from this that this method has stronger anti-noise interference ability and higher estimation accuracy than the traditional AEKF algorithm.

Claims

1. An AEKF battery state of charge estimation method combining MAP and fuzzy control, characterized in that: The following steps are involved: S1. Parameter initialization: Initialize the error innovation sequence array Ek1, error innovation sequence array Ek2 and variance value array θ respectively; S2. Error innovation sequence variance distribution fitting: The error innovation data e x Add them to the error innovation sequence array Ek1 and error innovation sequence array Ek2 respectively, and use the least squares algorithm to fit the variance value of the error innovation sequence array Ek1 And the fitted variance value Stored in the variance value array θ; S3. Use the maximum a posteriori probability estimation method to determine the change in the distribution of the error innovation sequence: Take the middle moment of the array Ek2 as the check point for the change in the distribution of the error innovation sequence, and calculate the maximum a posteriori probability function value PP1 when the distribution of the error innovation sequence stored in the error innovation sequence array Ek2 does not change and the maximum a posteriori probability function value PP2 when the distribution of the error innovation sequence changes, and use the difference δ between PP2 and PP1 and the relative change rate δ of the difference δ grad As a criterion for determining the change in the distribution of the error innovation sequence; S4. Update the size of the recognition window L in combination with the fuzzy control strategy: Determine the control amount as the starting position strat_index of the updated recognition window L, and the observation amount as the difference δ and the relative change rate δ grad ; Fuzzify the input and output; formulate fuzzy control rules; use Mamdani algorithm to calculate the fuzzy relationship R of each rule i And output U i And get the total output; The control quantity is defuzzified according to the maximum membership principle.

2. The AEKF battery state of charge estimation method combining MAP and fuzzy control as claimed in claim 1, characterized in that: The step S2 comprises: S21. Get error update data e x , the error information data e x They are added to the error innovation sequence array Ek1 and the error innovation sequence array Ek2 respectively, and used for error innovation variance fitting and error innovation sequence distribution change determination respectively; when the data length of the error innovation sequence array Ek1 reaches the set maximum length N k1 , use the least squares algorithm to fit the variance value of the error innovation sequence array Ek1 And the fitted variance value Stored in the variance value array θ; S22. When the data length of the variance value array θ reaches N θ , the least squares algorithm is used again to fit the variance distribution of the error innovation, and this is used as the prior probability distribution P(θ) of the variance value of the error innovation in the current calculation cycle. Where θ is the error innovation variance distribution array, is its variance value; S23. The probability density function of the error innovation sequence in the error innovation sequence array Ek2 is expressed as: In the formula, e k is the error innovation sequence data, is the variance value of the error innovation sequence array.

3. The AEKF battery state of charge estimation method combining MAP and fuzzy control as claimed in claim 1, characterized in that: The step S3 comprises: S31. Assuming that the distribution of the error innovation sequence stored in the error innovation sequence array Ek2 has not changed, the error innovation sequence stored therein is independent and identically distributed. Calculate the maximum posterior probability function value PP1 of the error innovation sequence distribution in the error innovation sequence array Ek2: In the formula, f e1 (ω) is the maximum a posteriori probability function of the error innovation sequence of the error innovation sequence array Ek2; k is the error innovation sequence data, is the error innovation variance value, P(θ) is the prior probability distribution of the error innovation variance value fitted in step S2; N k2 is the preset maximum data length of the error innovation sequence array Ek2; S32. Assume that the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed. Each time the identification window length is updated, the intermediate moments in the error innovation sequence array Ek2 are tested. If the tested time is the critical point of the error innovation sequence distribution change, the maximum a posteriori probability function of the error innovation sequence in the array Ek2 is equal to the product of the probability density functions of the random variable sequences on both sides of the critical point. The maximum a posteriori probability function value PP2 after the error innovation sequence distribution changes is calculated: In the formula, f e2 (ω) is the maximum a posteriori probability function of the error innovation sequence of the error innovation sequence array Ek2; and are the sequences e before the critical point in array Ek2 k1 and the sequence e after the critical point k2 The variance value of k2 is the preset maximum data length of the error innovation sequence array Ek2; P(θ) is the prior probability distribution of the error innovation variance value fitted in step S2; S33. Use the difference δ between PP2 and PP1 and the relative change rate δ of the difference δ grad As a criterion for determining the change in the distribution of the error innovation sequence: If δ>0, the distribution of the error innovation sequence in the error innovation sequence array Ek2 has changed; Define the relative rate of change of the difference δ: In the formula, δ grad Indicates the relative rate of change of δ value, δ - Indicates the effective delta value in the previous calculation cycle.

4. The AEKF battery state of charge estimation method combining MAP and fuzzy control as claimed in claim 1, characterized in that: The step S4 comprises: S41. Set an initial length of N init The recognition window L, initial length N init The same length as the error innovation sequence array Ek2, the identification window L is used to store the error innovation e x , identify the error information e within the window L x The addition is synchronized with the arrays Ek2 and Ek1; the starting position of the recognition window L that needs to be returned by the fuzzy controller is start_index=0; S42. Determine the observation and control amount: the control amount is the starting position strat_index of the updated recognition window L, and the observation amount is δ and δ grad ; S43. Fuzzify the input and output: The difference δ is divided into three fuzzy sets: zero, positive small, and positive large; according to the variation range of the difference δ, it is divided into five levels: δ1, δ2, δ3, δ4, and δ5; and the fuzzy relationship table of the variation of the difference δ is obtained; The relative rate of change δ grad Divided into five fuzzy sets: negative large, negative small, zero, positive small, positive large; according to the relative change rate δ grad The range of variation is divided into seven levels: grad1 , δ grad2 , δ grad3 , δ grad4 , δ grad5 , δ grad6 , δ grad7 ; Get the relative rate of change δ grad Changing fuzzy relationship table; The starting position strat_index of the updated recognition window is divided into three fuzzy sets: zero, positive small, and positive large; according to the range of control quantity, it is divided into four levels: 0, 1 / 8N init , 1 / 4N init , 1 / 2N init ; Obtain the fuzzy table of control quantity changes; S44. Formulate fuzzy control rules and obtain a fuzzy control table; S45. Use Mamdani algorithm to calculate the fuzzy relation R of each rule i And output U i , and get the total output; S46. Defuzzify the control quantity according to the maximum membership principle, and convert the fuzzy control action into the final control output u f (start_index) is input into the AEKF state of charge estimator for SOC calculation.

5. The AEKF battery state of charge estimation method combining MAP and fuzzy control as claimed in claim 4, characterized in that: The fuzzy rule formulated in step S44 is: Rule (1): If the difference δ is 0, the relative change rate δ grad If is 0, the control amount is 0; Rule (2): If the difference δ is positive and small, the relative rate of change δ grad If it is positively small, the control amount is positively small; Rule (3): If the difference δ is positive, the relative rate of change δ grad If it is positive, then the control amount is positive; Rule (4): If the difference δ is positive, the relative rate of change δ grad The more negative, the smaller the control amount; Rule (5): If the difference δ is positive, the relative rate of change δ grad The smaller the negative value, the larger the control amount; Rule (6): If the difference δ is positive and small, the relative rate of change δ grad The smaller the negative value is, the smaller the control amount is; Rule (7); if the difference δ is positive and small, the relative rate of change δ grad The more negative, the smaller the control amount; Rule (8): If the difference δ is positive, the relative rate of change δ grad If it is positively small, the control amount is positively small; Rule (9): If the difference δ is positive and small, the relative rate of change δ grad The larger the value, the smaller the control amount.

6. The AEKF battery state of charge estimation method combining MAP and fuzzy control as claimed in claim 5, characterized in that: The step S45 comprises: The Mamdani algorithm is used to calculate the fuzzy relation R of each rule. i And output U i : Get the total output:

Images