CN204666798U - 一种基于fpga的锂离子电池soc估计设备 - Google Patents

一种基于fpga的锂离子电池soc估计设备 Download PDF

Info

Publication number
CN204666798U
CN204666798U CN201520385559.0U CN201520385559U CN204666798U CN 204666798 U CN204666798 U CN 204666798U CN 201520385559 U CN201520385559 U CN 201520385559U CN 204666798 U CN204666798 U CN 204666798U
Authority
CN
China
Prior art keywords
centerdot
soc
fpga
lithium ion
ion battery
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201520385559.0U
Other languages
English (en)
Inventor
马彦
何鹏才
王留
周秀文
孙延帅
茹敬佩
李炳思
殷明月
陈虹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN201520385559.0U priority Critical patent/CN204666798U/zh
Application granted granted Critical
Publication of CN204666798U publication Critical patent/CN204666798U/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Tests Of Electric Status Of Batteries (AREA)

Abstract

本实用新型提供了一种基于FPGA的锂离子电池SOC估计设备,包括锂离子电池、充放电设备、接线端、FPGA板卡、上位机、报警电路、一分二串口线和信号线,充放电设备通过接线端与锂离子电池相连接,充放电设备通过一分二串口线的左端与FPGA板卡相连,上位机通过一分二串口线的右端与FPGA板卡相连,FPGA板卡通过信号线与报警电路相连,FPGA板卡为Altera公司的型号为DE2-70的FPGA板卡,其中包括UART通信电路。本实用新型采用FPGA估计电池的SOC,解决了现有电池管理系统所使用的处理器在估计SOC时存在内存小、运行速度慢的问题;提出的快速矩阵运算法降低了多维矩阵运算的复杂度,减少了系统在进行矩阵运算时的存储次数和计算次数,具有运行速度快、节约资源的优点。

Description

