电力系统分析PQ分解法法潮流计算的MATLAB程序

%PQ分解法潮流计算 %可以使用的元件模型：PQ节点、PV节点、平衡节点、输电线路、双绕组变压器、接地阻抗 %在运行此程序前，确保文件"system_parameters"和"flowpq_output"与此程序的源文件在同一个路径下 %在"system_parameters"中从每张工作表的第二行开始，按照格式输入系统参数，保存后运行此程序，计算的结果将存放在"flowcalc1_output"中 %为了保证每次计算后存放结果的表格中多余的内容能够被全部刷新，在运行此程序前尽量确保"flowcalc1_output"中的每张工作表的"A1"单元格都是非空的，或者整张工作表都是空的(若未人为编辑过该表格则无须注意) %所有数据均采用标幺值，允许输入初等算式，节点的功率正方向以节点流向系统为正，例如负荷节点的功率一般为负值，发电机节点的功率一般为正值 %若两个节点之间存在变压器对地支路，在输入系统参数时应当把与对地支路之间不经过阻抗直接相连的节点的序号填入变压器对地支路表格第一列的单元格里，另一个节点序号填入第二列 %若某节点i单独连接了对地导纳y，可以在"system_parameters"的变压器对地支路表格中添加一行[i i y]，则运行程序后输出的支路功率分布数据中第i行第i列表示从节点i流向大地的复功率 %本程序并未考虑节点的功率和电压上下限，若存在限制，可以先运行一遍，根据结果把达到上下限的节点在下次运行时更改为其他类型 %输电线路采用π型等值电路，变压器采用τ型等值电路，两个节点之间可以并联一个变压器和多条输电线路，但不可并联变压器 %如果某种类型的节点或支路数量为0，则无需输入 %本程序中电压相角以弧度制参与计算和收敛判断，但运行程序后输出的相角均为角度制 %% 输入并整理系统参数 filename='system_parameters.xlsx';%系统参数文件名称 %输电线路 [~,~,temp]=xlsread(filename,'输电线路参数'); n_line=size(temp,1)-1;%输电线路条数 transmission_line=zeros(n_line,4);%输电线路参数[从(序号) 到(序号) 阻抗 总导纳;] for a=1:n_line for b=1:4 if isa(temp{a+1,b},'double') transmission_line(a,b)=temp{a+1,b}; else%复数在excel中是字符串，要先转化为数值 transmission_line(a,b)=str2num(temp{a+1,b}); end end end transmission_line=rmmissing(transmission_line);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 n_line=size(transmission_line,1);%输电线路条数 %变压器串联支路 [~,~,temp]=xlsread(filename,'变压器串联支路参数'); n_series=size(temp,1)-1;%变压器串联支路条数 transformer_series=zeros(n_series,4);%变压器串联支路参数[从(序号) 到(序号) 变比(一次/二次) 阻抗;] for a=1:n_series for b=1:4 if isa(temp{a+1,b},'double') transformer_series(a,b)=temp{a+1,b}; else%excel中的数据可能是字符串形式，要先转化为数值 transformer_series(a,b)=str2num(temp{a+1,b}); end end end transformer_series=rmmissing(transformer_series);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 n_series=size(transformer_series,1);%变压器串联支路条数 %变压器对地支路 [~,~,temp]=xlsread(filename,'变压器对地支路参数'); n_shunt=size(temp,1)-1;%变压器对地支路条数 transformer_shunt=zeros(n_shunt,3);%变压器对地支路参数[从(序号) 到(序号) 总导纳;] for a=1:n_shunt for b=1:3 if isa(temp{a+1,b},'double') transformer_shunt(a,b)=temp{a+1,b}; else%复数在excel中是字符串，要先转化为数值 transformer_shunt(a,b)=str2num(temp{a+1,b}); end end end transformer_shunt=rmmissing(transformer_shunt);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 n_shunt=size(transformer_shunt,1);%变压器对地支路条数 %% 输入并整理节点信息 [~,~,temp]=xlsread(filename,'PQ节点'); n_PQ=size(temp,1)-1;%PQ节点个数 given_PQ=zeros(n_PQ,3);%PQ节点参数[序号 注入有功 注入无功;] for a=1:n_PQ for b=1:3 if isa(temp{a+1,b},'double') given_PQ(a,b)=temp{a+1,b}; else%复数在excel中是字符串，要先转化为数值 given_PQ(a,b)=str2num(temp{a+1,b}); end end end given_PQ=rmmissing(given_PQ);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 PQ_number=size(given_PQ,1);%PQ节点个数 [~,~,temp]=xlsread(filename,'PV节点'); n_PV=size(temp,1)-1;%PQ节点个数 given_PV=zeros(n_PV,3);%PQ节点参数[序号 注入有功 电压幅值;] for a=1:n_PV for b=1:3 if isa(temp{a+1,b},'double') given_PV(a,b)=temp{a+1,b}; else%复数在excel中是字符串，要先转化为数值 given_PV(a,b)=str2num(temp{a+1,b}); end end end given_PV=rmmissing(given_PV);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 PV_number=size(given_PV,1);%PV节点个数 [~,~,temp]=xlsread(filename,'平衡节点'); n_slack=size(temp,1)-1;%平衡节点个数 given_slack=zeros(n_slack,3);%平衡节点参数[序号 电压幅值 电压相角;] for a=1:n_slack for b=1:3 if isa(temp{a+1,b},'double') given_slack(a,b)=temp{a+1,b}; else%复数在excel中是字符串，要先转化为数值 given_slack(a,b)=str2num(temp{a+1,b}); end end end given_slack=rmmissing(given_slack);%若用户不小心更改了表格的单元格格式可能会出现nan,将其删除 slack_number=size(given_slack,1);%平衡节点个数 number=PQ_number+PV_number+slack_number;%节点总数 i_PQ=[]; i_PV=[]; i_slack=[]; if ~isempty(given_PQ)%若系统有PQ节点 i_PQ=given_PQ(:,1);%PQ节点序号 end if ~isempty(given_PV)%若系统有PV节点 i_PV=given_PV(:,1);%PV节点序号 end if ~isempty(given_slack)%若系统有平衡节点 i_slack=given_slack(:,1);%平衡节点序号 end real_power=zeros(PQ_number+PV_number,1); image_power=zeros(PQ_number,1); if ~isempty(given_PQ)%若系统中有PQ节点 real_power(1:PQ_number,1)=given_PQ(:,2);%PQ节点的有功按输入的顺序放入有功列向量中 image_power=given_PQ(:,3);%PQ节点的无功按输入的顺序放入无功列向量里 end if ~isempty(given_PV)%若系统中有PV节点 real_power(PQ_number+1:PQ_number+PV_number,1)=given_PV(:,2);%PV节点的有功按输入的顺序接着PQ节点放入有功列向量中 end %% 计算导纳矩阵 %生成节点导纳矩阵、对地导纳矩阵和串联支路导纳矩阵，后两者用于计算功率分布，不需考虑理想变压器，因为后续的电压矩阵V_matrix会考虑。Y0第i行第j列表示ij之间接在i上的导纳。 Y=zeros(number);%节点导纳矩阵 Y0=zeros(number);%对地导纳矩阵 Y1=zeros(number);%串联支路导纳矩阵 %将输电线路参数加入节点导纳矩阵、对地导纳矩阵和串联支路导纳矩阵 for a=1:n_line b1=transmission_line(a,1);%首端节点 b2=transmission_line(a,2);%末端节点 y_imp=1/transmission_line(a,3);%线路阻抗对应导纳值 y_con=transmission_line(a,4)/2;%线路导纳对应导纳值 %节点导纳矩阵 %-非对角元素 Y(b1,b2)=Y(b1,b2)-y_imp;%加上负的线路阻抗的倒数 Y(b2,b1)=Y(b2,b1)-y_imp;%对称位置值相等 %-对角元素 Y(b1,b1)=Y(b1,b1)+y_imp+y_con;%加上线路阻抗的倒数和二分之一线路导纳 Y(b2,b2)=Y(b2,b2)+y_imp+y_con;%输电线路是对称元件，所以对两端节点对应元素增加量一致 %对地导纳矩阵 Y0(b1,b2)=Y0(b1,b2)+y_con; Y0(b2,b1)=Y0(b2,b1)+y_con; %串联支路导纳矩阵 Y1(b1,b2)=Y1(b1,b2)+y_imp; Y1(b2,b1)=Y1(b2,b1)+y_imp; end %将变压器串联支路参数加入节点导纳矩阵 for a=1:n_series b1=transformer_series(a,1);%始端节点 b2=transformer_series(a,2);%末端节点 k=transformer_series(a,3);%变比（始端/末端） y_imp=1/transformer_series(a,4);%阻抗 %节点导纳矩阵 Y([b1,b2],[b1,b2])=Y([b1,b2],[b1,b2])+[1/k -1;-1 k]*y_imp/k; %串联支路导纳矩阵 Y1(b1,b2)=Y1(b1,b2)+y_imp; Y1(b2,b1)=Y1(b2,b1)+y_imp; end %将变压器对地支路参数加入节点导纳矩阵 for a=1:n_shunt b1=transformer_shunt(a,1);%所属变压器首端节点 b2=transformer_shunt(a,2);%所属变压器末端节点 y_con=transformer_shunt(a,3);%导纳值 n1=ismember(transformer_series(:,(1:2)),[b1,b2],'rows');%理想变压器在对地支路同侧的变压器串联支路在参数表中的行数 %节点导纳矩阵 if any(n1)%若节点经理想变压器连接对地导纳 Y(b1,b1)=Y(b1,b1)+y_con/transformer_series(n1,3)^2; else%对地导纳直接接在节点上 Y(b1,b1)=Y(b1,b1)+y_con; end %对地导纳矩阵 Y0(b1,b2)=Y0(b1,b2)+y_con; end disp('节点导纳矩阵') disp(Y) %% 定义列向量形式的功率方程和B矩阵 %-有功功率方程 输入(要计算有功的节点序号向量,全部节点电压向量[幅值 相角;]) real_equation=@(i,V)V(i,1).*sum(ones(size(i,1),1)*V(:,1)'.*(cos(V(i,2)*ones(1,number)-ones(size(i,1),1)*V(:,2)').*real(Y(i,:))+sin(V(i,2)*ones(1,number)-ones(size(i,1),1)*V(:,2)').*imag(Y(i,:))),2); %-无功功率方程 输入(要计