A public repository for BMS algorithms for a NUCLEO BOARD.
Hi Everyone!
Welcome to this repository from Howey's Research Group at the University of Oxford.
The code published here incorporates BMS algorithms for diagnosis functions such as SOC, SOH and Power estimation on a Kokam 53Ah Li-ion battery. This code was designed to work with a NUCLEO F401-RE board and to be tested with a dSPACE HIL Simulator. A short guide on how the set up works is available at https://bitbucket.org/ff95/bms .
The code is made up of three key parts. "Headers" and "Source" folders and the "main.cpp" file. As the code was generated by converting a Simulink model ( available on the BitBucket page), the headers and source code files generated by the conversion are in the corresponding "Headers" and "Source" folders. The "main.cpp" file sets up the ADC, the USB data transmission and starts the estimation (once a character "y" has been received by the computer it is connected to). It also transmits the data from the estimation via USB. Explanation on how to set up the communication with the board is available at BitBucket webpage, from where a MATLAB file can be downloaded which allows real time communication.
For any questions you can contact the author at federicomariaferrari@gmail.com .
The Simulink and Matlab files, together with a short guide, are all available at: https://bitbucket.org/ff95/bms.
Thanks for trying this out!
Federico
source/Simulink/EKF.cpp
- Committer:
- fmferrari
- Date:
- 2016-12-13
- Revision:
- 12:bf3e3b87224e
- Parent:
- 11:e39c490420d2
- Child:
- 13:831eab218c33
File content as of revision 12:bf3e3b87224e:
// // Academic License - for use in teaching, academic research, and meeting // course requirements at degree granting institutions only. Not for // government, commercial, or other organizational use. // // File: EKF.cpp // // Code generated for Simulink model 'EKF'. // // Model version : 1.83 // Simulink Coder version : 8.11 (R2016b) 25-Aug-2016 // C/C++ source code generated on : Tue Dec 13 11:11:00 2016 // // Target selection: ert.tlc // Embedded hardware selection: STMicroelectronics->ST10/Super10 // Code generation objectives: Unspecified // Validation result: Not run // #include "EKF.h" #include "EKF_private.h" // Block states (auto storage) DW_EKF_T EKF_DW; // External inputs (root inport signals with auto storage) ExtU_EKF_T EKF_U; // External outputs (root outports fed by signals with auto storage) ExtY_EKF_T EKF_Y; // Real-time model RT_MODEL_EKF_T EKF_M_; RT_MODEL_EKF_T *const EKF_M = &EKF_M_; // Forward declaration for local functions static int16_T EKF_ixamax(int16_T n, const real_T x[13], int16_T ix0); static void EKF_xswap(real_T x[273], int16_T ix0, int16_T iy0); static real_T EKF_xnrm2_i(int16_T n, const real_T x[273], int16_T ix0); static real_T EKF_xzlarfg(int16_T n, real_T *alpha1, real_T x[273], int16_T ix0); static void EKF_xzlarf(int16_T m, int16_T n, int16_T iv0, real_T tau, real_T C [273], int16_T ic0, real_T work[13]); static real_T EKF_xnrm2_iz(int16_T n, const real_T x[273], int16_T ix0); static real_T EKF_xnrm2(const real_T x[273], int16_T ix0); static void EKF_xgeqp3(real_T A[273], real_T tau[13], int16_T jpvt[13]); static void EKF_qrsolve(const real_T A[273], const real_T B[21], real_T Y[13], int16_T *rankR); static void EKF_polyfit(const real_T x[21], const real_T y[21], real_T p[13]); static real_T EKF_xnrm2_izo(int16_T n, const real_T x[9], int16_T ix0); static real_T EKF_xnrm2_izoj(int16_T n, const real_T x[3], int16_T ix0); static void EKF_xaxpy_it(int16_T n, real_T a, const real_T x[3], int16_T ix0, real_T y[9], int16_T iy0); static void EKF_xaxpy_i(int16_T n, real_T a, const real_T x[9], int16_T ix0, real_T y[3], int16_T iy0); static void EKF_xrotg(real_T *a, real_T *b, real_T *c, real_T *s); static real_T EKF_xdotc(int16_T n, const real_T x[9], int16_T ix0, const real_T y[9], int16_T iy0); static void EKF_xaxpy(int16_T n, real_T a, int16_T ix0, real_T y[9], int16_T iy0); static void EKF_svd(const real_T A[9], real_T U[3]); static real_T EKF_norm(const real_T x[9]); // Function for MATLAB Function: '<Root>/EKF_SLX' static int16_T EKF_ixamax(int16_T n, const real_T x[13], int16_T ix0) { int16_T idxmax; int16_T ix; real_T smax; real_T s; int16_T k; if (n < 1) { idxmax = 0; } else { idxmax = 1; if (n > 1) { ix = ix0 - 1; smax = std::abs(x[ix0 - 1]); for (k = 2; k <= n; k++) { ix++; s = std::abs(x[ix]); if (s > smax) { idxmax = k; smax = s; } } } } return idxmax; } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xswap(real_T x[273], int16_T ix0, int16_T iy0) { int16_T ix; int16_T iy; real_T temp; int16_T k; ix = ix0 - 1; iy = iy0 - 1; for (k = 0; k < 21; k++) { temp = x[ix]; x[ix] = x[iy]; x[iy] = temp; ix++; iy++; } } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xnrm2_i(int16_T n, const real_T x[273], int16_T ix0) { real_T y; real_T scale; int16_T kend; real_T absxk; real_T t; int16_T k; y = 0.0; if (!(n < 1)) { if (n == 1) { y = std::abs(x[ix0 - 1]); } else { scale = 2.2250738585072014E-308; kend = (ix0 + n) - 1; for (k = ix0; k <= kend; k++) { absxk = std::abs(x[k - 1]); if (absxk > scale) { t = scale / absxk; y = y * t * t + 1.0; scale = absxk; } else { t = absxk / scale; y += t * t; } } y = scale * std::sqrt(y); } } return y; } real_T rt_hypotd_snf(real_T u0, real_T u1) { real_T y; real_T a; a = std::abs(u0); y = std::abs(u1); if (a < y) { a /= y; y *= std::sqrt(a * a + 1.0); } else if (a > y) { y /= a; y = std::sqrt(y * y + 1.0) * a; } else { if (!rtIsNaN(y)) { y = a * 1.4142135623730951; } } return y; } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xzlarfg(int16_T n, real_T *alpha1, real_T x[273], int16_T ix0) { real_T tau; real_T xnorm; int16_T knt; int16_T b_k; int16_T c_k; tau = 0.0; if (!(n <= 0)) { xnorm = EKF_xnrm2_i(n - 1, x, ix0); if (xnorm != 0.0) { xnorm = rt_hypotd_snf(*alpha1, xnorm); if (*alpha1 >= 0.0) { xnorm = -xnorm; } if (std::abs(xnorm) < 1.0020841800044864E-292) { knt = 0; do { knt++; b_k = (ix0 + n) - 2; for (c_k = ix0; c_k <= b_k; c_k++) { x[c_k - 1] *= 9.9792015476736E+291; } xnorm *= 9.9792015476736E+291; *alpha1 *= 9.9792015476736E+291; } while (!(std::abs(xnorm) >= 1.0020841800044864E-292)); xnorm = rt_hypotd_snf(*alpha1, EKF_xnrm2_i(n - 1, x, ix0)); if (*alpha1 >= 0.0) { xnorm = -xnorm; } tau = (xnorm - *alpha1) / xnorm; *alpha1 = 1.0 / (*alpha1 - xnorm); b_k = (ix0 + n) - 2; for (c_k = ix0; c_k <= b_k; c_k++) { x[c_k - 1] *= *alpha1; } for (b_k = 1; b_k <= knt; b_k++) { xnorm *= 1.0020841800044864E-292; } *alpha1 = xnorm; } else { tau = (xnorm - *alpha1) / xnorm; *alpha1 = 1.0 / (*alpha1 - xnorm); knt = (ix0 + n) - 2; for (b_k = ix0; b_k <= knt; b_k++) { x[b_k - 1] *= *alpha1; } *alpha1 = xnorm; } } } return tau; } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xzlarf(int16_T m, int16_T n, int16_T iv0, real_T tau, real_T C [273], int16_T ic0, real_T work[13]) { int16_T lastv; int16_T lastc; int16_T coltop; int16_T ix; real_T c; int16_T iac; int16_T d; int16_T b_ia; int16_T jy; int32_T exitg1; boolean_T exitg2; if (tau != 0.0) { lastv = m; lastc = iv0 + m; while ((lastv > 0) && (C[lastc - 2] == 0.0)) { lastv--; lastc--; } lastc = n; exitg2 = false; while ((!exitg2) && (lastc > 0)) { coltop = (lastc - 1) * 21 + ic0; jy = coltop; do { exitg1 = 0L; if (jy <= (coltop + lastv) - 1) { if (C[jy - 1] != 0.0) { exitg1 = 1L; } else { jy++; } } else { lastc--; exitg1 = 2L; } } while (exitg1 == 0L); if (exitg1 == 1L) { exitg2 = true; } } } else { lastv = 0; lastc = 0; } if (lastv > 0) { if (lastc != 0) { for (coltop = 1; coltop <= lastc; coltop++) { work[coltop - 1] = 0.0; } coltop = 0; jy = (lastc - 1) * 21 + ic0; for (iac = ic0; iac <= jy; iac += 21) { ix = iv0; c = 0.0; d = (iac + lastv) - 1; for (b_ia = iac; b_ia <= d; b_ia++) { c += C[b_ia - 1] * C[ix - 1]; ix++; } work[coltop] += c; coltop++; } } if (!(-tau == 0.0)) { coltop = ic0 - 1; jy = 0; for (iac = 1; iac <= lastc; iac++) { if (work[jy] != 0.0) { c = work[jy] * -tau; ix = iv0; d = lastv + coltop; for (b_ia = coltop; b_ia + 1 <= d; b_ia++) { C[b_ia] += C[ix - 1] * c; ix++; } } jy++; coltop += 21; } } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xnrm2_iz(int16_T n, const real_T x[273], int16_T ix0) { real_T y; real_T scale; int16_T kend; real_T absxk; real_T t; int16_T k; y = 0.0; if (!(n < 1)) { if (n == 1) { y = std::abs(x[ix0 - 1]); } else { scale = 2.2250738585072014E-308; kend = (ix0 + n) - 1; for (k = ix0; k <= kend; k++) { absxk = std::abs(x[k - 1]); if (absxk > scale) { t = scale / absxk; y = y * t * t + 1.0; scale = absxk; } else { t = absxk / scale; y += t * t; } } y = scale * std::sqrt(y); } } return y; } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xnrm2(const real_T x[273], int16_T ix0) { real_T y; real_T scale; real_T absxk; real_T t; int16_T k; y = 0.0; scale = 2.2250738585072014E-308; for (k = ix0; k <= ix0 + 20; k++) { absxk = std::abs(x[k - 1]); if (absxk > scale) { t = scale / absxk; y = y * t * t + 1.0; scale = absxk; } else { t = absxk / scale; y += t * t; } } return scale * std::sqrt(y); } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xgeqp3(real_T A[273], real_T tau[13], int16_T jpvt[13]) { real_T work[13]; real_T vn1[13]; real_T vn2[13]; int16_T k; int16_T pvt; int16_T itemp; real_T temp2; real_T b_atmp; int16_T i_i; k = 1; for (i_i = 0; i_i < 13; i_i++) { jpvt[i_i] = 1 + i_i; work[i_i] = 0.0; b_atmp = EKF_xnrm2(A, k); vn2[i_i] = b_atmp; k += 21; vn1[i_i] = b_atmp; } for (k = 0; k < 13; k++) { i_i = k * 21 + k; pvt = (EKF_ixamax(13 - k, vn1, k + 1) + k) - 1; if (pvt + 1 != k + 1) { EKF_xswap(A, 1 + 21 * pvt, 1 + 21 * k); itemp = jpvt[pvt]; jpvt[pvt] = jpvt[k]; jpvt[k] = itemp; vn1[pvt] = vn1[k]; vn2[pvt] = vn2[k]; } b_atmp = A[i_i]; temp2 = EKF_xzlarfg(21 - k, &b_atmp, A, i_i + 2); tau[k] = temp2; A[i_i] = b_atmp; if (k + 1 < 13) { b_atmp = A[i_i]; A[i_i] = 1.0; EKF_xzlarf(21 - k, 12 - k, i_i + 1, tau[k], A, (k + (k + 1) * 21) + 1, work); A[i_i] = b_atmp; } for (i_i = k + 1; i_i + 1 < 14; i_i++) { if (vn1[i_i] != 0.0) { b_atmp = std::abs(A[21 * i_i + k]) / vn1[i_i]; b_atmp = 1.0 - b_atmp * b_atmp; if (b_atmp < 0.0) { b_atmp = 0.0; } temp2 = vn1[i_i] / vn2[i_i]; temp2 = temp2 * temp2 * b_atmp; if (temp2 <= 1.4901161193847656E-8) { vn1[i_i] = EKF_xnrm2_iz(20 - k, A, (21 * i_i + k) + 2); vn2[i_i] = vn1[i_i]; } else { vn1[i_i] *= std::sqrt(b_atmp); } } } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_qrsolve(const real_T A[273], const real_T B[21], real_T Y[13], int16_T *rankR) { real_T b_A[273]; real_T tau[13]; int16_T jpvt[13]; real_T tol; real_T b_B[21]; int16_T b_i; int16_T b_j; memcpy(&b_A[0], &A[0], 273U * sizeof(real_T)); EKF_xgeqp3(b_A, tau, jpvt); *rankR = 0; tol = 21.0 * std::abs(b_A[0]) * 2.2204460492503131E-16; while ((*rankR < 13) && (std::abs(b_A[21 * *rankR + *rankR]) >= tol)) { (*rankR)++; } memset(&Y[0], 0, 13U * sizeof(real_T)); memcpy(&b_B[0], &B[0], 21U * sizeof(real_T)); for (b_j = 0; b_j < 13; b_j++) { if (tau[b_j] != 0.0) { tol = b_B[b_j]; for (b_i = b_j + 1; b_i + 1 < 22; b_i++) { tol += b_A[21 * b_j + b_i] * b_B[b_i]; } tol *= tau[b_j]; if (tol != 0.0) { b_B[b_j] -= tol; for (b_i = b_j + 1; b_i + 1 < 22; b_i++) { b_B[b_i] -= b_A[21 * b_j + b_i] * tol; } } } } for (b_j = 0; b_j < 13; b_j++) { Y[jpvt[b_j] - 1] = b_B[b_j]; } for (b_j = 12; b_j >= 0; b_j += -1) { Y[jpvt[b_j] - 1] /= b_A[21 * b_j + b_j]; for (b_i = 0; b_i + 1 <= b_j; b_i++) { Y[jpvt[b_i] - 1] -= b_A[21 * b_j + b_i] * Y[jpvt[b_j] - 1]; } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_polyfit(const real_T x[21], const real_T y[21], real_T p[13]) { real_T V[273]; real_T p1[13]; int16_T k; int16_T c_k; for (k = 0; k < 21; k++) { V[252 + k] = 1.0; V[231 + k] = x[k]; } for (k = 0; k < 11; k++) { for (c_k = 0; c_k < 21; c_k++) { V[c_k + 21 * (10 - k)] = V[(11 - k) * 21 + c_k] * x[c_k]; } } EKF_qrsolve(V, y, p1, &k); memcpy(&p[0], &p1[0], 13U * sizeof(real_T)); } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xnrm2_izo(int16_T n, const real_T x[9], int16_T ix0) { real_T y; real_T scale; int16_T kend; real_T absxk; real_T t; int16_T k; y = 0.0; if (!(n < 1)) { if (n == 1) { y = std::abs(x[ix0 - 1]); } else { scale = 2.2250738585072014E-308; kend = (ix0 + n) - 1; for (k = ix0; k <= kend; k++) { absxk = std::abs(x[k - 1]); if (absxk > scale) { t = scale / absxk; y = y * t * t + 1.0; scale = absxk; } else { t = absxk / scale; y += t * t; } } y = scale * std::sqrt(y); } } return y; } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xnrm2_izoj(int16_T n, const real_T x[3], int16_T ix0) { real_T y; real_T scale; int16_T kend; real_T absxk; real_T t; int16_T k; y = 0.0; if (!(n < 1)) { if (n == 1) { y = std::abs(x[ix0 - 1]); } else { scale = 2.2250738585072014E-308; kend = (ix0 + n) - 1; for (k = ix0; k <= kend; k++) { absxk = std::abs(x[k - 1]); if (absxk > scale) { t = scale / absxk; y = y * t * t + 1.0; scale = absxk; } else { t = absxk / scale; y += t * t; } } y = scale * std::sqrt(y); } } return y; } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xaxpy_it(int16_T n, real_T a, const real_T x[3], int16_T ix0, real_T y[9], int16_T iy0) { int16_T ix; int16_T iy; int16_T k; if (!((n < 1) || (a == 0.0))) { ix = ix0 - 1; iy = iy0 - 1; for (k = 0; k < n; k++) { y[iy] += a * x[ix]; ix++; iy++; } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xaxpy_i(int16_T n, real_T a, const real_T x[9], int16_T ix0, real_T y[3], int16_T iy0) { int16_T ix; int16_T iy; int16_T k; if (!((n < 1) || (a == 0.0))) { ix = ix0 - 1; iy = iy0 - 1; for (k = 0; k < n; k++) { y[iy] += a * x[ix]; ix++; iy++; } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xrotg(real_T *a, real_T *b, real_T *c, real_T *s) { real_T roe; real_T absa; real_T absb; real_T scale; real_T ads; real_T bds; roe = *b; absa = std::abs(*a); absb = std::abs(*b); if (absa > absb) { roe = *a; } scale = absa + absb; if (scale == 0.0) { *s = 0.0; *c = 1.0; scale = 0.0; absa = 0.0; } else { ads = absa / scale; bds = absb / scale; scale *= std::sqrt(ads * ads + bds * bds); if (roe < 0.0) { scale = -scale; } *c = *a / scale; *s = *b / scale; if (absa > absb) { absa = *s; } else if (*c != 0.0) { absa = 1.0 / *c; } else { absa = 1.0; } } *a = scale; *b = absa; } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_xdotc(int16_T n, const real_T x[9], int16_T ix0, const real_T y[9], int16_T iy0) { real_T d; int16_T ix; int16_T iy; int16_T k; d = 0.0; if (!(n < 1)) { ix = ix0; iy = iy0; for (k = 1; k <= n; k++) { d += x[ix - 1] * y[iy - 1]; ix++; iy++; } } return d; } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_xaxpy(int16_T n, real_T a, int16_T ix0, real_T y[9], int16_T iy0) { int16_T ix; int16_T iy; int16_T k; if (!((n < 1) || (a == 0.0))) { ix = ix0 - 1; iy = iy0 - 1; for (k = 0; k < n; k++) { y[iy] += a * y[ix]; ix++; iy++; } } } // Function for MATLAB Function: '<Root>/EKF_SLX' static void EKF_svd(const real_T A[9], real_T U[3]) { real_T b_A[9]; real_T s[3]; real_T e[3]; real_T work[3]; int16_T q; boolean_T apply_transform; int16_T m; int16_T iter; real_T snorm; real_T ztest0; int16_T kase; int16_T qs; real_T ztest; real_T smm1; real_T emm1; real_T sqds; real_T shift; int16_T d_ii; real_T varargin_1[5]; boolean_T exitg1; boolean_T exitg2; int32_T exitg3; memcpy(&b_A[0], &A[0], 9U * sizeof(real_T)); e[0] = 0.0; work[0] = 0.0; e[1] = 0.0; work[1] = 0.0; e[2] = 0.0; work[2] = 0.0; apply_transform = false; snorm = EKF_xnrm2_izo(3, b_A, 1); if (snorm > 0.0) { apply_transform = true; if (b_A[0] < 0.0) { s[0] = -snorm; } else { s[0] = snorm; } if (std::abs(s[0]) >= 1.0020841800044864E-292) { snorm = 1.0 / s[0]; for (qs = 0; qs + 1 < 4; qs++) { b_A[qs] *= snorm; } } else { for (qs = 0; qs + 1 < 4; qs++) { b_A[qs] /= s[0]; } } b_A[0]++; s[0] = -s[0]; } else { s[0] = 0.0; } for (q = 1; q + 1 < 4; q++) { qs = 3 * q; if (apply_transform) { EKF_xaxpy(3, -(EKF_xdotc(3, b_A, 1, b_A, qs + 1) / b_A[0]), 1, b_A, qs + 1); } e[q] = b_A[qs]; } snorm = EKF_xnrm2_izoj(2, e, 2); if (snorm == 0.0) { e[0] = 0.0; } else { if (e[1] < 0.0) { snorm = -snorm; } e[0] = snorm; if (std::abs(snorm) >= 1.0020841800044864E-292) { snorm = 1.0 / snorm; for (iter = 1; iter + 1 < 4; iter++) { e[iter] *= snorm; } } else { for (iter = 1; iter + 1 < 4; iter++) { e[iter] /= snorm; } } e[1]++; e[0] = -e[0]; for (iter = 2; iter < 4; iter++) { work[iter - 1] = 0.0; } for (iter = 1; iter + 1 < 4; iter++) { EKF_xaxpy_i(2, e[iter], b_A, 3 * iter + 2, work, 2); } for (iter = 1; iter + 1 < 4; iter++) { EKF_xaxpy_it(2, -e[iter] / e[1], work, 2, b_A, 3 * iter + 2); } } apply_transform = false; snorm = EKF_xnrm2_izo(2, b_A, 5); if (snorm > 0.0) { apply_transform = true; if (b_A[4] < 0.0) { s[1] = -snorm; } else { s[1] = snorm; } if (std::abs(s[1]) >= 1.0020841800044864E-292) { snorm = 1.0 / s[1]; for (qs = 4; qs + 1 < 7; qs++) { b_A[qs] *= snorm; } } else { for (qs = 4; qs + 1 < 7; qs++) { b_A[qs] /= s[1]; } } b_A[4]++; s[1] = -s[1]; } else { s[1] = 0.0; } for (q = 2; q + 1 < 4; q++) { qs = 3 * q + 1; if (apply_transform) { EKF_xaxpy(2, -(EKF_xdotc(2, b_A, 5, b_A, qs + 1) / b_A[4]), 5, b_A, qs + 1); } e[q] = b_A[qs]; } m = 1; s[2] = b_A[8]; ztest = e[0]; if (s[0] != 0.0) { ztest0 = std::abs(s[0]); snorm = s[0] / ztest0; s[0] = ztest0; ztest = e[0] / snorm; } if (ztest != 0.0) { ztest0 = std::abs(ztest); snorm = ztest; ztest = ztest0; s[1] *= ztest0 / snorm; } e[0] = ztest; ztest = b_A[7]; if (s[1] != 0.0) { ztest0 = std::abs(s[1]); snorm = s[1] / ztest0; s[1] = ztest0; ztest = b_A[7] / snorm; } if (ztest != 0.0) { ztest0 = std::abs(ztest); snorm = ztest; ztest = ztest0; s[2] = ztest0 / snorm * b_A[8]; } e[1] = ztest; if (s[2] != 0.0) { s[2] = std::abs(s[2]); } e[2] = 0.0; iter = 0; snorm = 0.0; if ((s[0] >= e[0]) || rtIsNaN(e[0])) { ztest0 = s[0]; } else { ztest0 = e[0]; } if (!((0.0 >= ztest0) || rtIsNaN(ztest0))) { snorm = ztest0; } if ((s[1] >= ztest) || rtIsNaN(ztest)) { ztest = s[1]; } if (!((snorm >= ztest) || rtIsNaN(ztest))) { snorm = ztest; } if (s[2] >= 0.0) { ztest0 = s[2]; } else { ztest0 = 0.0; } if (!((snorm >= ztest0) || rtIsNaN(ztest0))) { snorm = ztest0; } while ((m + 2 > 0) && (!(iter >= 75))) { kase = m + 1; do { exitg3 = 0L; q = kase; if (kase == 0) { exitg3 = 1L; } else { ztest0 = std::abs(e[kase - 1]); if ((ztest0 <= (std::abs(s[kase - 1]) + std::abs(s[kase])) * 2.2204460492503131E-16) || ((ztest0 <= 1.0020841800044864E-292) || ((iter > 20) && (ztest0 <= 2.2204460492503131E-16 * snorm)))) { e[kase - 1] = 0.0; exitg3 = 1L; } else { kase--; } } } while (exitg3 == 0L); if (m + 1 == kase) { kase = 4; } else { qs = m + 2; d_ii = m + 2; exitg2 = false; while ((!exitg2) && (d_ii >= kase)) { qs = d_ii; if (d_ii == kase) { exitg2 = true; } else { ztest0 = 0.0; if (d_ii < m + 2) { ztest0 = std::abs(e[d_ii - 1]); } if (d_ii > kase + 1) { ztest0 += std::abs(e[d_ii - 2]); } ztest = std::abs(s[d_ii - 1]); if ((ztest <= 2.2204460492503131E-16 * ztest0) || (ztest <= 1.0020841800044864E-292)) { s[d_ii - 1] = 0.0; exitg2 = true; } else { d_ii--; } } } if (qs == kase) { kase = 3; } else if (m + 2 == qs) { kase = 1; } else { kase = 2; q = qs; } } switch (kase) { case 1: ztest0 = e[m]; e[m] = 0.0; for (qs = m; qs + 1 >= q + 1; qs--) { ztest = s[qs]; EKF_xrotg(&ztest, &ztest0, &sqds, &smm1); s[qs] = ztest; if (qs + 1 > q + 1) { ztest0 = -smm1 * e[0]; e[0] *= sqds; } } break; case 2: ztest0 = e[q - 1]; e[q - 1] = 0.0; while (q + 1 <= m + 2) { ztest = s[q]; EKF_xrotg(&ztest, &ztest0, &sqds, &smm1); s[q] = ztest; ztest0 = -smm1 * e[q]; e[q] *= sqds; q++; } break; case 3: varargin_1[0] = std::abs(s[m + 1]); varargin_1[1] = std::abs(s[m]); varargin_1[2] = std::abs(e[m]); varargin_1[3] = std::abs(s[q]); varargin_1[4] = std::abs(e[q]); qs = 1; ztest = varargin_1[0]; if (rtIsNaN(varargin_1[0])) { kase = 2; exitg1 = false; while ((!exitg1) && (kase < 6)) { qs = kase; if (!rtIsNaN(varargin_1[kase - 1])) { ztest = varargin_1[kase - 1]; exitg1 = true; } else { kase++; } } } if (qs < 5) { while (qs + 1 < 6) { if (varargin_1[qs] > ztest) { ztest = varargin_1[qs]; } qs++; } } ztest0 = s[m + 1] / ztest; smm1 = s[m] / ztest; emm1 = e[m] / ztest; sqds = s[q] / ztest; smm1 = ((smm1 + ztest0) * (smm1 - ztest0) + emm1 * emm1) / 2.0; emm1 *= ztest0; emm1 *= emm1; if ((smm1 != 0.0) || (emm1 != 0.0)) { shift = std::sqrt(smm1 * smm1 + emm1); if (smm1 < 0.0) { shift = -shift; } shift = emm1 / (smm1 + shift); } else { shift = 0.0; } ztest0 = (sqds + ztest0) * (sqds - ztest0) + shift; ztest = e[q] / ztest * sqds; for (qs = q + 1; qs <= m + 1; qs++) { EKF_xrotg(&ztest0, &ztest, &sqds, &smm1); if (qs > q + 1) { e[0] = ztest0; } ztest0 = s[qs - 1] * sqds + e[qs - 1] * smm1; e[qs - 1] = e[qs - 1] * sqds - s[qs - 1] * smm1; ztest = smm1 * s[qs]; s[qs] *= sqds; EKF_xrotg(&ztest0, &ztest, &sqds, &smm1); s[qs - 1] = ztest0; ztest0 = e[qs - 1] * sqds + smm1 * s[qs]; s[qs] = e[qs - 1] * -smm1 + sqds * s[qs]; ztest = smm1 * e[qs]; e[qs] *= sqds; } e[m] = ztest0; iter++; break; default: if (s[q] < 0.0) { s[q] = -s[q]; } iter = q + 1; while ((q + 1 < 3) && (s[q] < s[iter])) { ztest0 = s[q]; s[q] = s[iter]; s[iter] = ztest0; q = iter; iter++; } iter = 0; m--; break; } } U[0] = s[0]; U[1] = s[1]; U[2] = s[2]; } // Function for MATLAB Function: '<Root>/EKF_SLX' static real_T EKF_norm(const real_T x[9]) { real_T y; real_T s[3]; real_T absx; int16_T k; boolean_T exitg1; y = 0.0; k = 0; exitg1 = false; while ((!exitg1) && (k < 9)) { absx = std::abs(x[k]); if (rtIsNaN(absx)) { y = (rtNaN); exitg1 = true; } else { if (absx > y) { y = absx; } k++; } } if ((!rtIsInf(y)) && (!rtIsNaN(y))) { EKF_svd(x, s); y = s[0]; } return y; } // Model step function void EKF_step(void) { real_T p[13]; real_T C[3]; real_T KGain[3]; real_T OCV; real_T dOCV[21]; real_T work; real_T b_y1[20]; int16_T ixLead; int16_T iyLead; real_T tmp2; int16_T high_i; int16_T mid_i; static const real_T f[21] = { 2.1301, 2.38, 2.4114, 2.4468, 2.4765, 2.4991, 2.5282999999999998, 2.5524, 2.5675, 2.5841, 2.6041, 2.6292999999999997, 2.6601, 2.7028, 2.7662, 2.8141, 2.8638, 2.9168, 2.9734000000000003, 3.0340000000000003, 3.1 }; static const real_T b_y[9] = { 3.5707572690619888E-21, -3.8357517594746764E-22, -1.8896447467876043E-13, -3.8357517594746764E-22, 4.1204121287633946E-23, 2.0298798311139986E-14, -1.8896447467876043E-13, 2.0298798311139986E-14, 1.0E-5 }; static const real_T OCV_data[21] = { 2.1301, 2.38, 2.4114, 2.4468, 2.4765, 2.4991, 2.5282999999999998, 2.5524, 2.5675, 2.5841, 2.6041, 2.6292999999999997, 2.6601, 2.7028, 2.7662, 2.8141, 2.8638, 2.9168, 2.9734000000000003, 3.0340000000000003, 3.1 }; static const real_T SOC_spacing[21] = { 0.0, 0.05, 0.1, 0.15000000000000002, 0.2, 0.25, 0.30000000000000004, 0.35000000000000003, 0.4, 0.45, 0.5, 0.55, 0.6, 0.64999999999999991, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0 }; static const real_T g[21] = { 0.0, 0.05, 0.1, 0.15000000000000002, 0.2, 0.25, 0.30000000000000004, 0.35000000000000003, 0.4, 0.45, 0.5, 0.55, 0.6, 0.64999999999999991, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0 }; static const real_T a[9] = { 1.0, 0.0, 0.0, 0.0, 0.9999985539602041, 0.0, 0.0, 0.0, 1.0 }; static const real_T b_a[3] = { -1.8896447467876042E-8, 2.0298798311139983E-9, 0.0 }; real_T a_0[9]; real_T b_y1_0[21]; real_T b_y1_1[21]; real_T KGain_0; // MATLAB Function: '<Root>/EKF_SLX' incorporates: // Inport: '<Root>/Current' // Inport: '<Root>/Voltage' // MATLAB Function 'EKF_SLX': '<S1>:1' // '<S1>:1:3' coder.extrinsic('display') // CONVERT VOLTAGE_CURRENT MEASURMENT TO ACTUAL CURRENT // dSPACE_range=1; %[V] // I_range=1; %[A] // Current=V_Current*I_range/dSPACE_range; // current through resistor 100kOhm // '<S1>:1:12' dt=100e-6; // Q_total_t=[190800]; %Ampere*seconds // '<S1>:1:15' Q_total_t=3600*1470e-3; // R_0_t=0.0011622500; %Ohm // '<S1>:1:18' R1_t=0.00140375000; // Ohm // '<S1>:1:19' C1_t=49264; // '<S1>:1:20' OCV_data=[3.2141 // '<S1>:1:21' 3.4640 // '<S1>:1:22' 3.4954 // '<S1>:1:23' 3.5308 // '<S1>:1:24' 3.5605 // '<S1>:1:25' 3.5831 // '<S1>:1:26' 3.6123 // '<S1>:1:27' 3.6364 // '<S1>:1:28' 3.6515 // '<S1>:1:29' 3.6681 // '<S1>:1:30' 3.6881 // '<S1>:1:31' 3.7133 // '<S1>:1:32' 3.7441 // '<S1>:1:33' 3.7868 // '<S1>:1:34' 3.8502 // '<S1>:1:35' 3.8981 // '<S1>:1:36' 3.9478 // '<S1>:1:37' 4.0008 // '<S1>:1:38' 4.0574 // '<S1>:1:39' 4.1180 // '<S1>:1:40' 4.1840]'; // '<S1>:1:42' OCV_data=OCV_data-OCV_data(numel(OCV_data))+3.1; // '<S1>:1:45' z=0:0.05:1; // '<S1>:1:46' Order=12; // EKF ALLOWS FOR A NON LINEAR OCV FUNCTION // '<S1>:1:47' [p]=polyfit(z,OCV_data,Order); EKF_polyfit(SOC_spacing, OCV_data, p); // '<S1>:1:48' Q_total=Q_total_t; // '<S1>:1:49' R1=R1_t; // '<S1>:1:50' C1=C1_t; // '<S1>:1:51' t_RC=R1*C1; // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // INITIALISE SS VECTORS AND COV. MATRICES // '<S1>:1:60' if isempty(Current_prev) // '<S1>:1:65' if isempty(Xhat) // '<S1>:1:71' if isempty(P0) // '<S1>:1:79' sigmaW=1e-5; // COV. OF PROCESS NOISE (CURRENT) [Q] // '<S1>:1:81' SigmaV=1e-3; // COV. OF SENSOR NOISE (VOLTAGE) [R] // EXTENDED KALMA FILTER MATRICES // '<S1>:1:86' Ahat=[1 0 0 // '<S1>:1:87' 0 exp(-dt/t_RC) 0 // '<S1>:1:88' 0 0 1]; // '<S1>:1:90' Bhat_I=[-dt/Q_total // '<S1>:1:91' dt/C1 // '<S1>:1:92' 0]; // '<S1>:1:94' Bhat_W=[-dt/Q_total // '<S1>:1:95' dt/C1 // '<S1>:1:96' 1]; // '<S1>:1:99' Xhat=Ahat*Xhat+Bhat_I*Current_prev; for (ixLead = 0; ixLead < 3; ixLead++) { KGain[ixLead] = ((a[ixLead + 3] * EKF_DW.Xhat[1] + a[ixLead] * EKF_DW.Xhat[0]) + a[ixLead + 6] * EKF_DW.Xhat[2]) + b_a[ixLead] * EKF_DW.Current_prev; } // STATE PREDICTION TIME UPDATE // '<S1>:1:101' P0= Ahat*P0*Ahat'+ Bhat_W*sigmaW*Bhat_W'; for (ixLead = 0; ixLead < 3; ixLead++) { EKF_DW.Xhat[ixLead] = KGain[ixLead]; for (iyLead = 0; iyLead < 3; iyLead++) { a_0[ixLead + 3 * iyLead] = 0.0; a_0[ixLead + 3 * iyLead] += EKF_DW.P0[3 * iyLead] * a[ixLead]; a_0[ixLead + 3 * iyLead] += EKF_DW.P0[3 * iyLead + 1] * a[ixLead + 3]; a_0[ixLead + 3 * iyLead] += EKF_DW.P0[3 * iyLead + 2] * a[ixLead + 6]; } } for (ixLead = 0; ixLead < 3; ixLead++) { for (iyLead = 0; iyLead < 3; iyLead++) { EKF_DW.P0[iyLead + 3 * ixLead] = ((a[3 * ixLead + 1] * a_0[iyLead + 3] + a[3 * ixLead] * a_0[iyLead]) + a[3 * ixLead + 2] * a_0[iyLead + 6]) + b_y[3 * ixLead + iyLead]; } } // ERROR COV. TIME UPDATE // '<S1>:1:105' Ytrue=Voltage; // GET TRUE TERMINAL VOLTAGE READING // '<S1>:1:107' Yhat=OCV(p,Xhat(1))-Xhat(2)-Xhat(3)*Current; // '<S1>:1:151' OCV=polyval(p,SOC); OCV = p[0]; for (ixLead = 0; ixLead < 12; ixLead++) { OCV = EKF_DW.Xhat[0] * OCV + p[ixLead + 1]; } // ESTIMATE SYSTEM OUTPUT USING NON LINEAR OCV(Z) // '<S1>:1:109' C=[dOCVatZ(z,OCV_data,Xhat(1)) -1 -Current]; // FIND THE NUMERICAL DERIVATIVE AT A SPECIFIC POINT BY INTERPOLATING // '<S1>:1:145' dOCV=dOCV_Num(SOC_spacing,OCV_data); // THIS FUNCTION RETURNS THE NUMERICAL DERIVATIVE OF SOME DATA Y (OCV DATA) // AND ITS SPACING Z // '<S1>:1:158' OCV=OCV_data; // '<S1>:1:159' dz=z(2)-z(1); // '<S1>:1:160' dUdZ=diff(OCV)/dz; ixLead = 1; iyLead = 0; work = 2.1301; for (high_i = 0; high_i < 20; high_i++) { tmp2 = work; work = f[ixLead]; tmp2 = f[ixLead] - tmp2; ixLead++; b_y1[iyLead] = tmp2; iyLead++; } // FIND THE NUMERICAL DERIVATIVE // '<S1>:1:162' dOCV_data=([dUdZ(1) dUdZ]+ [dUdZ dUdZ(end)])/2; for (ixLead = 0; ixLead < 20; ixLead++) { work = b_y1[ixLead] / 0.05; b_y1_0[ixLead + 1] = work; b_y1_1[ixLead] = work; b_y1[ixLead] = work; } b_y1_0[0] = b_y1[0]; b_y1_1[20] = b_y1[19]; for (ixLead = 0; ixLead < 21; ixLead++) { dOCV[ixLead] = (b_y1_0[ixLead] + b_y1_1[ixLead]) / 2.0; } // GET THE WHOLE VECTOR OF THE DERIVATIVE // '<S1>:1:147' dOCVz=interp1(SOC_spacing,dOCV,z); work = (rtNaN); if ((!rtIsNaN(EKF_DW.Xhat[0])) && (!(EKF_DW.Xhat[0] > 1.0)) && (!(EKF_DW.Xhat [0] < 0.0))) { ixLead = 0; iyLead = 2; high_i = 21; while (high_i > iyLead) { mid_i = ((ixLead + high_i) + 1) >> 1; if (EKF_DW.Xhat[0] >= g[mid_i - 1]) { ixLead = mid_i - 1; iyLead = mid_i + 1; } else { high_i = mid_i; } } work = (EKF_DW.Xhat[0] - g[ixLead]) / (g[ixLead + 1] - g[ixLead]); if (work == 0.0) { work = dOCV[ixLead]; } else if (work == 1.0) { work = dOCV[ixLead + 1]; } else if (dOCV[ixLead + 1] == dOCV[ixLead]) { work = dOCV[ixLead]; } else { work = (1.0 - work) * dOCV[ixLead] + dOCV[ixLead + 1] * work; } } C[0] = work; C[1] = -1.0; C[2] = -EKF_U.Current; // FIND C MATRIX FOR GAIN CALCULATION // '<S1>:1:110' KGain=P0*C'*(C*P0*C'+SigmaV)^-1; tmp2 = 0.0; for (ixLead = 0; ixLead < 3; ixLead++) { tmp2 += (EKF_DW.P0[3 * ixLead + 2] * -EKF_U.Current + (-EKF_DW.P0[3 * ixLead + 1] + EKF_DW.P0[3 * ixLead] * work)) * C[ixLead]; } work = 1.0 / (tmp2 + 0.001); // KALMAN FILTER GAIN CALCULATION // STATE ESTIMATE MEASURMENT UPDATE // '<S1>:1:113' Xhat=Xhat+KGain*(Ytrue - Yhat); OCV = EKF_U.Voltage - ((OCV - EKF_DW.Xhat[1]) - EKF_DW.Xhat[2] * EKF_U.Current); // '<S1>:1:115' P0=P0-KGain*(C*P0*C'+SigmaV)*KGain'; tmp2 = 0.0; for (ixLead = 0; ixLead < 3; ixLead++) { KGain_0 = (EKF_DW.P0[ixLead + 6] * -EKF_U.Current + (-EKF_DW.P0[ixLead + 3] + EKF_DW.P0[ixLead] * C[0])) * work; EKF_DW.Xhat[ixLead] += KGain_0 * OCV; tmp2 += (EKF_DW.P0[3 * ixLead + 2] * -EKF_U.Current + (-EKF_DW.P0[3 * ixLead + 1] + EKF_DW.P0[3 * ixLead] * C[0])) * C[ixLead]; KGain[ixLead] = KGain_0; } for (ixLead = 0; ixLead < 3; ixLead++) { EKF_DW.P0[ixLead] -= (tmp2 + 0.001) * KGain[ixLead] * KGain[0]; EKF_DW.P0[ixLead + 3] -= (tmp2 + 0.001) * KGain[ixLead] * KGain[1]; EKF_DW.P0[ixLead + 6] -= (tmp2 + 0.001) * KGain[ixLead] * KGain[2]; } // '<S1>:1:118' Current_prev=Current; EKF_DW.Current_prev = EKF_U.Current; // MAKE SURE THERE IS NO NEGATIVE SOC OUTPUT // '<S1>:1:121' if Xhat(1)<0 if (EKF_DW.Xhat[0] < 0.0) { // '<S1>:1:122' Xhat(1)=0; EKF_DW.Xhat[0] = 0.0; } // '<S1>:1:125' if Xhat(1)>1 if (EKF_DW.Xhat[0] > 1.0) { // '<S1>:1:126' Xhat(1)=1; EKF_DW.Xhat[0] = 1.0; } // '<S1>:1:129' if isnan(norm(P0)) OCV = EKF_norm(EKF_DW.P0); if (rtIsNaN(OCV)) { // '<S1>:1:130' P0=zeros(3,3)+diag([0.5, 0.05, 0.5]); memset(&EKF_DW.P0[0], 0, 9U * sizeof(real_T)); EKF_DW.P0[0] = 0.5; EKF_DW.P0[4] = 0.05; EKF_DW.P0[8] = 0.5; // '<S1>:1:131' Xhat=[1 0.05 1]'; EKF_DW.Xhat[0] = 1.0; EKF_DW.Xhat[1] = 0.05; EKF_DW.Xhat[2] = 1.0; } // '<S1>:1:134' if isinf(norm(P0)) OCV = EKF_norm(EKF_DW.P0); if (rtIsInf(OCV)) { // '<S1>:1:135' P0=zeros(3,3)+diag([0.5, 0.05, 0.5]); memset(&EKF_DW.P0[0], 0, 9U * sizeof(real_T)); EKF_DW.P0[0] = 0.5; EKF_DW.P0[4] = 0.05; EKF_DW.P0[8] = 0.5; // '<S1>:1:136' Xhat=[1 0.05 1]'; EKF_DW.Xhat[0] = 1.0; EKF_DW.Xhat[1] = 0.05; EKF_DW.Xhat[2] = 1.0; } // Outport: '<Root>/SOC' incorporates: // MATLAB Function: '<Root>/EKF_SLX' // '<S1>:1:140' R0=Xhat(3); // '<S1>:1:141' SOC=Xhat(1); EKF_Y.SOC = EKF_DW.Xhat[0]; // Outport: '<Root>/R0' incorporates: // MATLAB Function: '<Root>/EKF_SLX' EKF_Y.R0 = EKF_DW.Xhat[2]; } // Model initialize function void EKF_initialize(void) { // Registration code // initialize non-finites rt_InitInfAndNaN(sizeof(real_T)); // initialize error status rtmSetErrorStatus(EKF_M, (NULL)); // states (dwork) (void) memset((void *)&EKF_DW, 0, sizeof(DW_EKF_T)); // external inputs (void)memset((void *)&EKF_U, 0, sizeof(ExtU_EKF_T)); // external outputs (void) memset((void *)&EKF_Y, 0, sizeof(ExtY_EKF_T)); { static const real_T tmp[9] = { 0.1, 0.0, 0.0, 0.0, 0.05, 0.0, 0.0, 0.0, 0.5 }; // SystemInitialize for MATLAB Function: '<Root>/EKF_SLX' EKF_DW.Xhat[0] = 0.5; EKF_DW.Xhat[1] = 0.1; EKF_DW.Xhat[2] = 1.0; memcpy(&EKF_DW.P0[0], &tmp[0], 9U * sizeof(real_T)); // '<S1>:1:61' Current_prev=0; EKF_DW.Current_prev = 0.0; // '<S1>:1:66' Xhat=[0.5 // '<S1>:1:67' 0.1 // '<S1>:1:68' 1]; // '<S1>:1:72' P0=[.1 0 0 // '<S1>:1:73' 0 0.05 0 // '<S1>:1:74' 0 0 0.5]; } } // Model terminate function void EKF_terminate(void) { // (no terminate code required) } // // File trailer for generated code. // // [EOF] //