一种基于FPGA的锂离子电池SOC估计设备
技术领域
本发明属于电动汽车动力电池技术领域,具体涉及一种基于FPGA的锂离子电池SOC估计方法及估计设备。
背景技术
对于电动汽车电池管理系统(Battery Management System,简称BMS)而言,实时、高效、准确地估计出电池的荷电状态(State of Charge,简称SOC)至关重要,可以从SOC估计算法和处理器结构两方面进行改进。
针对锂离子电池的非线性模型,运用扩展卡尔曼滤波(Extended Kalman Filter,简称EKF)算法可以准确地估计出电池的SOC值。该方法不仅可以克服安时法无法确定SOC初值、累计误差大的问题,还能解决开路电压法不能在线估计SOC的问题。然而,扩展卡尔曼滤波算法(EKF)存在大量的多维矩阵加、减、乘法运算,运算过程复杂,占用内存较大,并且现有的单片机等处理器已无法满足系统对运行速度和内存的要求。
现有的电池管理系统中的处理器大多以单片机、ARM或DSP为核心。单片机采用分布式结构,其内存小、处理速度慢、计算能力差,导致电池管理系统的效率低,可靠性差。ARM和DSP较单片机在内存和计算速度方面有所改善,但是其本质还是基于指令的串行执行方式,在处理速度的提高上仍有限制。在实际应用中,为了减小处理器的压力,也常常采用多片单片机、ARM或DSP相结合,组成主控制器和从控制器的方式来实现电池管理系统不同部分的工作,然而这些方法造成了调试困难和系统的分割,降低了系统的稳定性。
现场可编程门阵列(Field Programmable Gate Array,简称FPGA)作为一种内部硬件结构可编程的信号处理器件,具有灵活且可配置、执行速度快、工作性能稳定可靠、存储容量大等优点。此外,Altera和Xilinx公司还提出将软核处理器嵌入到FPGA的可编程逻辑器件中,不仅能在软核中方便地实现算法,而且通过总线实现了微处理器和逻辑器件的通信。基于这些优点,FPGA适合实现大型和复杂的算法,能满足SOC估计算法对硬件的需求。
发明内容
本发明克服了现有电池管理系统中所用的处理器在完成含有多维矩阵运算的复杂算法时存在内存小、运行速度慢,内存和时钟资源消耗大的问题,提供了一种基于FPGA的锂离子 电池SOC估计方法和估计设备。
为解决上述技术问题,本发明是采用如下技术方案实现的:
一种基于FPGA的锂离子电池SOC估计方法,其特征在于步骤如下:
第一步、建立锂离子电池的二阶等效电路模型:
锂离子电池模型采用二阶RC等效电路模型,其离散化的状态方程和输出方程为:
SOC ( k + 1 ) V s ( k + 1 ) V d ( k + 1 ) = 1 0 0 0 1 - T τ s 0 0 0 1 - T τ d · SOC ( k ) V s ( k ) V d ( k ) + - T Q N T C s T C d · I ( k ) - - - ( 1 )
VD(k)=Voc(SOC(k))-Vs(k)-Vd(k)-Ri·I(k)  (2)
其中,选择电池内部的荷电状态SOC、极化电压Vs和Vd作为状态变量,则在k时刻,x(k)=[SOC(k) Vs(k) Vd(k)]T;τs=Rs·Cs,τd=Rd·Cd,Rs、Rd是电池的极化内阻,Cs、Cd是电池的极化电容;Ri是电池的欧姆内阻;T是采样间隔时间;QN为电池额定容量;Voc(SOC(k))为k时刻电池的开路电压,VD(k)为k时刻电池的端电压,I(k)为k时刻电池的充放电电流;
引入过程噪声w(k)和测量噪声v(k),式(1)、(2)可表示为
x(k+1)=A·x(k)+B·I(k)+w(k)  (3)
y(k)=VD(k)+v(k)=C·x(k)+v(k)  (4)
其中,w(k)和v(k)是互不相关的零均值高斯白噪声y(k)表示实际测量的电池端电压;
系数矩阵: A = 1 0 0 0 1 - T τ s 0 0 0 1 - T τ d , B = - T Q N T C s T C d T , C = ∂ y ( k ) ∂ x ( k ) T ;
第二步、结合二阶等效电路模型,建立基于EKF的SOC估计算法
结合式(3)、式(4)和第一步中的锂离子电池二阶等效电路模型,建立基于EKF的SOC估计算法估计电池的SOC,具体为:
(1)初始化
k=0时确定过程噪声的方差Q=E[w(k),w(k)T]、测量噪声的方差R=E[v(k),v(k)T]、初始状态估计值初始状态的协方差 P ( 0 ) = E [ ( x ‾ - ( 0 ) - x ^ ( 0 ) ) ( x ^ - ( 0 ) - x ^ ( 0 ) ) T ] ;
(2)滚动更新
在k=0,1,2…时刻,
状态预测: x ^ - ( k + 1 ) = A · x ^ ( k ) + B · I ( k ) - - - ( 5 )
协方差预测:P-(k+1)=A·P(k)·AT+Q  (6)
卡尔曼增益:L(k+1)=P-(k+1)·C(k)T·[Ck+1·P-(k+1)·C(k)T+R]-1  (7)
状态更新: x ^ ( k + 1 ) = x ^ - ( k + 1 ) + L ( k + 1 ) · [ y ( k + 1 ) - V d ( k + 1 ) ] - - - ( 8 )
协方差更新:P(k+1)=[E-L(k+1)·Ck]·P-(k+1)  (9)
本步骤过程(2)中的公式(5)至公式(9)循环迭代,所要求解的SOC值为的状态之一,分离的状态量即可得到滚动更新的锂电池SOC值;
第三步、分析快速矩阵运算法的原理,并把第二步中的EKF算法进行快速矩阵运算分解,其过程为:
(1)快速矩阵运算法原理
设矩阵A、B、C和D分别满足
A m × s · B s × n = a 11 · b 11 + . . . + a 1 q · b q 1 + . . . + a 1 s · b s 1 . . . a 11 · b 1 e + . . . + a 1 q · b qe . . . + a 1 s · b se . . . a 11 · b 1 n + . . . + a 1 q · b qn . . . + a 1 s · b sn . . . . . . . . . . . . . . . a p 1 · b 11 + . . . + a pq · b q 1 + . . . + a ps · b s 1 . . . a p 1 · b 1 e + . . . + a pq · b qe . . . + a ps · b se . . . a p 1 · b 1 n + . . . + a pq · b qn . . . + a ps · b sn . . . . . . . . . . . . . . . a m 1 · b 11 + . . . + a mq · b q 1 + . . . + a ms · b s 1 . . . a m 1 · b 1 e + . . . + a mq · b qe . . . + a ms · b se . . . a m 1 · b 1 n + . . . + a mq · b qn . . . + a ms · b sn - - - ( 10 )
C m × d · D d × n = c 11 · d 11 + . . . + c 1 f · d f 1 + . . . + c 1 d · d d 1 . . . c 11 · d 1 e + . . . + c 1 f · d fe . . . + c 1 d · d de . . . c 11 · d 1 n + . . . + c 1 f · d fn . . . + c 1 d · d dn . . . . . . . . . . . . . . . c p 1 · d 11 + . . . + c pf · d f 1 + . . . + c pd · d d 1 . . . c p 1 · d 1 e + . . . + c pf · d fe . . . + c pd · d de . . . c p 1 · d 1 n + . . . + c pf · d fn . . . + c pd · d dn . . . . . . . . . . . . . . . c m 1 · d 11 + . . . + c mf · d f 1 + . . . + c md · d d 1 . . . c m 1 · d 1 e + . . . + c mf · d fe . . . + c md · d de . . . c m 1 · d 1 n + . . . + c mf · d fn . . . + c md · d dn - - - ( 11 )
所以,任意维数的矩阵多项式运算可表示为:
把式(10)、(11)带入式(12)中,得到快速矩阵运算法的分解式为:
g 11 = ( a 11 · b 11 + . . . + a 1 q · b q 1 + . . . + a 1 s · b s 1 ) ± ( c 11 · d 11 + . . . + c 1 f · d f 1 + . . . + c 1 d · d d 1 ) ± . . . g 1 e = ( a 11 · b 1 e + . . . + a 1 q · b qe + . . . + a 1 s · b se ) ± ( c 11 · d 1 e + . . . + c 1 f · d fe + . . . + c 1 d · d de ± ) . . . g pe = ( a p 1 · b 1 e + . . . + a pq · b qe + . . . + a ps · b se ) ± ( c p 1 · d 1 e + . . . + c pf · d fe + . . . + c pd · d de ) ± . . . g mn = ( a m 1 · b 1 n + . . . + a mq · b qn + . . . + a ms · b sn ) ± ( c m 1 · d 1 n + . . . + c mf · d fn + . . . + c md · d dn ) ± . . . - - - ( 13 )
快速矩阵运算法只需要按式(13)计算矩阵Gm×n内各元素即可;
(2)EKF算法分解
按照快速矩阵运算法可将EKF算法中式(5)至式(9)按式(13)进行分解,其过程如下:
1)分解状态预测值
设状态预测值 x ^ - ( k + 1 ) = x ^ 1 - ( k + 1 ) x ^ 2 - ( k + 1 ) x ^ 3 - ( k + 1 ) T ;
状态量 x ^ ( k ) = x ^ 1 ( k ) x ^ 2 ( k ) x ^ 3 ( k ) T = x ( k ) ; 系数矩阵 A = a 11 0 0 0 a 22 0 0 0 a 33 , B = b 1 b 2 b 3 ; 则状态预测式(5)的分解式为:
x ^ 1 - ( k + 1 ) = a 11 · x ^ 1 ( k ) + b 1 · I ( k ) x ^ 2 - ( k + 1 ) = a 22 · x ^ 2 ( k ) + b 2 · I ( k ) x ^ 3 - ( k + 1 ) = a 33 · x ^ 3 ( k ) + b 3 · I ( k ) - - - ( 14 )
2)分解协方差矩阵
设协方差预测值 P - ( k + 1 ) = p 11 - ( k + 1 ) p 12 - ( k + 1 ) p 13 - ( k + 1 ) p 21 - ( k + 1 ) p 22 - ( k + 1 ) p 23 - ( k + 1 ) p 31 - ( k + 1 ) p 32 - ( k + 1 ) p 33 - ( k + 1 ) ;
协方差 P ( k ) = p 11 ( k ) p 12 ( k ) p 13 ( k ) p 21 ( k ) p 22 ( k ) p 23 ( k ) p 31 ( k ) p 32 ( k ) p 33 ( k ) ; 过程噪声的方差 Q = Q 11 0 0 0 Q 22 0 0 0 Q 33 ; 则协方差预测式(6)的分解式为:
p 11 - ( k + 1 ) = a 11 2 · p 11 ( k ) + Q 11 p 12 - ( k + 1 ) = a 11 · a 22 · p 12 ( k ) p 13 - ( k + 1 ) = a 11 · a 33 · p 13 ( k ) p 21 - ( k + 1 ) = a 11 · a 22 · p 21 ( k ) p 22 - ( k + 1 ) = a 22 2 · p 22 ( k ) + Q 22 p 23 - ( k + 1 ) = a 22 · a 33 · p 23 ( k ) p 31 - ( k + 1 ) = a 33 · a 11 · p 31 ( k ) p 32 - ( k + 1 ) = a 22 · a 33 · p 32 ( k ) p 33 - ( k + 1 ) = a 33 2 · p 33 ( k ) + Q 33 - - - ( 15 )
3)分解卡尔曼增益
设卡尔曼增益L(k+1)=[L11(k+1) L21(k+1) L31(k+1)]T;系数矩阵C=[c11 c12 c13];则卡尔曼增益式(7)的分解式为:
L 11 ( k + 1 ) = ( p 11 - ( k + 1 ) · c 11 + p 12 - ( k + 1 ) · c 12 + p 13 - ( k + 1 ) · c 13 ) · ( g ) - 1 L 21 ( k + 1 ) = ( p 21 - ( k + 1 ) · c 11 + p 22 - ( k + 1 ) · c 12 + p 23 - ( k + 1 ) · c 13 ) · ( g ) - 1 L 31 ( k + 1 ) = ( p 31 - ( k + 1 ) · c 11 + p 32 - ( k + 1 ) · c 12 + p 33 - ( k + 1 ) · c 13 ) · ( g ) - 1 - - - ( 16 )
其中, g = c 11 2 · p 11 - + p 12 - · c 11 · c 12 + c 11 · c 13 · p 13 - + p 21 - · c 11 · c 12 + c 12 2 · p 22 - + c 12 · c 13 · p 23 - + c 11 · c 13 · p 31 - + c 12 · c 13 · p 32 - + c 13 2 · p 33 - + R ;
4)分解状态更新值
状态更新值式(7)的分解式为: 
x ^ 11 ( k + 1 ) = x ^ 11 - ( k + 1 ) + L 11 ( k + 1 ) · h x ^ 21 ( k + 1 ) = x ^ 21 - ( k + 1 ) + L 21 ( k + 1 ) · h x ^ 31 ( k + 1 ) = x ^ 31 - ( k + 1 ) + L 31 ( k + 1 ) · h - - - ( 17 )
其中,h=yk+1-VD(k+1);
5)分解协方差矩阵更新值
协方差矩阵更新值式(9)的分解式为: 
P 11 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 11 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 21 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 31 - ( k + 1 ) P 12 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 12 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 22 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 32 - ( k + 1 ) P 13 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 13 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 23 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 33 - ( k + 1 ) P 21 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 11 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 21 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 31 - ( k + 1 ) P 22 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 11 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 22 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 32 - ( k + 1 ) P 23 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 13 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 23 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 33 - ( k + 1 ) P 31 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 11 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 21 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 31 - ( k + 1 ) P 32 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 12 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 22 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 32 ( k + 1 ) P 33 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 13 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 23 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 33 - ( k + 1 ) - - - ( 18 )
最终,在FPGA中实现EKF算法时,只需要计算出分解式(14)至式(18)即可;
第四步、给锂离子电池充放电,同时采集其电压和电流数据:
第一步至第三步完成以后,在充放电设备中设置锂离子电池的充放电电流工况,开始为电池进行充放电,并将该工况下采集到的电压y(k)和电流I(k)的数据通过UART传输给FPGA板卡;
第五步、建立UART通信,读取第四步中的电压电流数据,并传输第六步输出的SOC值:
在FPGA中建立UART数据通信协议,一方面为第六步中的SOC估计器传送第四步中输出的电压y(k)和电流I(k)的数据,另一方面将第六步中所实时估计的SOC值传输给上位机显示;
第六步、按第四步的EKF分解算法在FPGA中建立SOC估计器,把第五步中的电压电流数据作为输入,估计出锂离子电池的SOC:
建立SOC估计器,嵌入到FPGA中,并将第三步中的EKF算法分解式(14)至式(18)在该估计器中实现,然后实时读取第四步中的电压y(k)和电流I(k)的数据,估计出电池的SOC值,并在SOC值超出设定的安全范围时由FPGA的I/O口向报警电路发送报警信号;
第七步、上位机监测与报警提示:
上位机接收第六步所传来的SOC值,并实时显示出SOC的变化曲线,同时第六步若发出报警信号,报警电路报警提示驾驶员。
一种基于FPGA的锂离子电池SOC估计设备,其特征在于,包括锂离子电池、充放电设备、接线端、FPGA板卡、上位机、报警电路、一分二串口线,信号线,充放电设备通过接线端与锂离子电池相连接,充放电设备通过一分二串口线的左端与FPGA板卡相连,上位机通过一分二串口线的右端与FPGA板卡相连,FPGA板卡通过信号线与报警电路相连,FPGA板卡为Altera公司的型号为DE2-70的FPGA板卡,其中包括UART通信电路。
与现有技术相比本发明的有益效果是:
本发明可实时监控锂离子电池SOC的变化情况,在SOC值超出安全范围时作出准确预警,这样驾驶员就可以清楚得知电池的使用情况,合理安排出行计划,并能良好的管理电池,延长电池的使用寿命;
快速矩阵运算法,可将复杂的多维矩阵的运算转换成矩阵内元素之间的运算,在FPGA内易于实现,减少了矩阵运算的存储次数和计算次数,在FPGA内节约了大量的存储空间和时钟资源,提高了系统完成锂离子电池SOC估计的效率;
在FPGA中采用EKF算法估计锂离子电池的SOC,解决了现有电池管理系统所使用的处理器采用EKF算法估计电池的SOC存在内存小、运行速度慢的问题。
附图说明
下面结合附图对本发明作进一步的说明:
图1是一种基于FPGA的锂离子电池SOC估计方法的流程图;
图2是一种基于FPGA的锂离子电池SOC估计设备的示意图;
图3是一种基于FPGA的锂离子电池SOC估计方法中的EKF算法流程图;
图4是一种基于FPGA的锂离子电池SOC估计设备的UART通信电路图;
图5是一种基于FPGA的锂离子电池SOC估计方法的基于Qsys搭建的Nios II系统;
图6是一种基于FPGA的锂离子电池SOC估计设备的报警电路的电路图;
图7是3.7V 300mAh锰酸锂电池SOC估计运行结果图。
图中:1.锂离子电池,2.充放电设备,3.接线端,4.FPGA板卡,5.上位机,6.报警电路,7.一分二串口线,8.信号线。
具体实施方式
下面结合附图对本发明作详细的描述:
所述的一种基于FPGA的锂离子电池SOC估计设备,如图2所示,包括锂离子电池1、充放电设备2、接线端3、FPGA板卡4、上位机5、报警电路6、一分二串口线7,信号线8,充放设备2通过接线端3与锂离子电池1相连接,充放电设备2通过一分二串口线7的左端与FPGA板卡4相连,上位机5通过一分二串口线7的右端与FPGA板卡4相连,FPGA板卡4通过信号线8与报警电路6相连。FPGA板卡4为Altera公司的型号为DE2-70的FPGA板卡,其中包括UART通信电路。
所述的一种基于FPGA的锂离子电池SOC估计方法,如图1所示,包括七个步骤:建立 锂离子电池的二阶等效电路模型;建立基于EKF的SOC估计算法;快速矩阵运算法的原理分析与EKF算法的快速矩阵运算分解;电压电流数据采集;UART通信;在FPGA中建立SOC估计器;上位机监测与报警提示。
所述的一种基于FPGA的锂离子电池SOC估计方法的具体步骤为
第一步、建立锂离子电池二阶等效电路模型:
锂离子电池模型采用二阶RC等效电路模型,其离散化的状态方程和输出方程为:
SOC ( k + 1 ) V s ( k + 1 ) V d ( k + 1 ) = 1 0 0 0 1 - T τ s 0 0 0 1 - T τ d · SOC ( k ) V s ( k ) V d ( k ) + - T Q N T C s T C d · I ( k ) - - - ( 1 )
VD(k)=Voc(SOC(k))-Vs(k)-Vd(k)-Ri·I(k)  (2)
其中,选择电池内部的荷电状态SOC、极化电压Vs和Vd作为状态变量,则在k时刻,x(k)=[SOC(k) Vs(k) Vd(k)]T;τs=Rs·Cs,τd=Rd·Cd,Rs、Rd是电池的极化内阻,Cs、Cd是电池的极化电容;Ri是电池的欧姆内阻;T是采样间隔时间;QN为电池额定容量;Voc(SOC(k))为k时刻电池的开路电压,VD(k)为k时刻电池的端电压,I(k)为k时刻电池的充放电电流;
引入过程噪声w(k)和测量噪声v(k),式(1)、(2)可表示为:
x(k+1)=A·x(k)+B·I(k)+w(k)  (3)
y(k)=VD(k)+v(k)=C·x(k)+v(k)  (4)
其中,w(k)和v(k)是互不相关的零均值高斯白噪声;y(k)表示实际测量的电池端电压; 系数矩阵 A = 1 0 0 0 1 - T τ s 0 0 0 1 - T τ d , B = - T Q N T C s T C d T , C = ∂ y ( k ) ∂ x ( k ) T ;
第二步、结合二阶等效电路模型,建立基于EKF的SOC估计算法
本发明选用EKF算法完成电池的SOC估计,EKF算法可滤除噪声与干扰,得到动态的系统状态最优解,适用于锂离子电池的非线性模型,可准确地估计出电池的SOC;
结合式(3)、式(4)和第一步中的锂离子电池二阶等效电路模型,建立基于EKF的SOC估计算法估计电池的SOC,如图3所示,具体过程为:
(1)初始化
k=0时确定过程噪声的方差Q=E[w(k),w(k)T]、测量噪声的方差R=E[v(k),v(k)T]、初始状态估计值初始状态的协方差 P ( 0 ) = E [ ( x ‾ - ( 0 ) - x ^ ( 0 ) ) ( x ^ - ( 0 ) - x ^ ( 0 ) ) T ] ;
(2)滚动更新
在k=0,1,2…时刻,
状态预测: x ^ - ( k + 1 ) = A · x ^ ( k ) + B · I ( k ) - - - ( 5 )
协方差预测:P-(k+1)=A·P(k)·AT+Q  (6)
卡尔曼增益:L(k+1)=P-(k+1)·CT·[C·P-(k+1)·CT+R]-1  (7) 
状态更新: x ^ ( k + 1 ) = x ^ - ( k + 1 ) + L ( k + 1 ) · [ y ( k + 1 ) - V d ( k + 1 ) ] - - - ( 8 )
协方差更新:P(k+1)=[E-L(k+1)·Ck]·P-(k+1)  (9)
(3)本步骤过程(2)中的公式(5)至公式(9)循环迭代,所要求解的SOC值为的状态之一,分离的状态量即可得到滚动更新的锂电池SOC值;
由图3可知在计算电池SOC的每次循环初始时,输入变量不仅需要电压y(k)和电流I(k)的值,同时也需要初值SOC(k),极化电压初值Vs(k)和Vd(k),协方差初值P(k)(k=0,1,…1823),每完成一次输入,就能估计出电池的一个SOC值,并更新状态值和协方 差值,将其作为下一次循环的初始值,进入到下一次循环计算SOC中去。就这样在循环体内不断运行EKF算法,直到完成实验所设定的1823组数据的采样,得到相应的电池SOC估计曲线。
第三步、快速矩阵运算法的原理分析与EKF算法的快速矩阵运算分解
(1)快速矩阵算法原理
本发明提出了快速矩阵运算法,把复杂的多维矩阵的加、减、乘法运算转换成矩阵内元素间的运算,节省了大量的系统资源,加快了系统运行速度,快速矩阵运算法原理如下:
设矩阵A、B、C和D分别满足:
A m × s · B s × n = a 11 · b 11 + . . . + a 1 q · b q 1 + . . . + a 1 s · b s 1 . . . a 11 · b 1 e + . . . + a 1 q · b qe . . . + a 1 s · b se . . . a 11 · b 1 n + . . . + a 1 q · b qn . . . + a 1 s · b sn . . . . . . . . . . . . . . . a p 1 · b 11 + . . . + a pq · b q 1 + . . . + a ps · b s 1 . . . a p 1 · b 1 e + . . . + a pq · b qe . . . + a ps · b se . . . a p 1 · b 1 n + . . . + a pq · b qn . . . + a ps · b sn . . . . . . . . . . . . . . . a m 1 · b 11 + . . . + a mq · b q 1 + . . . + a ms · b s 1 . . . a m 1 · b 1 e + . . . + a mq · b qe . . . + a ms · b se . . . a m 1 · b 1 n + . . . + a mq · b qn . . . + a ms · b sn - - - ( 10 )
C m × d · D d × n = c 11 · d 11 + . . . + c 1 f · d f 1 + . . . + c 1 d · d d 1 . . . c 11 · d 1 e + . . . + c 1 f · d fe . . . + c 1 d · d de . . . c 11 · d 1 n + . . . + c 1 f · d fn . . . + c 1 d · d dn . . . . . . . . . . . . . . . c p 1 · d 11 + . . . + c pf · d f 1 + . . . + c pd · d d 1 . . . c p 1 · d 1 e + . . . + c pf · d fe . . . + c pd · d de . . . c p 1 · d 1 n + . . . + c pf · d fn . . . + c pd · d dn . . . . . . . . . . . . . . . c m 1 · d 11 + . . . + c mf · d f 1 + . . . + c md · d d 1 . . . c m 1 · d 1 e + . . . + c mf · d fe . . . + c md · d de . . . c m 1 · d 1 n + . . . + c mf · d fn . . . + c md · d dn - - - ( 11 )
所以,任意维数的矩阵多项式运算可表示为:
把式(10)、式(11)带入式(12)中,得到快速矩阵运算法的分解式:
g 11 = ( a 11 · b 11 + . . . + a 1 q · b q 1 + . . . + a 1 s · b s 1 ) ± ( c 11 · d 11 + . . . + c 1 f · d f 1 + . . . + c 1 d · d d 1 ) ± . . . g 1 e = ( a 11 · b 1 e + . . . + a 1 q · b qe + . . . + a 1 s · b se ) ± ( c 11 · d 1 e + . . . + c 1 f · d fe + . . . + c 1 d · d de ± ) . . . g pe = ( a p 1 · b 1 e + . . . + a pq · b qe + . . . + a ps · b se ) ± ( c p 1 · d 1 e + . . . + c pf · d fe + . . . + c pd · d de ) ± . . . g mn = ( a m 1 · b 1 n + . . . + a mq · b qn + . . . + a ms · b sn ) ± ( c m 1 · d 1 n + . . . + c mf · d fn + . . . + c md · d dn ) ± . . . - - - ( 13 )
因此按式(13)计算矩阵Gm×n内各元素即可,而常规的矩阵运算方法与快速矩阵运算法有所不同。设a=[a0 … ai … an],b=[b0 … bi … bn]T(0≤i≤n),ci=ai·bi,则采用常规方法完成行向量a与列向量b相乘的思路是:ci=ci+ai·bi(c0=0),即每完成一次ai·bi的运算,需要存储一次中间结果,然后结果与上一个ci值相加,并赋值给ci存储。常规方法存在两个明显的缺点:一是增加了矩阵内元素相乘的中间结果的存储次数;二是若矩阵内含有0元素,对于含0项的乘法、加法和减法,公式ci=ci+ai·bi依然会执行,增加了计算次数。与常规方法相比,快速矩阵运算法能有效的减少矩阵运算的计算次数和数据的存储次数,缩短了系统运行时间,节约了系统计算内存,从而提高了系统运行效率,对于实现类似EKF这种含有大量矩阵运算的算法非常有益。
表1分别给出了采用常规方法和快速矩阵运算法完成公式(10)所用的存储次数和计算次数。从表1可以看出,常规方法的存储次数是快速矩阵运算法的s倍,而快速矩阵运算法的计算次数比常规方法少(视矩阵内0元素的多少而定,0元素越多,节省的计算次数越多)。
表1常规方法与快速矩阵运算法的存储次数和计算次数
(2)EKF算法分解
由于在FPGA中实现多维矩阵的运算比较复杂,所以本发明首先采用快速矩阵运算法对矩阵加、减、乘法运算进行分解,转换成矩阵元素之间的运算。
按照快速矩阵运算法可将EKF算法中式(5)至式(9)按式(13)进行分解,其过程如下:
6)分解状态预测值
设状态预测值 x ^ - ( k + 1 ) = x ^ 1 - ( k + 1 ) x ^ 2 - ( k + 1 ) x ^ 3 - ( k + 1 ) T ; 状态量  x ^ ( k ) = x ^ 1 ( k ) x ^ 2 ( k ) x ^ 3 ( k ) T = x ( k ) ; 系数矩阵 A = a 11 0 0 0 a 22 0 0 0 a 33 , B = b 1 b 2 b 3 ; 则状态预测式(5)的分解式为:
x ^ 1 - ( k + 1 ) = a 11 · x ^ 1 ( k ) + b 1 · I ( k ) x ^ 2 - ( k + 1 ) = a 22 · x ^ 2 ( k ) + b 2 · I ( k ) x ^ 3 - ( k + 1 ) = a 33 · x ^ 3 ( k ) + b 3 · I ( k ) - - - ( 14 )
7)分解协方差矩阵
设协方差预测值 P - ( k + 1 ) = p 11 - ( k + 1 ) p 12 - ( k + 1 ) p 13 - ( k + 1 ) p 21 - ( k + 1 ) p 22 - ( k + 1 ) p 23 - ( k + 1 ) p 31 - ( k + 1 ) p 32 - ( k + 1 ) p 33 - ( k + 1 ) ;
协方差 P ( k ) = p 11 ( k ) p 12 ( k ) p 13 ( k ) p 21 ( k ) p 22 ( k ) p 23 ( k ) p 31 ( k ) p 32 ( k ) p 33 ( k ) ; 过程噪声的方差 Q = Q 11 0 0 0 Q 22 0 0 0 Q 33 ; 则协方差预测式(6)的分解式为
p 11 - ( k + 1 ) = a 11 2 · p 11 ( k ) + Q 11 p 12 - ( k + 1 ) = a 11 · a 22 · p 12 ( k ) p 13 - ( k + 1 ) = a 11 · a 33 · p 13 ( k ) p 21 - ( k + 1 ) = a 11 · a 22 · p 21 ( k ) p 22 - ( k + 1 ) = a 22 2 · p 22 ( k ) + Q 22 p 23 - ( k + 1 ) = a 22 · a 33 · p 23 ( k ) p 31 - ( k + 1 ) = a 33 · a 11 · p 31 ( k ) p 32 - ( k + 1 ) = a 22 · a 33 · p 32 ( k ) p 33 - ( k + 1 ) = a 33 2 · p 33 ( k ) + Q 33 - - - ( 15 )
8)分解卡尔曼增益
设卡尔曼增益L(k+1)=[L11(k+1) L21(k+1) L31(k+1)]T;系数矩阵C=[c11 c12 c13];则卡尔曼增益式(7)的分解式为:
L 11 ( k + 1 ) = ( p 11 - ( k + 1 ) · c 11 + p 12 - ( k + 1 ) · c 12 + p 13 - ( k + 1 ) · c 13 ) · ( g ) - 1 L 21 ( k + 1 ) = ( p 21 - ( k + 1 ) · c 11 + p 22 - ( k + 1 ) · c 12 + p 23 - ( k + 1 ) · c 13 ) · ( g ) - 1 L 31 ( k + 1 ) = ( p 31 - ( k + 1 ) · c 11 + p 32 - ( k + 1 ) · c 12 + p 33 - ( k + 1 ) · c 13 ) · ( g ) - 1 - - - ( 16 )
其中, g = c 11 2 · p 11 - + p 12 - · c 11 · c 12 + c 11 · c 13 · p 13 - + p 21 - · c 11 · c 12 + c 12 2 · p 22 - + c 12 · c 13 · p 23 - + c 11 · c 13 · p 31 - + c 12 · c 13 · p 32 - + c 13 2 · p 33 - + R ;
9)分解状态更新值
状态更新值式(7)的分解式为: 
x ^ 11 ( k + 1 ) = x ^ 11 - ( k + 1 ) + L 11 ( k + 1 ) · h x ^ 21 ( k + 1 ) = x ^ 21 - ( k + 1 ) + L 21 ( k + 1 ) · h x ^ 31 ( k + 1 ) = x ^ 31 - ( k + 1 ) + L 31 ( k + 1 ) · h - - - ( 17 )
其中,h=yk+1-VD(k+1);
10)分解协方差矩阵更新值
协方差矩阵更新值式(9)的分解式为: 
P 11 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 11 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 21 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 31 - ( k + 1 ) P 12 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 12 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 22 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 32 - ( k + 1 ) P 13 ( k + 1 ) = ( 1 - L 11 ( k + 1 ) · c 11 ) · P 13 - ( k + 1 ) - L 11 ( k + 1 ) · c 12 P 23 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 33 - ( k + 1 ) P 21 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 11 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 21 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 31 - ( k + 1 ) P 22 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 11 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 22 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 32 - ( k + 1 ) P 23 ( k + 1 ) = - L 21 ( k + 1 ) · c 11 · P 13 - ( k + 1 ) + ( 1 - L 21 ( k + 1 ) · c 12 ) P 23 - ( k + 1 ) - L 11 ( k + 1 ) · c 13 P 33 - ( k + 1 ) P 31 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 11 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 21 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 31 - ( k + 1 ) P 32 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 12 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 22 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 32 ( k + 1 ) P 33 ( k + 1 ) = - L 13 ( k + 1 ) · c 11 P 13 - ( k + 1 ) - L 31 ( k + 1 ) · c 12 P 23 - ( k + 1 ) + ( 1 - L 31 ( k + 1 ) · c 13 ) · P 33 - ( k + 1 ) - - - ( 18 )
最终,在FPGA中实现EKF算法时,只需要计算出分解式(14)至式(18)即可。
特别说明的是式(7)中有一个求逆运算,通过分析发现C·P-(k+1)·C+R的运算结果是一个数,因此式(7)的求逆运算只需要按一个数的除法运算即可。
表2常规方法与快速矩阵运算法的存储次数和计算次数(EKF算法)
在FPGA板卡中采用式(14)至(18)的快速矩阵运算法和常规方法分别完成EKF算法中所需要的计算次数(包括加减法次数、乘法次数和求逆次数)和存储次数如表2所示。从表2可以看出,完成一次锂离子电池SOC估计,快速矩阵运算法比常规方法少存储130次,少计算50次,如果完成多次SOC估计,所要存储的次数和计算次数将会减少更多。
第四步、给锂离子电池充放电,同时采集其电压和电流数据
第一步至第三步完成以后,在充放电设备中设置锂离子电池的充放电电流工况,开始为电池进行充放电,并将该工况下采集到的电压y(k)和电流I(k)的数据通过UART传输给FPGA板卡4。
表3自定义电流工况
本实施例选4节3.7V 300mAh锰酸锂单体电池串联而成,电池的充放设备选用的是新威尔电池测试系统,自定义的电流工况为表3所示。设置完成后,开启充电设备,为电池充放电,并将电压和电流数据传输出去。
第五步、建立UART通信,读取第四步中的电压电流数据,并传输第六步输出的SOC值
本实施例选用的是一分二串口线,一端连接FPGA板卡4,另外两端分别连接充放电设备1和上位机5,一方面可为第六步中的SOC估计器传送第四步中输出的电压y(k)和电流I(k)的数据,另一方面可把第六步中所实时估计的SOC值传输给上位机显示,如图2所示。
本实施例使用ADM3202收发芯片和一个九针D-SUB连接器用于构成FPGA板卡4的UART串口通信,图4是其硬件原理图。其中UART_RXD是接收端,接收FPGA处理单元输出的SOC值,与FPGA芯片D21引脚相连接;UART_TXD是发送端,将接收端所接收的SOC值发送给上位机,与FPGA芯片G22引脚相连接;UART_RTS/UART_CTS是请求发送/允许发送信号,分别与FPGA芯片E21和D21引脚相连接。本实施例中传输速率设为19200bit/s。图5中接收和发送端与发光二极管相连接,以提示是否正在发送和接收信息。
本实施例中根据UART通信协议,用Verilog HDL设计了UART通信,不但提高了执行速度,而且将UART核心功能集成到了FPGA板卡4中,可以实现可靠、稳定、紧凑的数据传输。
第六步、按第四步的EKF分解算法在FPGA中建立SOC估计器,把第五步中的电压电流数据作为输入,估计出锂离子电池的SOC
FPGA为锂离子电池SOC估计算法提供了一个Qsys(SOPC)平台,在该平台上建立Nios II系统,嵌入到FPGA中,然后在该系统中建立EKF估计器估计出锂离子电池的SOC值。
本发明选用FPGA板卡4是Alera公司的DE2-70,其芯片型号为Cyclone II EP2C70F896C6N。在FPGA中配置Qsys硬件平台,并产生Nios II系统,嵌入到FPGA板卡4中,然后在该系统中建立EKF估计器估计出锂离子电池的SOC,通过FPGA板 卡4的I/O口和UART通信传输给上位机5和报警电路6。详细过程如下:
(1)在Quartus II软件中建立工程。
(2)建立所需要的Qsys系统。
Qsys工具是新一代SOPC builder开发工具,是Quartus II中的一个组件,具有很好的用户定制特性,可通过添加和删除标准组件或者配置不同特性的组件来构建具有不同性能、满足不同要求的处理器。配置出的Qsys系统如图5所示,包括Nios II CPU、代表身份的System ID、下载程序的JTAG UART、内部存储器OnChip RAM、输出SOC值的PIO、分析算法的运行时间和时钟消耗的Performance counter等内核,各内核通过Avalon总线与CPU连接起来。
(3)建立系统的顶层文件,通过顶层列化将Qsys系统与外部硬件电路连接起来。
(4)产生Nios II系统,并在Nios II IDE中完成EKF算法的编写与调试。
按照上述过程完成以后,将EKF算法分解式(14)至式(18)在该估计器中实现。最后,读取采集的电压y(k)和电流I(k)的数据,作为SOC估计器的输入,在FPGA板卡4中实时估计出锂离子电池的SOC值,并在SOC值超出设定的安全范围时由FPGA板卡4的I/O口向报警电路6发送报警信号。
第七步、上位机监测与报警提示
上位机接收第六步所传来的SOC值,并实时显示出SOC的变化曲线。 
同时第六步若发出报警信号,报警电路6报警提示驾驶员。报警电路6如图6所示,EN_SOC表示FPGA芯片发出的报警信号,与FPGA板卡4的I/O口相连接。当电池的SOC低于15%或达到100%时,EN_SOC变为高电平,蜂鸣器报警,同时红灯亮,提醒驾驶员及时充电或停止充电。
按照本发明完成实验,上位机5显示的结果如图7所示。其中曲线1是基于FPGA的锂离子电池SOC估计方法的锂离子电池SOC估计曲线,曲线2是基于安时法(所用的电流为设置工况的理论电流,不存在积分累计误差)的锂离子电池SOC估计参考曲线。从图7中可以看出,在0~1823s内,锂离子电池的SOC由89%下降到71%,估计曲线与参考曲线相比的误差在合理范围内,系统估计精度较高。
同时本发明通过Nios II系统中Performance counter内核仿真分析了使用快速矩阵运算法和常规方法实现EKF算法估计电池SOC的优劣,结果如表4。通过表4可以看出,利用快速矩阵运算法实现的锂离子电池SOC估计与常规方法相比,完成一次SOC估计所用的时间明 显减少了,如常规方法完成一次SOC估计所用的时间是0.500526ms,所用的时钟个数为25026,而快速矩阵运算法完成一次SOC估计所用的时间是0.086594ms,所用的时钟个数为4330,时间减少了近5倍。
表4 EKF算法性能分析 

