本发明涉及电力系统分析技术领域，具体地，涉及一种基于bpa和matlab的大型受端电网轨迹灵敏度批处理计算方法。
背景技术：
轨迹灵敏度指标可用于衡量电力系统某节点上安装动态无功补偿设备后系统暂态电压稳定的提升程度。对于大型受端电网中的某一负荷节点j，它的轨迹灵敏度tsi可采用下列计算式：
式中：nl为考察故障总数；wfl为故障fl的权重，fl(l＝1,2,...,nl)为考察故障；wdi为母线i的权重，n为负荷节点总数；nk为总的时间采样点数；wk为tk时刻的权重；vi,l(tk,qj0+△qj(comp))为节点j安装额定容量为scomp的动态无功补偿设备后，在故障fl下，节点电压vi在t＝tk时刻的值；△qj(comp)为在节点j安装动态无功补偿设备后，该设备释放的无功功率；vi,l(tk,qj0)为未安装动态无功补偿设备情况下，在故障fl下，vl在t＝tk时刻的值。
可以看出，计算j节点的轨迹灵敏度，关键是需要获得整个电力系统中所有负荷节点在j节点是否安装动态无功补偿设备两种情况下，每种故障发生后的节点电压时域分布。
对于电网节点故障电压的计算和仿真，国内外有不少成熟的软件。中国电力科学院开发的机电暂态软件包psd-bpa(简称bpa)主要由潮流程序和稳定程序构成，具有计算规模大、计算速度快、数值稳定性好、功能强等特点，已在我国电力系统规划、调度、生产运行及科研部门得到了广泛应用。
利用bpa对大型受端电网的负荷节点进行轨迹灵敏度计算时，需要首先根据网络数据通过潮流程序得到系统所处运行方式下的潮流结果文件(.bse格式)，然后在稳定程序中导入该潮流结果文件，并手动填写故障参数卡片和动态无功补偿设备参数卡片。在计算受端电网所有负荷节点的轨迹灵敏度从而寻找全网最优动态无功补偿设备安装点时，需要遍历每一种典型故障与补偿设备安装点(也就是负荷节点)的组合，利用bpa稳定程序对每一种组合进行稳定计算，并在每个组合的计算中得到所有负荷节点上的电压变化时域分布。对于大型受端电网，其典型故障数与负荷节点数都具有一定规模，手动修改故障参数卡片和动态无功补偿设备参数卡片将大大增加工作量，降低计算效率；并且稳定程序的输出数据分散在不同的.swx格式结果文件中，没有统一的命名管理，不利于统计分析。
技术实现要素：
针对现有技术中存在的上述不足，本发明的目的是提供一种基于bpa和matlab的大型受端电网轨迹灵敏度批处理计算方法，该方法利用matlab脚本遍历每种故障与负荷节点的组合，自动修改bpa稳定程序所用到的故障参数卡片与动态无功补偿设备卡片，从而进行稳定计算；并自动对输出数据文件按规律进行命名管理，实现对每种组合电压变化曲线的自动读取，提高计算效率，最终得到受端电网每一个负荷节点上的轨迹灵敏度。
为了达到上述目的，本发明采用如下技术方案。
一种大型受端电网轨迹灵敏度批处理计算方法，基于bpa和matlab实现，包括如下步骤：
s1.利用matlab根据bpa潮流文件所提供的受端电网的网络数据分别生成故障集文件、动态无功补偿设备安装点集文件和负荷节点电压输出语句集文件，并记录下故障个数nf和动态无功补偿设备安装点个数nc。
优选地，上述文件均为txt格式。其中，故障集文件包括所有典型故障对应的bpa稳定程序卡片；动态无功补偿设备安装点集文件包括在所要计算轨迹灵敏度指标的所有节点上安装动态无功补偿设备的bpa稳定程序卡片；负荷节点电压输出语句集文件包括全网所有负荷节点上电压时域分布信息的输出指令语句。
s2.根据s1生成的三个文件(故障集文件、动态无功补偿设备安装点集文件和负荷节点电压输出语句集文件)，形成原始的.swi格式稳定程序输入数据文件，并记录下故障参数卡和动态无功补偿设备参数卡的起始行信息。
优选地，上述原始的.swi格式稳定程序输入数据文件中，既包括所要遍历所有组合间相同的部分，如控制数据卡、发电机模型卡、励磁系统模型卡、电力系统稳定器(pss)模型卡、调速器和原动机模型卡、直流控制系统模型卡、负荷模型卡、输出控制卡等；也包括不同组合间产生变化的部分，如故障操作模型卡、动态无功补偿(电力电子)设备模型卡。所有故障操作模型卡和动态无功补偿设备模型卡均需要打上注释符号，以便求解到某种组合时能以去掉注释符号的方式形成相应的.swi格式稳定程序输入数据文件。
s3.以批处理的方式遍历nf×nc种不同故障搭配不同动态无功补偿设备安装点的组合，根据稳定程序输出结果文件得到每种组合下全网所有负荷节点的电压时域分布(每种组合下全网所有负荷节点的电压时域分布记录在.swx格式的稳定程序输出结果文件中)。
优选地，包括如下子步骤：
s3.1，定义故障编号变量faultnum和动态无功补偿设备安装点编号变量compnum，利用matlab的for循环功能遍历不同故障与安装点的组合。其中，变量faultnum和compnum的格式均为整数。
s3.2，利用matlab读取原始的.swi格式稳定程序输入数据文件，根据故障编号变量faultnum、动态无功补偿设备安装点编号变量compnum以及之前记录下的故障参数卡、补偿设备参数卡起始行信息，删除相应故障操作模型卡和动态无功补偿设备模型卡首列的注释符号，使故障操作模型卡和动态无功补偿设备模型卡生效。
s3.3，利用matlab生成每种组合中稳定计算需用到的.swi格式稳定程序输入数据文件，并统一按照fault_faultnum_comp_compnum的规则命名。文件名fault_faultnum_comp_compnum的格式为字符串，其中，变量faultnum和compnum的格式由整数通过num2str函数转换成字符串，再与固定字符串fault_和_comp_组合。
按照上述规则统一命名有两个好处：一是方便之后调用稳定程序时传递稳定程序输入数据文件(.swi文件)；二是由于稳定程序输出结果文件(.swx文件)的名称与稳定程序输入数据文件相同，所以统一命名能方便读取稳定计算结果。
s3.4，在matlab里利用dos命令调用bpa稳定程序(swnt.exe)，并将s3.3中生成的稳定程序输入数据文件(.swi文件)作为传递参数导入bpa稳定程序中，执行稳定计算。在所遍历的每种组合中执行稳定计算将得到相应的全网所有负荷节点的电压时域分布结果，这些结果以所遍历的组合为单位，记录在不同的稳定程序输出结果文件(.swx文件)中。由于s3.3的操作，所得稳定程序输出结果文件将按fault_faultnum_comp_compnum的命名规则排列。
s4.利用matlab读取s3中的稳定程序输出结果文件(.swx文件)，得到当节点j(j＝1,2,...,n)安装动态无功补偿设备后，各负荷节点在不同故障情况下的电压时域分布。并根据权重系数加权求和，得到中间变量r1j。r1j的计算式如下式所示：
式中，nl为考察故障总数；nk为总的时间采样点数；wfl为故障fl的权重；wdi为母线i的权重；wk为tk时刻的权重；△qj(comp)为在节点j安装动态无功补偿设备后，该设备释放的无功功率；vi,l(tk,qj0+△qj(comp))为节点j安装额定容量为scomp的动态无功补偿设备后，在故障fl下，节点电压vi在t＝tk时刻的值。
s5.不考虑安装动态无功补偿设备，以批处理的方式遍历nf种不同故障的情况，得到每种情况下全网所有负荷节点的电压时域分布，并记录在稳定程序输出结果文件中。
优选地，包括如下子步骤：
s5.1，根据s3.1中定义的故障编号变量faultnum，利用matlab的for循环功能遍历不同故障情况。
s5.2，利用matlab读取原始的.swi格式稳定程序输入数据文件，根据faultnum以及之前记录下的故障参数卡起始行信息，删除相应故障操作模型卡首列的注释符号，使故障操作模型卡生效。
s5.3，利用matlab生成每种故障情况中稳定计算需用到的.swi格式稳定程序输入数据文件，并统一按照no_comp_fault_faultnum的规则命名。文件名no_comp_fault_faultnum的格式为字符串，其中，变量faultnum的格式由整数通过num2str函数转换成字符串，再与固定字符串no_comp_fault_组合。
s5.4，在matlab里利用dos命令调用bpa稳定程序(swnt.exe)，并将s5.3生成的稳定程序输入数据文件(.swi文件)作为传递参数导入bpa稳定程序中，执行稳定计算。在所遍历的每种故障情况中执行稳定计算将得到相应的全网所有负荷节点的电压时域分布结果，这些结果以所遍历的故障情况为单位，记录在不同的稳定程序输出结果文件(.swx文件)中。由于s5.3的操作，所得稳定程序输出结果文件将按no_comp_fault_faultnum的命名规则排列。
s6.利用matlab读取s5中的稳定程序输出结果文件(.swx文件)，得到不安装动态无功补偿设备时，各负荷节点在不同故障情况下的电压时域分布。并根据权重系数加权求和，得到中间变量r2。r2的计算式如下式所示：
式中，vi,l(tk,qj0)为未安装动态无功补偿设备情况下，在故障fl下，节点电压vl在t＝tk时刻的值。
s7.由下式计算得到负荷节点j上的轨迹灵敏度tsij：
tsij＝r1j-r2。
与现有技术相比，本发明具有如下有益效果：
现有的手动计算大型受端电网所有负荷节点的轨迹灵敏度指标涉及到很大的工程量，而采用本发明提供的方法，所有过程均可由matlab脚本实现，自动调用bpa稳定程序进行批处理计算，避免重复而枯燥的手动操作，不仅大大增加了计算效率，还避免了出错。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1是基于bpa和matlab的大型受端电网轨迹灵敏度批处理计算方法流程图；
图2是matlab利用dos命令调用bpa稳定程序，并传递输入数据文件示意图。
具体实施方式
下面对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。
实施例
本实施例提供了一种基于bpa和matlab的大型受端电网轨迹灵敏度批处理计算方法，其流程如图1所示，包括下列步骤：
第s1步.利用matlab根据潮流文件所提供的受端电网的网络数据生成故障集文件、动态无功补偿设备安装点集文件和负荷节点电压输出语句集文件，并记录下故障个数nf和动态无功补偿设备安装点个数nc。
具体地：
第s1.1步，在matlab中读取.dat格式的bpa潮流文件，并识别文件中的支路数据卡，提取支路首末端节点名称及其基准电压信息，据此生成flt支路三相短路接地故障卡。将所有故障卡依次排列，组成.txt格式的故障集文件。
第s1.2步，在matlab中读取.dat格式的bpa潮流文件，并识别文件中的节点数据卡，提取负荷节点的节点名称及其基准电压信息，据此生成动态无功补偿设备安装卡。以statcom为例，可生成多张不同安装点的vg/vg+卡。将其依次排列，组成.txt格式动态无功补偿设备安装点集文件。
第s1.3步，根据第s1.2步中提取的负荷节点的节点名称及其基准电压信息，生成bh/b母线输出卡。将其依次排列，组成.txt格式负荷节点电压输出语句集文件。
第s2步.根据第1步生成的三个文件，形成原始的.swi格式稳定程序输入数据文件，并记录下故障(参数)卡和动态无功补偿设备(参数)卡的起始行。
具体地：利用matlab读写文件的功能，将第s1步生成的三个文件中的卡片(故障操作模型卡和动态无功补偿设备模型卡)统一在列首打上注释符号，写入同一文件中，并加上所要遍历的不同组合间相同的部分，如控制数据卡、发电机模型卡、励磁系统模型卡、电力系统稳定器(pss)模型卡、调速器和原动机模型卡、直流控制系统模型卡、负荷模型卡、输出控制卡等，生成原始的.swi格式稳定程序输入数据文件。记录下所生成的.swi格式文件中，flt卡的起始行startline_f和vg/vg+卡的起始行startline_c。
第s3步.以批处理的方式遍历nf×nc种不同故障搭配不同动态无功补偿设备安装点的组合，得到每种组合下全网所有负荷节点的电压时域分布。
具体地：
第s3.1步，定义故障编号变量“faultnum”和动态无功补偿设备安装点编号变量“compnum”，利用matlab的for循环功能遍历faultnum×compnum种不同故障与安装点的组合。
第s3.2步，在每一个循环中，利用matlab读取原始的.swi格式的稳定程序输入数据文件，删除原始文件中第startline_f+faultnum-1行和第startline_c+(compnum-1)*2行首列的注释符号，使相应卡片生效。
第s3.3步，将第s3.2步中修改后的.swi格式文件统一按照“fault_faultnum_comp_compnum”的规则重命名，生成每种组合中稳定计算需用到的.swi格式稳定程序输入数据文件。由于文件名称为字符串格式，而faultnum和compnum的格式为整数，故需要通过num2str函数将数值变量faultnum和compnum转换成的字符串，再与固定字符串fault_和_comp_组合成文件名。
“faultnum”和“compnum”在最初定义的时候为数值类型的变量(数据格式为整数)，用来给故障和补偿设备安装点编号。本步骤中，“fault_faultnum_comp_compnum”为字符串格式，在命名的时候，需要将整数型的数值变量(faultnum和compnum)转换成字符串，于是用到num2str函数将数值变量faultnum和compnum转换成字符串。num2str函数是matlab函数库里边现成的函数，直接调用即可。
按照上述规则统一命名有两个好处：一是方便之后调用稳定程序时可根据当前循环中faultnum和compnum的值传递输入数据文件(.swi文件)；二是由于稳定程序输出结果文件(.swx文件)的名称与输入数据文件(.swi文件)相同，所以统一命名能方便读取稳定计算结果。
第s3.4步，在matlab里利用dos命令调用bpa稳定程序(swnt.exe)，并将第s3.3步中生成的稳定程序输入数据文件(.swi文件)作为传递参数导入稳定程序中，执行稳定计算，如图2所示。图中，path_swnt为字符串格式的bpa稳定程序的路径；path_workspace为字符串格式的当前matlab脚本所在文件夹路径；
filename为第s3.3步中所生成的.swi格式稳定程序输入数据文件名。在所遍历的每种组合中执行稳定计算将得到相应的全网负荷节点的电压时域分布结果，这些结果以所遍历的组合为单位，记录在不同的.swx稳定计算结果文件(稳定程序输出结果文件)中。由于第s3.3步的操作，所得计算结果文件将按“fault_faultnum_comp_compnum”的命名规则排列。
在步骤s3中，根据稳定程序输出结果文件得到每种组合下全网所有负荷节点的电压时域分布。在稳定程序输入数据文件(.swi文件)填好后执行稳定计算，系统会生成稳定程序输出结果文件(.swx文件)，该.swx文件中具有当前条件(组合)下全网所有负荷节点的电压时域分布。
第s4步.利用matlab读取第s3.4步中的bpa稳定程序输出结果文件(.swx文件)，得到当节点j(j＝1,2,...,n)安装动态无功补偿设备后，各负荷节点在不同故障情况下的电压时域分布。并根据权重系数加权求和，得到中间变量r1j。r1j的计算式如下式所示：
第s5步.不考虑安装动态无功补偿设备，以批处理的方式遍历nf种不同故障的情况，得到每种情况下全网所有负荷节点的电压时域分布。
具体地：
第s5.1步，根据第s3步中定义的故障编号变量“faultnum”，利用matlab的for循环功能遍历不同故障情况。
第s5.2步，在每一个循环中，利用matlab读取原始的.swi格式稳定程序输入数据文件，删除原始文件中第startline_f+faultnum-1行首列的注释符号，使相应故障卡生效。
第s5.3步，将第s5.2步中修改后的.swi格式文件统一按照“no_comp_fault_faultnum”的规则重命名，生成每种故障情况中稳定计算需用到的.swi格式稳定程序输入数据文件。由于文件名称为字符串格式，而faultnum的格式为整数，故需要通过num2str函数将数值变量faultnum转换成的字符串，再与固定字符串no_comp_fault_组合成文件名。
第s5.4步，在matlab里利用dos命令调用bpa稳定程序(swnt.exe)，并将第s5.3步生成的稳定程序输入数据文件(.swi文件)作为传递参数导入稳定程序中，执行稳定计算。调用过程除了filename为第s5.3步中所生成的.swi格式稳定程序输入数据文件名，其余与第s3.4步相似。在所遍历的每种故障情况中执行稳定计算将得到相应的全网负荷节点的电压时域分布结果，这些结果以所遍历的故障情况为单位，记录在不同的.swx稳定计算结果文件(稳定程序输出结果文件)中。由于第s5.3步的操作，所得计算结果文件将按“no_comp_fault_faultnum”的命名规则排列。
第s6步.利用matlab读取第s5.4步中的bpa稳定程序输出结果文件(.swx文件)，得到不安装动态无功补偿设备时，各负荷节点在不同故障情况下的电压时域分布。并根据权重系数加权求和，得到中间变量r2。r2的计算式如下式所示：
第s7步.由下式计算得到负荷节点j上的轨迹灵敏度tsij：
tsij＝r1j-r2。
由于大型受端电网中不止一个负荷节点，所以最终可得到一个n维数组tsi[n]，用来记录大型受端电网中所有负荷节点上的轨迹灵敏度。
需要说明的是：上述对大型受端电网轨迹灵敏度批处理计算方法的描述，不受执行步骤的限制，根据实际应用情况，步骤顺序可以调整。
本实施例提供的大型受端电网轨迹灵敏度批处理计算方法，利用matlab脚本遍历每种故障与负荷节点的组合，自动修改bpa稳定程序所用到的故障参数卡片与动态无功补偿设备卡片，从而进行稳定计算；并自动对输出数据文件按规律进行命名管理，实现对每种组合电压变化曲线的自动读取，提高计算效率，最终得到受端电网每一个负荷节点上的轨迹灵敏度。现有的手动方法计算大型受端电网所有负荷节点的轨迹灵敏度指标涉及到很大的工程量，而采用本实施例的好处在于提供的方法，所有过程均可由matlab脚本实现，自动调用bpa稳定程序进行批处理计算，避免重复而枯燥的手动操作，不仅大大增加了计算效率，还避免了出错。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。