接下來的Matlab數值計算，將以電池建模的foster型等效電路為基礎，在連續時間內對電池進行系統同定。

foster等效電路

首先，我們利用輸入輸出的數據來推測電池模型的參數。推測的參數為：

θ = [ R0 Rd Cd FCC SOC0 ] T

輸入輸出數據與SOC

電池的建模（Matlab)

數值實驗用輸入出數據的生成：

[ex_batt.m]1

%%假想實驗數據的作成

Te = 2*3600; %實驗結束時刻[sec]

Ts = 1; %采樣周期[sec]

t = ( 0: Ts: Te );

N = length(t);

%輸入（電流）的設定

u = 40sawtooth( tsqrt(2)) + 10*sin(t) - 18;

%輸出（電壓）的生成

SOC0 = 95; %初期SOC[%]

FCC = 40*3600; %滿充電容量[C]

R0 = 0.450e-3; %電池內阻[ohm]

Rd = 0.500e-3; %擴散電阻[ohm]

Cd = 82000; %擴散電容[F]

th0 = [ R0, Rd, Cd, FCC ];

Nd = 3; %foster等效電路的近似級數

[ f, h, A, C, ymodel ] = batterymodel_foster(Nd);

%Simulation

[ y0, x ] = ymodel( u, t, th0, [SOC0; zeros(Nd, 1)] );

SOC = x(1, :);

%含有誤差的采樣值的模擬

um = u + randn(1, N) * 0.1; %電流傳感器的誤差

ym = y0 + randn(1, N) * 0.01; %電壓傳感器的誤差

%波形

figure(3), hold on

subplot(3,1,1); plot(t/60, u)

ylabel('Current [A]'), xlim([0 Te/60]), ylim([-100 50)]

subplot(3,1,2); plot(t/60, y0)

ylabel('Voltage [V]'), xlim([0 Te/60]), ylim([3.2 4.2)]

subplot(3,1,3); plot(t/60, SOC)

ylabel('SOC [%]'), xlim([0 Te/60]), ylim([0 100)]

