《电力系统稳态分析》课程设计.doc

电力系统稳态分析课程设计学生姓名： 郑立娟 学 号： 010800734 专业班级： 电力系统2008级7班 指导教师： 刘丽军 二 一一 年 十二 月 四 日目 录1、 设计任务2、 牛顿拉夫逊算法原理和步骤3、 牛顿拉夫逊算法和P-Q解耦法的比较4、 牛顿拉夫逊算法程序流程图5、 MATLAB程序6、 运行结果一、设计任务 ；。1,2,3号节点为PQ节点，4号为PV节点，5号为平衡节点，；允许误差选用牛顿拉夫逊直角坐标算法计算上图九节点网络中各节点的电压和相角以及各支路的潮流分布。二、牛顿拉夫逊潮流算法原理及计算步骤牛顿型潮流计算的核心问题是修正方程式的建立和求解。为说明这一修正方程式的建立过程，先对网络中各类节点的编号作如下约定：（1） 网络中共有n个节点，编号为1，2，3.,n其中包含一个平衡节点，编号为s;（2） 网络中有（m-1）个节P-Q点，编号为1，2，3.,m，其中包含编号为s的平衡节点；（3） 网络中有n-m个PV节点，编号为m+1,m+2,.,n. 据此，节点电压用直角坐标表示时，可建立如下修正方程式： 式中的、分别为为注入功率和节点电压平方的不平衡量。式中雅克比矩阵的各个元素则分别为为求取这些偏导数，可将、分别展开如下时，由于对特定的j，只有该特定节点的和是变量，当i=j时牛顿拉夫逊潮流算法计算基本步骤（1） 形成节点矩阵；（2） 设各节点电压的初值、和；（3） 将各节点电压的初值代入公式，求修正方程式中的不平衡量（4） 将各节点电压的初值带入公式，求修正方程式的系数矩阵雅可比矩阵的各个元素（5） 解修正方程式，求各节点电压的变化量，即修正量（6） 计算各节点电压的新值，即修正后的值 (7)运用各节点电压的新值自第三步开始进入下一次迭代； （8）计算平衡节点功率和线路功率。三、牛顿拉夫逊算法和P-Q解耦法的比较P-Q分解法的修正方程式具有如下的特点：（1） 以一个（n-1）阶和一个（m-1）阶系数矩阵替代原有的（n+m-2）阶矩阵J，提高了计算速度，降低了大队贮存容量的要求；（2） 以迭代过程中保持不变的系数矩阵替代起变化的系数矩阵J，显著提高了计算速度；（3） 以对称的系数矩阵替代不对称的系数矩阵J,使求逆等运算和所需的储存容量都大大减少；但应强调指出，到处这修正方程式式所作的种种简化毫不影响用这种方法计算的精确度。因采用这种方法时，迭代收敛的判据仍是。 一般情况下采用P-Q分解法计算时要求的迭代次数较采用牛顿拉夫逊法时多，但每次迭代的时间较较采用牛顿拉夫逊法时少，以致总的计算速度仍是P-Q分解法快。四、牛顿拉夫逊算法程序流程图启动输入原始数据 形成节点导纳矩阵YB设节点电压ei(0),fi(0);i=1,2,n,is置迭代次数 k=0对PQ节点:计算Pi(k),Qi(k)对PV节点:计算Pi(k),Ui(k)2置节点号i=1雅克比矩阵J是否已全部形成,in?计算雅克比矩阵元素 否 是 增大节点号 ii+1 增大迭代次数,kk+1解修正方程式,求各节点电压变量ei(k),fi(k);i=1,2,n,is计算各节点电压新值求出|e(k)|max,|f (k)|max迭代是否收敛, |e(k)|max,|f (k)|max? 否计算平衡节点功率和线路功率 是 停止 五、MATLAB程序 function bus_res,S_res = PowerFlow_NR % 牛顿-拉夫逊法解潮流方程的主程序% % % bus = 1 1.05 0 0 0 3 ; 2 1 0 -2 -1 1 ; 3 1 0 -1.8 -0.4 1 ; 4 1 0 -1.6 -0.8 1 ; 5 1 0 -3.7 -1.3 1 ; 6 1.05 0 5 0 2; line = 2 1 0 0.03 0 0 1.05 ;2 3 0.06 0.025 0 0 0 ;3 4 0.08 0.3 0 0 0 ; 4 5 0.1 0.35 0 0 0;6 5 0 0.015 0 0 1.05;2 5 0.04 0.25 0 0 0 ; % 读取bus（节点数据）和line（线路数据）nb,mb=size(bus);nl,ml=size(line); % 读取bus和line的行数和列数bus,line,nPQ,nPV,nodenum = Num_(bus,line); % 对bus和line重新排列 Y = y_(bus,line); % 计算节点导纳矩阵myf = fopen(Result.m,w);fprintf(myf,- by longdinhohe - -nn);fclose(myf); % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵format longEPS = 0.00001; % 设定误差精度for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出OrgS = Or_(Y,bus,nb); % 计算初始功率参数DetaS = dPQ_(Y,bus,nb,nPQ,OrgS); % 计算功率偏差dP和dQ的子程序 J = Jac_(bus,Y,OrgS); % 计算雅克比矩阵的子程序DetaU=qiuni(J,DetaS); % 计算修正量j=0; % 计算新值for i=1:nb %对PQ节点处理 if bus(i,6)=1 j=j+1; bus(i,2)=bus(i,2)+DetaU(2*j,1)+DetaU(2*j-1,1)*(sqrt(-1); endendfor i=1:nb %对PV节点的处理 if bus(i,6)=2 j=j+1; bus(i,2)=bus(i,2)+DetaU(2*j,1)+DetaU(2*j-1,1)*(sqrt(-1); endend if (max(abs(DetaU) 3 k4=k-1; for k3=3:k4% 用k3行从第三行开始到当前行的前一行k4行消去 for k2=k1:N1% k3行后各行上三角元素 b(k3,k2)=b(k3,k2)-b(k3,k)*b(k,k2);%消去运算（当前行k列元素消为0） end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素 b(k3,k)=0; %当前行第k列元素已消为0 end if k=N0 %若已到最后一行 break; end %= 前代运算 = for k3=k1:N0 % 从k+1行到2*n最后一行 for k2=k1:N1 % 从k+1列到扩展列消去k+1行后各行下三角元素 b(k3,k2)=b(k3,k2)-b(k3,k)*b(k,k2);%消去运算 end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素 b(k3,k)=0; %当前行第k列元素已消为0 end else %是第三行k=3 %= 第三行k=3的前代运算 = for k3=k1:N0 %从第四行到2n行（最后一行） for k2=k1:N1 %从第四列到2n+1列（即扩展列） b(k3,k2)=b(k3,k2)-b(k3,k)*b(k,k2);%消去运算（当前行3列元素消为0） end %用当前行K2列元素减去当前行3列元素乘以第三行K2列元素 b(k3,k)=0; %当前行第3列元素已消为0 end end endfor i=1:N0 DetaU(i,1)=b(i,N);endfunction bus,line = ReNum_(bus,line,nodenum)nb,mb=size(bus);nl,ml=size(line);bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus矩阵的数据k = 1;for i = 1 :nb for j = 1 : nb if bus(j,1) = k bus_temp(k,:) = bus(j,:); k = k + 1; end endend % 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号for i=1:nl for j=1:2 for k=1:nb if line(i,j)=nodenum(k,1) line(i,j)=nodenum(k,2); break end end end end % 恢复line的编号function S_res = S_res_(bus,line,YtYm)nl,ml=size(line);S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果S_res(:,1:2) = line(:,1:2); % 前两列为节点编号for k=1:nl I = S_res(k,1); J = S_res(k,2); if (J=0)&(I=0) S_res(k,3)=bus(I,2)*conj(bus(I,2)*(conj(YtYm(k,3)+conj(YtYm(k,4)-bus(I,2)*conj(bus(J,2)*conj(YtYm(k,3); S_res(k,4)=bus(J,2)*conj(bus(J,2)*(conj(YtYm(k,3)+conj(YtYm(k,5)-bus(J,2)*conj(bus(I,2)*conj(YtYm(k,3); S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流 elseif J=0 S_res(k,3)=bus(I,2)2*conj(YtYm(k,4); S_res(k,5)=S_res(k,3)+S_res(k,4); else S_res(k,4)=bus(J,2)2*conj(YtYm(k,5); S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流 endendfunction Y = y_(bus,line)nb,mb=size(bus);nl,ml=size(line);Y=zeros(nb,nb); % 对导纳矩阵赋初值0for k=1:nl I=line(k,1); J=line(k,2); Zt=line(k,3)+j*line(k,4); % 读入线路参数 if Zt=0 Yt=1/Zt; % 接地支路不计算Yt end Ym=line(k,5)+j*line(k,6); K=line(k,7); if (K=0)&(J=0) % 普通线路: K=0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); end if (K=0)&(J=0) % 对地支路: K=0,J=0,R=X=0 Y(I,I)=Y(I,I)+Ym; end if K0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/K/K; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J); end if K0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧 YtYm(k,3) = Yt/K; YtYm(k,4) = Ym+Yt*(K-1)/K; YtYm(k,5) = Yt*(1-K)/K/K; end if K0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧 YtYm(k,3) = -Yt*K; YtYm(k,4) = Ym+Yt*(1+K); YtYm(k,5) = Yt*(K2+K); endend六、运行结果- by longdinhohe - -牛顿拉夫逊法潮流计算结果- 节点计算结果：节点 节点电压 节点相角（角度） 节点注入功率 1 1.050000 + j 0.000000 0.000000 14.245276 + j 8.117797 2 0.756466 + j -0.427358-29.463862 -2.000000 + j -1.000000 3 0.124086 + j -0.130306-46.400557 -1.800000 + j -0.400000 4 0.214453 + j -0.775900-74.549595 -1.600000 + j -0.800000 5 1.027332 + j -1.902526-61.631703 -3.700000 + j -1.300000 6 1.084974 + j -1.932619-60.690108 5.000000 + j 8.117797线路计算结果：节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J) 2 1 -14.245276 + j -0.053001 14.245276 + j 8.117797 0.000000 + j 8.064796 2 3 8.865770 + j 2.935026 -1.933539 + j -0.046596 6.932231 + j 2.888429 3