一种gpu加速的电力潮流雅可比矩阵的qr分解方法
【专利摘要】本文公开了一种GPU加速的电力潮流雅可比矩阵的QR分解方法，所述方法包括：CPU中对雅可比矩阵J进行QR符号分解，得到Household变换矩阵V和上三角阵R阵的稀疏结构；根据R阵的稀疏结构，对矩阵J各列进行并行化分层；GPU中按层次递增的顺序计算分层QR分解内核函数SparseQR。本发明利用CPU控制程序的流程并处理基础数据和GPU处理密集的浮点运算相结合的模式提高了电力潮流雅可比矩阵QR分解的效率，解决了电力系统静态安全性分析中潮流计算耗时大的问题。
【专利说明】
一种GPU加速的电力潮流雅可比矩阵的QR分解方法
技术领域
[0001]本发明属于电力系统高性能计算应用领域，尤其涉及一种GPU加速的电力潮流雅可比矩阵的QR分解的GPU线程设计方法。
【背景技术】
[0002]潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，在系统规划设计和安排系统的运行方式时，采用离线潮流计算;在电力系统运行状态的实时监控中，则采用在线潮流计算。
[0003]而实际生产过程中，无论离线潮流和在线潮流计算都对潮流的计算速度有这比较高的要求。在涉及规划设计和安排运行方式的离线潮流中，因设备落地方案等情况复杂，需要仿真运行的种类多，潮流计算量大，单个潮流计算时间影响整体仿真时长;而在电力系统运行中进行的在线潮流计算对计算时间敏感度高，需要实时给出潮流计算结果，如在预想事故、设备退出运行对静态安全的影响的潮流计算中，系统需要计算大量预想事故下潮流分布，并实时地做出预想的运行方式调整方案。
[0004]传统的牛顿拉夫逊法潮流计算中，修正方程组求解占潮流计算时间的70%，修正方程组求解的计算速度影响程序的整体性能。而随着CPU计算速度提升的放缓，现阶段的单个潮流计算计算时间已经达到一个瓶颈。目前对潮流计算的加速方法主要集中在使用集群和多核服务器对多潮流进行粗粒度加速，实际成产中对单个潮流内部运算加速的研究较少。
[0005]GPU是一种众核并行处理器，在处理单元的数量上要远远超过CPU。传统上的GPU只负责图形渲染，而大部分的处理都交给了 CPU。现在的GPU已经法阵为一种多核，多线程，具有强大计算能力和极高存储器带宽，可编程的处理器。在通用计算模型下，GPU作为CPU的协处理器工作，通过任务合理分配分解完成高性能计算。
[0006]稀疏线性方程组求解计算具有并行性。对方程组系数矩阵进行QR符号分解后，得到Household变换矩阵V和上三角阵R阵的稀疏结构，根据R阵的稀疏结构，对矩阵J各列进行并行化分层。其中每层中的列的计算相互独立，没有依赖关系，天然可以被并行的计算处理，适合GPU加速。因此通过CPU和GPU之间合理的调度可以快速完成方程组系数矩阵进行QR分解，并求解稀疏线性方程组，国内外学者已经开始对GPU进行稀疏线性方程组加速求解的方法进行了研究，但是没有深入的优化线程设计，单纯从计算量的分配上研究计算线程设计，对线程计算方式，数据索引方式没有进行深入研究，无法使程序充分发挥GPU的优势。
[0007]因此，亟待解决上述问题。
【发明内容】
[0008]发明目的:针对现有技术的不足，本发明提供一种能大幅减少电力潮流雅可比矩阵QR分解计算时间并能提升潮流计算速度的一种GHJ加速的电力潮流雅可比矩阵的QR分解方法。
[0009]潮流计算:电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。
[0010]并行计算:相对于串行运算，是一种一次可执行多个指令的算法，目的是提高计算速度，及通过扩大问题求解规模，解决大型而复杂的计算问题。
[0011 ] GPU:图形处理器(英语:6瓜口11；^8?1'00688；[呢1]11；[1:，缩写:GPU)。
[0012]本发明公开了一种GPU加速的电力潮流雅可比矩阵的QR分解方法，所述方法包括:
[0013](I)CPU中对雅可比矩阵J进行QR符号分解，得到Household变换矩阵V和上三角阵R阵的稀疏结构;根据R阵的稀疏结构，对矩阵J各列进行并行化分层；
[0014](2)GPU中按层次递增的顺序计算分层QR分解内核函数SparseQR。
[0015]其中，所述步骤(I)中，并行化分层将矩阵J的η列归并到MaxLevel层中，属于同一层中的列并行进行QR分解;每层包含的列的数量为LeVelnum(k)，k表示层号;存储第k层中所有列号至映射表Mapk。
[0016]优选的，所述步骤(2)中，分层(^分解内核函数定义为3口3186(^〈他1。。1?，队111^如>，其线程块大小Nthread固定为128，当对k层进行计算时，线程块数量Nb1cks = Levelnum(k)，总线程数量为:Nb1cks XNthreads;按照层次递增的顺序，调用内核函数SparseQR〈Levelnum(k)，Nthreads〉分解属于第k层的所有列；SparseQR〈Levelnum(k)，Nthreads>的计算流程为:
[0017]((2.1)⑶DA自动为每个线程分配线程块索引blockID和线程块中的线程索引threadID；
[0018](2.2)将blockID和threadID赋值给变量bid和t，之后通过bid和t来索引bid号线程块中的t号线程；
[0019](2.3)第1^(1号线程块负责01?分解雅可比矩阵1的第」=1^?1{(1^(1)列；
[0020](2.4)第1^(1号线程块中，变量1从1递增到」-1，如果1?(1，」)#0，执行以下计算:
[0021]首先，计算中间变量0= 2￥(:[:11，；01'.J(1:n，j)，其中V(1:n，i)是Household变换矩阵V中第i列的第i?η行元素构成的向量，J(1:n，j)是雅克比矩阵J中第j列的第i?η行元素构成的向量;然后，采用公式氕1:11，」)=氕1:11，」)-|^(1:11，1)更新雅可比矩阵1的第」列；
[0022](2.5)计算第j列的Household变换向量；
[0023](2.6)更新了矩阵第」列:了(」，」)=&，了(」+1:11，」)=0。
[0024]优选的，所述步骤(2.4)中:
[0025]首先，计算中间变量0= 2￥(:[:11，；01'.J(1:n，j)，具体计算步骤为，
[0026]I)判断线程编号t是否小于n-1，否则该线程结束执行；
[0027]2)线程中的变量Temp+ = 2V( i+t，i) X J( i+t，j);
[0028]3)t = t+128，返回 I)执行；
[0029]4)线程中Temp的值赋给共享内存变量Cache[threadID]；
[0030]5)对线程块中的线程进行同步；
[0031]6)变量 M= 128/2;
[0032]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0033]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0034]9)对线程块中的线程进行同步；
[0035]10)M/ = 2；
[0036]11)返回7)执行；
[0037]12)P = Cache(0)；
[0038]然后，采用公式J(1: η，j) = J( 1: η，j)-βν(1: η，i)更新雅可比矩阵J的第j列，具体步骤为，
[0039]I)判断线程编号t是否小于n-1，否则线程结束执行；
[0040]2)J(t+i,j)=J(t+i,j)-PV(t+i,i);
[0041]3)t = t+128，返回 I)。
[0042]进一步，所述步骤(2.5)中，
[0043]首先，计算变量a2= J( j:n，j)T.J( j:n，j)，具体计算步骤如下:
[0044]I)判断线程编号t是否小于n-j，否则该线程结束执行；
[0045]2)线程中的变量Temp，Temp+ = J( j+t, j) X J( j+t, j);;
[0046]3)t = t+128，返回 I)执行；
[0047]4)线程中Temp的值赋给共享内存变量Cache[threadID]；；
[0048]5)对线程块中的线程进行同步；
[0049]6)变量 M= 128/2;
[0050]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0051 ]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0052 ]9)对线程块中的线程进行同步；
[0053]10)M/ = 2；
[0054]11)返回7)执行；
[0055]12)a2 = Cache(0)；
[0056]其次，计算，v(j:η，j) = J( j:η，j )-aej(j:n)，其中是ej是第j个元素为I的单位向量；
[0057]具体步骤如下:
[0058]I)判断线程编号t是否小于n-j，否则线程结束执行；
[0059]2)V(t+j, j)=J(t+j, j)-aej(t+j)；；
[0060]3)t = t+128，返回 I)。
[0061 ]然后，计算变量b2 = V( j:η，j )τ.V( j:η，j);
[0062]具体计算步骤如下:
[0063]I)判断线程编号t是否小于n-j，否则该线程结束执行；
[0064]2)线程中的变量为Temp，Temp+ = V( j+t，j) XV( j+t, j);
[0065]3)t = t+128，返回 I)执行；
[0066]4)线程中Temp的值赋给共享内存变量Cache[threadID]；
[0067]5)对线程块中的线程进行同步；
[0068]6)变量 M= 128/2;
[0069]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0070]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0071]9)对线程块中的线程进行同步；
[0072]10)M/ = 2；
[0073]11)返回7)执行；
[0074]12)b2 = Cache(0)；
[0075]最后，计算V (j:n, j)=V(j:n, j)/b；
[0076]具体步骤如下:
[0077]I)判断线程编号t是否小于n-j，否则线程结束执行；
[0078]2)V(t+j，j)=V(t+j，j)/b;
[0079]3)t = t+128，返回 I)。
[0080]有益效果:与现有技术比，本发明的有益效果为:首先本发明采用CPU对电力潮流的雅可比矩阵J进行QR符号分解，根据R阵的稀疏格式，可以减少不必要的浮点计算;其次根据R阵的稀疏结构将J阵的各列分到可以并行计算的不同层次，并将分层结果传给GPU;再者在GPU中按层次递增的顺序计算分层QR分解内核函数SparseQR，且GPU线程中采用了归约算法处理点积运算，提高了计算效率;最后本发明利用CPU控制程序的流程并处理基础数据和GPU处理密集的浮点运算相结合的模式提高了电力潮流雅可比矩阵QR分解的效率，解决了电力系统静态安全性分析中潮流计算耗时大的问题。
【附图说明】
:
[0081 ]图1为本发明的实例计算时间；
[0082]图2为本发明的实例分层结果；
[0083]图3为本发明的流程示意图。
【具体实施方式】
:
[0084]如图3所示，本发明一种GPU加速的电力潮流雅可比矩阵的QR分解方法，所述方法包括:
[0085](I)CPU中对雅可比矩阵J进行QR符号分解，得到Household变换矩阵V和上三角矩阵R矩阵的稀疏结构，符号分解之后的J的稀疏结构等于V+R;根据R阵的稀疏结构，对矩阵J各列进行并行化分层。
[0086](2)GPU中按层次递增的顺序计算分层QR分解内核函数SparseQR。
[0087]一、CPU中对电力潮流雅可比矩阵J进行QR符号分解方法
[0088]首先，在CPU中对雅可比矩阵J进行QR符号分解，得到Household变换矩阵V和上三角阵R阵的稀疏结构，符号分解之后的J的稀疏结构等于V+R;然后，并行化分层将矩阵J的η列归并到MaxLevel层中，属于同一层中的列并行进行QR分解；每层包含的列的数量为LeVelnum(k)，k表示层号;存储第k层中所有列号至映射表Mapk。最后，CPU将GPU计算所需数据传输给GPU，数据包括:雅可比矩阵J，其维度η，上三角阵R阵，层数MaxLeveI，每层包含的列数Levelnum以及映射表Map。
[0089]其中QR符号分解原理和并行化分层原理参见
[0090]aDirectMethodsforSparseLinearSystemsTimothyA.Davi s, SI AM,Philadelphia，2006。本专利使用的QR符号分解和并行化分层程序参见CSparse:
[0091 ] aConciseSparseMatrixpackage.VERS10N3.1.4,Copyright(c)2006-2014,T imothyA.Davis，OctlO，2014 Aouseho Id变换原理参见文献:胡冰新，李宁，吕俊.一种采用Household变换递归实现的复矩阵QR分解算法[J].系统仿真学报，2004，16(11):2432-2434。
[0092]二、GPU中按层次递增的顺序启动分层QR分解内核函数SparseQR
[0093]分层QR分解内核函数定义为SparseQIKNblocks，Nthreads>，其线程块大小Nthread固定为128，当对k层进行计算时，线程块数量Nb1cks = Levelnum(k)，总线程数量为:Nb1cks XNthreads;按照层次递增的顺序，调用内核函数SparseQR〈Levelnum(k)，Nthreads〉来分解属于第k层的所有列。
[0094]SparseQR〈Levelnum(k)，Nthreads〉的计算流程为:
[0095](I)CUDA自动为每个线程分配线程块索弓IblockID和线程块中的线程索引threadID；
[0096](2)将blockID和threadID赋值给变量bid和t，之后通过bid和t来索引bid号线程块中的t号线程；
[0097](3)第bid号线程块负责QR分解雅可比矩阵J的第j=Mapk(bid)列；
[0098](4)第bid号线程块中，变量i从I递增到j-Ι，如果R(i，j)#0，执行以下计算:
[0099]首先，计算中间变量0= 2￥(:[:11，；01'.J(1:n，j)，其中V(1:n，i)是Household变换矩阵V中第i列的第i?η行元素构成的向量，J(1:n，j)是雅克比矩阵J中第j列的第i?η行元素构成的向量;然后，采用公式氕1:11，」)=氕1:11，」)-|^(1:11，1)更新雅可比矩阵1的第」列；
[0100]首先，计算中间变量0= 2￥(:[:11，；01'.J(1:n，j)，具体计算步骤为，
[0101]I)判断线程编号t是否小于n-1，否则该线程结束执行；
[0102]2)线程中的变量Temp+ = 2V( i+t，i) X J( i+t，j);
[0103]3)t = t+128，返回 I)执行；
[0104]4)线程中Temp的值赋给共享内存变量Cache[threadID]；
[0105]5)对线程块中的线程进行同步；
[0106]6)变量 M=128/2;
[0107]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0108]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0109]9)对线程块中的线程进行同步；
[0110]10)M/ = 2；
[0111]11)返回7)执行;
[0112]12)P = Cache(0)；
[0113]然后，采用公式J(1: η，j) = J( 1: η，j)-βν(1: η，i)更新雅可比矩阵J的第j列，具体步骤为，
[0114]I)判断线程编号t是否小于n-1，否则线程结束执行；
[0115]2)J(t+i，j)=J(t+i，j)_PV(t+i，i);
[0116]3)t = t+128，返回 I)。
[0117](5)计算第j列的Household变换向量:
[0118]首先，计算变量a2= J( j:n，j)T.J( j:n，j)，具体计算步骤如下:
[0119]I)判断线程编号t是否小于n-j，否则该线程结束执行；
[0120]2)线程中的变量Temp，Temp+ = J( j+t, j) X J( j+t, j);;
[0121]3)t = t+128，返回 I)执行；
[0122]4)线程中Temp的值赋给共享内存变量Cache[threadID];;
[0123]5)对线程块中的线程进行同步；
[0124]6)变量 M= 128/2;
[0125]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0126]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0127]9)对线程块中的线程进行同步；
[0128]10)M/ = 2；
[0129]11)返回7)执行；
[0130]12)a2 = Cache(0)；
[0131 ]其次，计算，v( j:η，j) = J( j:η，j )-aej(j:n)，其中是ej是第j个元素为I的单位向量；
[0132]具体步骤如下:
[0133]I)判断线程编号t是否小于n-j，否则线程结束执行；
[0134]2)V(t+j, j)=J(t+j, j)-aej(t+j)；；
[0135]3)t = t+128，返回I)。
[0136]然后，计算变量b2= V( j:η，j )τ.V( j:η，j);
[0137]具体计算步骤如下:
[0138]I)判断线程编号t是否小于n-j，否则该线程结束执行；
[0139]2)线程中的变量为Temp，Temp+ = V( j+t，j) XV( j+t, j);
[0140]3)t = t+128，返回 I)执行；
[0141 ]4)线程中Temp的值赋给共享内存变量Cache[threadID]；
[0142]5)对线程块中的线程进行同步；
[0143]6)变量 M=128/2;
[0144]7)当M不等于O时执行8);否则结束计算，跳转到12);
[0145]8)判断 thread ID 小于M执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行；
[0146]9)对线程块中的线程进行同步；
[0147]10)M/ = 2；
[0148]11)返回7)执行；
[0149]12)b2 = Cache(0)；
[0150]最后，计算V (j:n, j)=V(j:n, j)/b；
[0151]具体步骤如下:
[0152]I)判断线程编号t是否小于n-j，否则线程结束执行；
[0153]2)V(t+j, j)=V(t+j, j)/b；
[0154]3)t = t+128，返回 I)。
[0155](6)更新J 阵第 j列:J( j，j)=a，J( j+1:n，j)=0。
[0156]本发明所使用的GPU计算平台配备一张TeslaK20CGPU卡和IntelXeonE5-2620CPU，GPU的峰值带宽可达208GB/s，单精度浮点计算量峰值可达3.52Tf lops，CPU主频为2GHz XPU计算平台配备IntelCorei7-3520M2.90GHz的CPU。在CPU和GPU计算平台上分别对16557维稀疏雅可比矩阵实例进行了测试，图1为16557维稀疏雅可比矩阵在不同计算库上的计算时间。在GPU上使用两级并行策略的QR算法求解时间为67ms，落后于KLU的数值分解时间，略快与UMPACK的计算时间。GPU加速算法落后于KLU计算时间的原因主要是每层包含的列数随着层数增加迅速减少，第I层包含1642列，前30层每层大于100列，但到了 207层之后仅剩I?2列。前面层次中可并行计算的列数量多，后面层次中的列需要串行执行，无法充分发挥GPU的计算性能。图2为16557维稀疏雅可比矩阵的并行度变化曲线(消去树每层包含的列数，纵坐标以10为底的对数表示)。
【主权项】
1.一种GRJ加速的电力潮流雅可比矩阵的QR分解方法，其特征在于:所述方法包括: (I)CPU中对雅可比矩阵J进行QR符号分解，得到Household变换矩阵V和上三角阵R阵的稀疏结构;根据R阵的稀疏结构，对矩阵J各列进行并行化分层； (2 )GPU中按层次递增的顺序计算分层QR分解内核函数Spar seQR。2.根据权利要求1所述的GPU加速的电力潮流雅可比矩阵的QR分解方法，其特征在于:所述步骤(I)中，并行化分层将矩阵J的η列归并到MaxLevel层中，属于同一层中的列并行进行QR分解;每层包含的列的数量为Levelnum(k)，k表不层号;存储第k层中所有列号至映射表Mapk ο3.根据权利要求1所述的GPU加速的电力潮流雅可比矩阵的QR分解方法，其特征在于:所述步骤(2)中，分层QR分解内核函数定义为SparseQtKNbIcicks，Nthreads>，其线程块大小^1^固定为128，当对1^层进行计算时，线程块数量他1。。1? = 1^611111111(10，总线程数量为:Nb1cks XNthreads;按照层次递增的顺序，调用内核函数SparseQR〈Levelnum(k)，Nthreads>分解属于第k层的所有列；SparseQR〈Levelnum(k)，Nthreads>的计算流程为: (2.1)0]0六自动为每个线程分配线程块索引1310(^10和线程块中的线程索引访^&(110; (2.2)将blockID和threadID赋值给变量bid和t，之后通过bid和t来索引bid号线程块中的t号线程； (2.3)第1^(1号线程块负责01?分解雅可比矩阵1的第」=1^?1{(1^(1)列； (2.4)第1^(1号线程块中，变量1从1递增到」-1，如果1?(1，」)#0，执行以下计算: 首先，计算中间变量β = 2￥(:[:11，；01'*了(;[:11，」)，其中￥(;[:11，；0是Househo I d变换矩阵V中第i列的第i?η行元素构成的向量，J(1:n，j)是雅克比矩阵J中第j列的第i?η行元素构成的向量;然后，采用公式J(1:n，j)=J(1:n，j)-0V(1:n，i)更新雅可比矩阵J的第j列； (2.5)计算第j列的Household变换向量； (2.6)更新了矩阵第」列:了(」，」)=3，了(」+1:]1，」)=0。4.根据权利要求3所述的GPU加速的电力潮流雅可比矩阵的QR分解方法，其特征在于:所述步骤(2.4)中: 首先，计算中间变量0 = 2￥(:[:11，；01'.J(1:n，j)，具体计算步骤为， 1)判断线程编号t是否小于n-1，否则该线程结束执行； 2)线程中的变量1611^1+= 2￥(1+1:，；0 X J(i+t，j); 3)t= t+128，返回 I)执行； 4)线程中Temp的值赋给共享内存变量Cache[threadID]； 5)对线程块中的线程进行同步； 6)变量M= 128/2; 7)当M不等于O时执行8);否则结束计算，跳转到12); 8)判断thread ID 小于 M 执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行； 9)对线程块中的线程进行同步；10)M/= 2； 11)返回7)执行；12)P= Cache(0)；然后，采用公式J(1: η，j) = J(1: η，j)-βν(1: η，i)更新雅可比矩阵J的第j列，具体步骤为， 1)判断线程编号t是否小于n-1，否则线程结束执行； 2)J(t+i，j)=J(t+i，j)-0V(t+i，i); 3)t= t+128，返回 I)。5.根据权利要求3所述的GPU加速的电力潮流雅可比矩阵的QR分解方法，其特征在于:所述步骤(2.5)中， 首先，计算变量a2 = J( j:n，j)T.J( j:n，j)，具体计算步骤如下: 1)判断线程编号t是否小于n-j，否则该线程结束执行； 2)线程中的变量Temp，Temp+= J( j+t, j) X J( j+t, j)；； 3)t= t+128，返回 I)执行； 4)线程中Temp的值赋给共享内存变量Cache[threadID]；； 5)对线程块中的线程进行同步； 6)变量M= 128/2; 7)当M不等于O时执行8);否则结束计算，跳转到12);8)判断thread ID 小于 M 执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行； 9)对线程块中的线程进行同步；10)M/= 2； 11)返回7)执行； 12)a2= Cache(0)； 其次，计算，V( j:n，j)=J( j:n，j)-aej( j:n)，其中是ej是第j个元素为I的单位向量； 具体步骤如下: 1)判断线程编号t是否小于n-j，否则线程结束执行； 2)V(t+j，j)=J(t+j，j)_aej(t+j);； 3)t= t+128，返回 I)。 然后，计算变量b2 = V( j:η，j )τ.V( j:η，j); 具体计算步骤如下: 1)判断线程编号t是否小于n-j，否则该线程结束执行； 2)线程中的变量为Temp，Temp+= V( j+t，j) XV( j+t, j); 3)t= t+128，返回 I)执行； 4)线程中Temp的值赋给共享内存变量Cache[threadID]； 5)对线程块中的线程进行同步； 6)变量M= 128/2; 7)当M不等于O时执行8);否则结束计算，跳转到12);8)判断thread ID 小于 M 执行:Cache [threadID] + = Cache [threadID+M]，否则线程不执行； 9)对线程块中的线程进行同步；10)M/= 2；11)返回7)执行；12)b2= Cache(0)；最后，计算 V (j: η，j) = V( j: η，j )/b;具体步骤如下:1)判断线程编号t是否小于n-j，否则线程结束执行；2)V(t+j，j)=V(t+j，j)/b;3)t= t+128，返回 I)。
【文档编号】H02J3/06GK106026107SQ201610592223
【公开日】2016年10月12日
【申请日】2016年7月26日
【发明人】周赣, 孙立成, 张旭, 柏瑞, 冯燕钧, 秦成明, 傅萌
【申请人】东南大学