xlabel('Time [min])，

上述 “m file” 里面使用的函數 “batterymodel_foster" ， 用以下方式記述。

[batterymodel_foster.m]

function [ f, h, A, C, ymodel ] = batterymodel_foster(Nd)

n = 1:Nd;

％參數矢量與物理常數的關系定義

R0 = @(th) th(1);

Rd = @(th) th(2);

Cd = @(th) th(3);

FCC = @(th) th(4);



A = @th blkdiag(0, -pi^2/4 * diag((2*n-1).^2)/Rd(th)/Cd(th();

B = @th [100/FCC(th); 2*ones(Nd,1)/Cd(th)];



f = @(x, u, th) A(th)*x + B(th)*u;

h = @(x, u, th) SOC2OCV(x(1,:)) + [0, ones(1, Nd)]*x+R0(th)*u;



%利用數值微分來計算SOC的雅可比數列

C = @(x) [numdiff(@SOC2OCV, x(1)), ones(1, Nd)];

%yhat:預測輸出的函數

function [y, x] = proto_ymode(u, t, th, x0)

    x = lsim(ss(A(th), B(th), eye(Nd+1), 0), u, t, x0, 'zoh')';

    y = h(x, u, th);

end

ymodel = @proto_ymodel;

end

上述的雅可比數列A，C為狀態推算時使用。

我們再來看看偏微分的數值計算。

[numdif.m]

function [dfdx] = numdiff(f, x)

n = length(x);

h = eye(n)*1e-5;

dfdx = arrayfun(@(k) (f(x+h(:, k)) - f(x-h(:,k)))/(2*h(k, k)), 1:n);

end

另外，SOC-OCV特性的簡易模型我們可以用以下m file來記述。

[SOC2OCV.m]

function ocv=SOC2OCV(SOC)

%系數設定

E0 = 4.14;    

K1 = 0.237;

K2 = -0.0516;

K3 = 1.05e-3;

K4 = 0.183;



mode = @(soc) E0 + K1.*log(soc/100) + K2.*log(1-soc/100) - K3./soc*100 - K4.*soc/100;

ocv = model(SOC);



ocv(SOC >98) = numdiff(model, 98)*(SOC(SOC >98)-98) + model(98);

ocv(SOC 2) = numdiff(model, 2)*(SOC(SOC< 2)-2) + model(2);

end

下面繼續［ex_batt.m]來進行連續時間里的電池系統同定，轉而來推算模型的參數。

[ex_batt.m]2

％％連續時間系統同定

yhat = @(th) ymodel(um, t, th(1:4), [th(5); zeros(Nd, 1)]);

%初期推算值

thhat0 = [1e-3, 1e-3, 1e5, 30*3600, 80];

%輸出誤差最小化的參數推算

thhat = Isqnonlin(...

@(th) ym-yhat(th.*thhat0), thhat0./thhat0, [ ], [ ],...

            optimoptions(@Isqnonlin, 'Display', 'iter',...

            'Algorithm', 'levenberg-marquardt')) .*thhat0

綜如上述所訴，在連續時間內利用系統同定的方法，從實驗得到的輸入輸出數據，可以定義內含物理現象的電池模型的參數。

電池狀態的推定

眾所周知，電池的狀態主要包括SOC（剩余電量）和SOH（剩余壽命）。

SOC的推定方法可以總結為：

1. 依據放電試驗的SOC測定
2. 基于電池端子電壓的SOC推定
3. 基于OCV測試的SOC推定
4. 依據安時積分法的SOC推定
5. 基于模型的SOC推定

SOH的推定方法可以總結為：

1. 利用完全充放電測試容量的SOH推定
2. 依據內部電阻查表的SOH推定
3. Bookkeeping的SOH推定
4. 基于模型的SOH推定

內部電阻的推定方法可以總結為：

1. 依據I-V特性（電流電壓特性）的線性回歸的內部電阻推定
2. 從step跳躍應答的內部電阻推定
3. 依據電阻測試的內部電阻推定
4. 基于模型的內部電阻推定

雖然有各種推算方法，但是我們需要根據：要求精度，使用環境，硬件條件，在線連續推定的必要性等等來具體判斷用什么方法。

下面我們就給出比較復雜的利用卡爾曼濾波以及安時積分推算SOC的Matlab建模方法。

[ex_batt.m]3

%%利用EKF的狀態推定

％離散時間的狀態空間建模

f_ct = @(x, u) f(x, u, th0); %模型參數的固定

f_dt = c2d_euler(f_ct, Ts); ％euler法的離散化

%fd的雅可比數列

Ad = expm(A(th0)*Ts); ％f的雅可比解析計算

Q = diag([0.1Ts)^2, 1e-6ones(1, Nd)]); %系統外亂

R = 0.01^2; %觀測外亂

xhat = zeros(Nd+1, N); %狀態推算值

P = zeros(Nd+1, Nd+1, N); %推定分散

%初期推算值

SOChat0 = OCV2SOC(ym(1));

xhat(:,1) = [SOChat0; zeros(Nd,1)];

P(:,:,1) = diag([1e2, 1e-4*ones(1, Nd)]);

%更新時間

for k=2:N

[xhat(:,k),P(:,:,k)] =...

ekf(@(x) f_dt(x, um(k-1)), @x h(x, um(k), th0),...

@(x) Ad, C, Q, R, ym(k), xhat(:,k-1), P(:,:,k-1));

end

%表示結果

figure, plot_SOC(t, xhat(1,:), squeeze(P(1,1,:))',SOC);

EKF_RMSE = sqrt(mean((xhat(1,:) - SOC).^2,2));

%%安時積分法

SOC_cc = zeros(1, N);

SOC_cc(1) =SOChat(0);

for k=1:N-1

SOC_cc(k+1) = SOC_cc(k) + um(k)*Ts/FCC*100;

end

%表示結果

figure, plot_SOC(t, SOC_cc, [ ], SOC);

CC_RMSE = sqrt(mean((SOC_cc - SOC).^2,2));

上述euler的離散時間的狀態方程式變化方法為：

[c2d_euler.m]

function fd = c2d_euler(fc, h)

fuction x_new = fproto(x, u)

    x_new = x + h*fc(x, u);

end

fd = @fproto;

end

SOC推算結果的表示

[plot_SOC.m]

function [ ] = plot_SOC(t, SOChat, P, SOC)

subplot(2,1,1), hold on

plot(t/60, SOC, 'r', 'LineWidth', 2)

plot(t/60, SOChat, 'b', 'LineWidth', 0.5)

xlim([0 t(end)/60]), ylim([0 100]);

xlabel('Time [min]'), ylabel('SOC [%]')

legend('True', 'Estimated', 'Location', 'Best')

box on



subplot(2,1,2), hold on

if ~isempty(P)

    hold on

    patch([t, fliplr(t)]/60,...

        [(SOChat-SOC-2*sqrt(P)), fliplr(SOChat-SOC+2*sqrt(P))],...

        ' ','FaceColor', [0.5, 0.5, 1], 'FaceAlpha', 0.5,...

        'EdgeAlpha', 0)

end

plot(t/60, SOChat-SOC, 'b')

plot([0 t(end)/60], [0 0], 'r', 'LineWidth', 2)

xlabel('Time [min]'), ylabel('Error [SOC%]')

xlim([0 t(end)/60]), ylim([-5 3]);

box on

end

上述為線性卡爾曼濾波的推算方法，非線性卡爾曼濾波的matlab建模方式基本相同，這里就不做介紹了。另外，本文是以EV為對象所建模，以HEV為對象時，由于電流波形的差異，推定的結果會有變故。HEV的剎車回生充電，內燃機回生充電，SOC都是在50%附近控制，所以對于輸入的電流波形，不會有低頻率的成分。

對于實際使用電池的系統來說，有著電池組內的誤差，生產誤差，環境溫度誤差等各項影響因素，所以未來我們使用模型推算電池狀態時，不僅要考慮到這些誤差因素的影響，更要考慮到鋰電池材料進化的影響。未必復雜的算法適用于每個產品，我們需要根據實際的條件以及性價比等各種綜合因素來挑選推算的方法，希望讀者朋友們互相切磋，共同進步。