Claims (1)

1.一种基于FPGA的锂离子电池SOC估计设备,其特征在于,包括锂离子电池(1)、充放电设备(2)、接线端(3)、FPGA板卡(4)、上位机(5)、报警电路(6)、一分二串口线(7)和信号线(8),充放电设备(2)通过接线端(3)与锂离子电池(1)相连接,充放电设备(2)通过一分二串口线(7)的左端与FPGA板卡(4)相连,上位机(5)通过一分二串口线(7)的右端与FPGA板卡(4)相连,FPGA板卡(4)通过信号线(8)与报警电路(6)相连,FPGA板卡(4)为Altera公司的型号为DE2-70的FPGA板卡,其中包括UART通信电路。
CN201520385559.0U 2015-06-05 2015-06-05 一种基于fpga的锂离子电池soc估计设备 Active CN204666798U (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201520385559.0U CN204666798U (zh) 2015-06-05 2015-06-05 一种基于fpga的锂离子电池soc估计设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201520385559.0U CN204666798U (zh) 2015-06-05 2015-06-05 一种基于fpga的锂离子电池soc估计设备

Publications (1)

Publication Number Publication Date
CN204666798U true CN204666798U (zh) 2015-09-23

Family

ID=54137168

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201520385559.0U Active CN204666798U (zh) 2015-06-05 2015-06-05 一种基于fpga的锂离子电池soc估计设备

Country Status (1)

Country Link
CN (1) CN204666798U (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104865535A (zh) * 2015-06-05 2015-08-26 吉林大学 一种基于fpga的锂离子电池soc估计方法和估计设备
CN109061520A (zh) * 2018-10-25 2018-12-21 杭州神驹科技有限公司 一种动力电池健康与功率状态在线估算方法及系统
CN109490611A (zh) * 2018-10-29 2019-03-19 宁波三星智能电气有限公司 一种嵌入式设备的时间统计方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104865535A (zh) * 2015-06-05 2015-08-26 吉林大学 一种基于fpga的锂离子电池soc估计方法和估计设备
CN109061520A (zh) * 2018-10-25 2018-12-21 杭州神驹科技有限公司 一种动力电池健康与功率状态在线估算方法及系统
CN109061520B (zh) * 2018-10-25 2021-04-23 杭州神驹科技有限公司 一种动力电池健康与功率状态在线估算方法及系统
CN109490611A (zh) * 2018-10-29 2019-03-19 宁波三星智能电气有限公司 一种嵌入式设备的时间统计方法
CN109490611B (zh) * 2018-10-29 2021-03-05 宁波三星智能电气有限公司 一种嵌入式设备的时间统计方法

Similar Documents

Publication Publication Date Title
CN106569136B (zh) 一种电池健康状态在线估计方法及系统
CN105676134B (zh) 一种车用锂离子动力电池的soh估算方法
US11526639B2 (en) Fractional-order KiBaM battery model considering nonlinear capacity characteristics and parameter identification method
CN107576919A (zh) 基于armax模型的动力电池荷电状态估算系统及方法
CN104101838B (zh) 动力电池系统及其荷电状态、最大充放电功率估算方法
CN104865535A (zh) 一种基于fpga的锂离子电池soc估计方法和估计设备
CN103926538A (zh) 基于aic准则的变阶数rc等效电路模型及实现方法
CN102937704A (zh) 一种动力电池rc等效模型的辨识方法
CN104977545A (zh) 一种动力电池的荷电状态估计方法及系统
CN204666798U (zh) 一种基于fpga的锂离子电池soc估计设备
CN103197251A (zh) 一种动力锂电池二阶rc等效模型的辨识方法
CN103744026A (zh) 基于自适应无迹卡尔曼滤波的蓄电池荷电状态估算方法
CN104392080A (zh) 一种锂电池分数阶变阶等效电路模型及其辨识方法
CN106338695A (zh) 一种基于粒子群算法的电池模型参数辨识方法
CN110045292A (zh) 基于大数据和bp神经网络的锂离子电池SOC预测方法
CN102944848A (zh) 一种动力电池剩余电量实时评估方法及其装置
CN106842056A (zh) 一种基于双步在线智能优化算法动力电池峰值功率估计方法
CN106772081A (zh) 基于扩展等效电路模型的电池极限充放电电流估计方法
CN112816875A (zh) 电动车辆电池云端管理系统、方法、介质及云端服务器
CN105203968A (zh) 一种铅酸蓄电池剩余电量的在线测量系统
CN104537166A (zh) 一种动力电池的等效电路模型方法
CN111123136A (zh) 一种锂离子电池健康状态在线测量装置及在线测量方法
CN110135527A (zh) 一种高效能的无人机锂离子电池荷电状态估计系统及方法
CN112946480B (zh) 一种提高soc估计实时性的锂电池电路模型简化方法
CN107255786B (zh) 一种磷酸铁锂电池loc模型

Legal Events

Date Code Title Description
C14 Grant of patent or utility model
GR01 Patent grant