复杂网络N-R法潮流分析与计算的设计_课程设计.doc

电气工程及其自动化专业课程设计复杂网络n-r法潮流分析与计算的设计学生学号： 学生姓名： 班 级： 指导教师： 起止日期： 哈尔滨工程大学自动化学院课程设计报告撰写内容一、设计要求（宋体，小四号字，加黑）用matlab编程，n_r法计算潮流分布具体要求为：（1）给出程序，并给出注释（2）输出迭代次数，各节点电压，各支路电流（3）在图中标明功率流向节点数据如下表所示（标幺值）123456p31.80.63.55q10.50.81.3v11.05支路及变压器数据线路t1t2l2l3l4l5阻抗j0.04j0.020.06+j0.0250.01+j0.20.06+j0.50.05+j0.3导纳/2j0.25j0.25j0.25j0.25变比1.05:11.05:1精度要求：0.0001二、设计方案（要求给出详细的设计思路及其必要的论证）（1.）潮流计算的方法（1）高斯雅克比迭代法（2）高斯-塞得尔法（对初值要求底，迭代次数多）（3）牛顿-拉夫逊法（使用广泛）（4）pq 快速分解法（提升运算速度）目前广泛应用的潮流计算方法都是基于节点电压法的，以节点导纳矩阵y作为电力网络的数学模型。节点电压ui 和节点注入电流ii 由节点电压方程yv=i （1）根据s=vi(i为i 的共轭)可得非线性的节点方程yv=i=(s/v) （2）在实际的电力系统中，已知的运行条件不是节点的注入电流，而是负荷和发电机的功率，而且这些功率一般不随节点电压的变化而变化。由于各节点注入功率与注入电流的关系为sipijqi =viii，因此可将式（2）改写为ii=si/vi=pi+jqi/vi (i= 1, 2,3 n) （3）式中，pi 和qi 分别为节点i 向网络注入的有功功率和无功功率，当i 为发电机节点时pi0；当i 为负荷节点时pi0；当i 为无源节点pi0，qi0；vi 和ii 分别为节点电压相量vi 和节点注入电流相量ii 的共轭。式（3）亦即潮流计算的基本方程式,它可以在直角坐标也可以在极坐标上建立2n 个实数形式功率方程式。发电机pi、qi 为正，负荷pi、qi 为负。展开yv=i 为ii=yijvj=yiivi+yijvi( i=1 2 3 n) (4)将式（4）代入式（3），得n 维的非线性复数的电压方程组潮流计算的基本方程为 (pi-jqi)/ vi= yiivi+yijvi （i=1, 2, n） (5)（2.）变量的分类假设系统中有n 个节点，构成n 个复数方程，2n 个实数方程，变量总数为6n 个。a)不可控变量（2n 个）：负荷消耗的有功功率li p 和无功功率li q .由于该类变量无法控制，取决于用户，而且出现事先没有预计的变动，使系统偏离原始运行状态，因此又称为不可控变量或扰动变量。b)控制变量（2n 个）:发电机发出的有功功率gi p 和无功功率gi q ，因为该类变量可控。也称独立变量。c)状态变量（2n 个）：母线电压或节点电压的幅值大小i v 与相角大小i ，又称依从变量或因变量。并且i v 受gi p 控制， i 受gi q 控制。其中2n 个扰动变量是给定的，2n 个控制变量和2n 个状态变量中给定两个，求另外两个。（3.）变量的约束条件a)扰动变量没有约束条件。b)控制变量约束条件：为满足发电机的技术经济特性指标。c)状态变量的i v 的约束条件：保证良好的电能质量。d) 状态变量的i 的约束条件：保证系统的稳定运行。(4.)系统节点的分类，根据给定的控制变量和状态变量进行分类如下：（1）pq 节点（即负荷节点）：给定pi、qi，求vi 和i （ i i e , f ）。通常变电所都是这一类型的节点，由于没有发电设备，因而发电功率为零电力系统中的绝大多数节点属于这一节点。其包含变电站节点（即联络节点或浮游节点）。（2）pv 节点（即调节节点、电压控制节点）：给定pi 和vi，求qi 和i （ i i e , f ）。这类节点必须有足够的可调无功容量，用以维持给定的电压幅值。一般时选择有一顶武功储备的发电厂和具有可调无功电源设备的变电所作为pv 节点。在电力系统中，这类节点数很少。（3）平衡节点（即松弛节点、参考节点、基准节点）：给定vi 和i （ i =0），求pi 和qi。（只有一个）有功功率不能给定，这个节点承担了系统的有功功率平衡。同时其电压幅值也是给定的，相位为零。(5. )p-q 分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时，修正方程式展开为：p = h + nv/ vq = j + lv /v以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特点。电力系统中有功功率主要与各节点电压向量的角度有关，无功功率则主要受各节点电压幅值的影响。大量运算经验也告诉我们，矩阵n 及j 中各元素的数值相对是很小的，因此对牛顿法的第一步简化就是把有功功率和无功功率分开来进行迭代，即将式(4)化简为：p = hq = lv /v (5)这样，由于我们把2n 阶的线性方程组变成了二个n 阶的线性方程组，对牛顿法的第二个化简，也是比较关键的一个化简，即把式(5)中的系数矩阵简化为在迭代过程中不变的对称矩阵。众所周知，一般线路两端电压的相角差是不大的(通常不超过1020 度)，因此可以认为： (6)此外，与系统各节点无功功率相应的导纳li b 必定远远小于该节点自导纳的虚部，即：因此， (7)考虑到以上关系后，式(5)中系数矩阵中的元素表达式可以化简为： （8）这样，式(5)中系数矩阵可以表示为： （9）进一步可以把它们表示为以下矩阵的乘积： （10）将它代入(5)中，并利用乘法结合率，我们可以把修正方程式变为：将以上两式的左右两侧用以下矩阵左乘就可得到以上两式就是p-q 分解法达到修正方程式，其中系数矩阵只不过是系统导纳矩阵的虚部，因而是对称矩阵，而且在迭代过程中维持不变。它们与功率误差方程式构成了p-q 分解法迭代过程中基本计算公式，其迭代步骤大致是：(1)给定各节点电压向量的电压初值v i (0) ,i (0);(2)根据(12)计算各节点有功功率误差pi ，并求出;pi/vi(3)解修正方程式(11)，并进而计算各节点电压向量角度的修正量i (4)修正各节点电压向量角度i ；(5)根据式(16)计算各节点无功功率误差i q ，并求出 / ; i i q v(6)解修正方程式(11)，求出各节点电压幅值的修正量i v(7)修正各节点电压幅值i v (i ) (i 1) (i 1)i i i v =v v (18)(8)返回(2)进行迭代，直到各节点功率误差及电压误差都满足收敛条件。p-q 分解法与牛顿法潮流程序的主要差别表现在它们的修正方程式上。p-q分解法通过对电力系统具体特点的分析，对牛顿法修正方程式的雅可比矩阵进行了有效的简化和改进，有以下三个特点：（1）在提高计算速度和减少内存方面的作用是明显的，不再叙述。（2）使我们得到以下好处。首先，因为修正方程式的系数矩阵就是导纳矩阵的虚部，因此在迭代过程中不必象牛顿法那样进行形成雅可比矩阵的计算，这样不仅是仅减少了运算量，而且也大大简化了程序。其次，由于系数矩阵在迭代过程中维持不变，因此在求解修正方程式时，可以迅速求得修正量，从而显著提高了迭代速度。（3）可以使我们减少形成因子表时的运算量，而且由于对称矩阵三角分解后，其上三角矩阵和下三角矩阵有非常简单的关系，所以在计算机中可以只储存上三角矩阵或下三角矩阵，从而也进一步节约了内存。三、设计内容%本程序的功能是用牛顿拉夫逊法进行潮流计算% b1矩阵：1、支路首端号； 2、末端号； 3、支路阻抗； 4、线路对地电纳 (或变压器导纳)；% 5、支路的变比； 6、支路首端处于k侧为1，1侧为0；% 7、线路/变压器标识（0/1）变压器参数当支路首端处于k侧标识为1时归算至末端侧，0归算至首端侧% b2矩阵：1、该节点发电机功率； 2、该节点负荷功率；% 3、pq节点电压初始值，或平衡节点及pv节点电压的给定值% 4、节点所接无功补偿并联电容（感）的电纳% 5、节点分类标号：1为平衡节点（应为1号节点）；2为pq节点；3为pv节点；clear;isb=1; %input(请输入平衡母线节点号：isb=);pr=0.0001; %input(请输入误差精度：pr=);%-n=6;%input(请输入节点数：n=);nl=6;%input(请输入支路数：nl=);b1=1 2 0+0.04i 0 1.05 1 1; 2 3 0.06+0.025i 0+0.5i 1 0 0; 2 5 0.01+0.2i 0+0.5i 1 0 0; 3 4 0.06+0.50i 0+0.5i 1 0 0; 4 5 0.05+0.3i 0 1 0 0; 6 5 0+0.02i 0 1.05 1 1b2=0 0 1 0 1; 0 3+1i 1.00 0 2; 0 1.8+0.50i 1.00 0 2; 0 0.6+0.8i 1.00 0 2; 0 3.5+1.3i 1.00 0 2; 0 -5+0i 1.05 0 3 %input(请输入各节点参数形成的矩阵： b2=);%x=1 0;2 0;3 0;4 0;5 0;6 0%-%n=4;%input(请输入节点数：n=);nl=4;%input(请输入支路数：nl=);%b1=1 2 4+16i 0 1 0 0;1 3 4+16i 0 1 0 0;2 3 2+8i 0 1 0 0;2 4 1.49+48.02i 0 11/110 0 1 %input(请输入由支路参数形成的矩阵： b1=);%b2=0 0 115 0 1;0 0 110 0 2;0 20+4i 110 0 2;0 10+6i 10 0 2 %input(请输入各节点参数形成的矩阵： b2=);%-y=zeros(n);e=zeros(1,n);f=zeros(1,n);v=zeros(1,n);sida=zeros(1,n);s1=zeros(nl);% % %-求导纳矩阵-%for i=1:n % if x(i,2)=0; % p=x(i,1); % y(p,p)=1/x(i,2); %end%endfor i=1:nl%从1到n1（总支路数）if b1(i,7)=1 %-如果是变压器支路- if b1(i,6)=0%左节点(首端)处于1侧 p=b1(i,1);q=b1(i,2); else %左节点(首端)处于k侧 p=b1(i,2);q=b1(i,1); end y(p,q)=y(p,q)-1./(b1(i,3)*b1(i,5);%非对角元 y(q,p)=y(p,q); %非对角元 y(q,q)=y(q,q)+1./(b1(i,3)*b1(i,5)2);%对角元k侧 y(p,p)=y(p,p)+1./b1(i,3)+b1(i,4);%对角元1侧+励磁导纳 else %-否则为线路支路-p=b1(i,1);q=b1(i,2);y(p,q)=y(p,q)-1./b1(i,3);%非对角元 y(q,p)=y(p,q); %非对角元 y(q,q)=y(q,q)+1./b1(i,3)+b1(i,4)./2.0000;%对角元j侧+线路电纳的一半 y(p,p)=y(p,p)+1./b1(i,3)+b1(i,4)./2.0000;%对角元i侧+线路电纳的一半endenddisp(导纳矩阵 y=);disp(y);%-给定各节点初始电压及给定各节点注入功率-g=real(y);b=imag(y);%分解出导纳阵的实部和虚部 for i=1:n%给定各节点初始电压的实部和虚部 e(i)=real(b2(i,3); f(i)=imag(b2(i,3); v(i)=abs(b2(i,3);%pv、平衡节点及pq节点电压模值 endfor i=1:n%给定各节点注入功率 s(i)=b2(i,1)-b2(i,2); %i节点注入功率sg-sl b(i,i)=b(i,i)+b2(i,4);%i节点无功补偿量(电纳值)end%=用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）=p=real(s);q=imag(s); %分解出各节点注入的有功和无功功率ict1=0;it2=1;n0=2*n;n1=n0+1;a=0; %迭代次数ict1、a；不满足收敛要求的节点数it2while it2=0 % n0=2*n 雅可比矩阵的阶数；n1=n0+1扩展列 it2=0;a=a+1; jz=jacobi矩阵第(,num2str(a),)次消去运算;jz1=jacobi矩阵第(,num2str(a),)次回代运算;jz0=功率方程第(,num2str(a),)次差值：; %-求取各个节点的功率及功率偏差及pv节点的电压偏差- for i=1:n %n个节点2n行(每节点两个方程p和q或u) p=2*i-1;m=p+1;c(i)=0;d(i)=0; for j1=1:n %第i行共n列(n个节点间互导纳及节点电压相乘即电流) c(i)=c(i)+g(i,j1)*e(j1)-b(i,j1)*f(j1);%(gij*ej-bij*fj) d(i)=d(i)+g(i,j1)*f(j1)+b(i,j1)*e(j1);%(gij*fj+bij*ej) end %求i节点有功和无功功率p,q的计算值 p1=c(i)*e(i)+f(i)*d(i);%节点功率p计算ei(gij*ej-bij*fj)+fi(gij*fj+bij*ej) q1=c(i)*f(i)-e(i)*d(i);%节点功率q计算fi(gij*ej-bij*fj)-ei(gij*fj+bij*ej) v2=e(i)2+f(i)2;%电压模平方 %=求取功率差及pv节点电压模平方差 = if i=isb%非平衡节点(pq或pv节点) if b2(i,5)=3%非pv节点(只能是pq节点) j(m,n1)=p(i)-p1;%pq节点有功功率差j(m,n1)扩展列p j(p,n1)=q(i)-q1; %pq节点无功功率差j(p,n1)扩展列q else %pv节点= j(m,n1)=p(i)-p1;%pv节点有功功率差j(m,n1)扩展列p j(p,n1)=v(i)2-v2;%pv节点电压模平方差j(p,n1)扩展列u end end %(if i=isb)非平衡节点(pq或pv节点) end %(for i=1:n) n个节点2n行(每节点两个方程p和q或u) for m=1:n0 jjn1(m)=j(m,n1); end disp(jz0);disp(jjn1); %-判断功率偏差量及pv节点的电压偏差量是否满足要求- for k=3:n0 %除去平衡节点1、2号以外的所有节点 det=abs(j(k,n1); if det=pr; %pq节点的功率偏差量及pv节点的电压偏差量是否满足要求 it2=it2+1; %不满足要求的节点数加1 end end ict2(a)=it2; %不满足要求的节点数；a为迭代次数 ict1=ict1+1; %迭代次数 if ict2(a)=0; %当前不满足要求的节点数为零 break %退出迭代运算 end %-以上为求取各个节点的功率及功率偏差及pv节点的电压偏差- %= 求取jacobi矩阵形成修正方程 = for i=2:n %n个节点2n行(每节点两个方程p和q或u) if i=isb%非平衡节点(pq或pv节点) if b2(i,5)=3 %下面是针对pq节点来求取jacobi矩阵的元素 = c(i)=0;d(i)=0; for j1=1:n %第i行共n列(n个节点间互导纳及节点电压相乘即电流) c(i)=c(i)+g(i,j1)*e(j1)-b(i,j1)*f(j1);%(gij*ej-bij*fj) d(i)=d(i)+g(i,j1)*f(j1)+b(i,j1)*e(j1);%(gij*fj+bij*ej) end for j1=2:n %第i行共n列(2n个jacobi矩阵元素dp/de及dp/df或dq/de及dq/df) if j1=isb&j1=i%非平衡节点&非对角元 x1=-g(i,j1)*e(i)-b(i,j1)*f(i);% x1=dp/de=-dq/df=-x4 x2=b(i,j1)*e(i)-g(i,j1)*f(i);% x2=dp/df=dq/de=x3 x3=x2; % x2=dp/df x3=dq/de x4=-x1; % x1=dp/de x4=dq/df p=2*i-1;q=2*j1-1; j(p,q)=x3;m=p+1; % x3=dq/de j(p,n)=dq节点无功功率差 j(p,n)=dq; j(m,q)=x1;q=q+1; % x1=dp/de j(m,n)=dp节点有功功率差 j(m,n)=dp; j(p,q)=x4;j(m,q)=x2; % x4=dq/df x2=dp/df elseif j1=i&j1=isb%非平衡节点&对角元 x1=-c(i)-g(i,i)*e(i)-b(i,i)*f(i);% dp/de x2=-d(i)+b(i,i)*e(i)-g(i,i)*f(i);% dp/df x3=d(i)+b(i,i)*e(i)-g(i,i)*f(i); % dq/de x4=-c(i)+g(i,i)*e(i)+b(i,i)*f(i);% dq/df p=2*i-1;q=2*j1-1;j(p,q)=x3;%扩展列q j(p,n)=dq; m=p+1; j(m,q)=x1;q=q+1;j(p,q)=x4;%扩展列p j(m,n)=dp; j(m,q)=x2; end end else %if b2(i,5)=3 % 否则（即为pv节点） %= 下面是针对pv节点来求取jacobi矩阵的元素 = for j1=1:n if j1=isb&j1=i%非平衡节点&非对角元 x1=-g(i,j1)*e(i)-b(i,j1)*f(i); % dp/de x2=b(i,j1)*e(i)-g(i,j1)*f(i); % dp/df x5=0;x6=0; p=2*i-1;q=2*j1-1;j(p,q)=x5; % pv节点电压误差j(p,n)=dv; m=p+1; j(m,q)=x1;q=q+1;j(p,q)=x6; % pv节点有功误差j(m,n)=dp; j(m,q)=x2; elseif j1=i&j1=isb %非平衡节点&对角元 x1=-c(i)-g(i,i)*e(i)-b(i,i)*f(i);% dp/de x2=-d(i)+b(i,i)*e(i)-g(i,i)*f(i);% dp/df x5=-2*e(i); x6=-2*f(i); p=2*i-1;q=2*j1-1;j(p,q)=x5; % pv节点电压误差j(p,n)=dv; m=p+1; j(m,q)=x1;q=q+1;j(p,q)=x6; % pv节点有功误差j(m,n)=dp; j(m,q)=x2; end end end %(if b2(i,5)=3 else) end %(if i=isb) end %(for i=1:n)n个节点2n行(每节点两个方程p和q或u) jz0=形成的第(,num2str(a),)次jacobi矩阵：; disp(jz0);disp(j); %= 以上为形成完整的jacobi矩阵 = %=下面用高斯消去法对由jacobi矩阵形成的修正方程进行求解(按列消去、回代) = for k=3:n0 % n0=2*n （从第三行开始，第一、二行是平衡节点） for k1=k+1:n1% 从k+1列的jacobi元素到扩展列的p、q 或 u j(k,k1)=j(k,k1)./j(k,k);% 用k行k列对角元素去除k行k列后的非对角元素进行规格化 end j(k,k)=1; % 对角元规格化k行k列对角元素赋1 %= 按列消去运算 = for k2=k+1:n0 % 从k+1行到2*n最后一行 for k3=k+1:n1 %从k2+1列到扩展列消去k+1行后各行下三角元素 j(k2,k3)=j(k2,k3)-j(k2,k)*j(k,k3);%消去运算 end %用当前行k3列元素减去当前行k列元素乘以第k行k3列元素 j(k2,k)=0; %当前行第k列元素已消为0 end end jz=jacobi矩阵第(,num2str(a),)次消去运算;jz1=jacobi矩阵第(,num2str(a),)次回代运算; disp(jz);disp(j); %= 按列回代运算 = for k=n0:-1:3 for k1=k-1:-1:3 j(k1,n1)=j(k1,n1)-j(k1,k)*j(k,n1); j(k1,k)=0; end end for m=1:n0 jjn1(m)=j(m,n1); end disp(jz1);disp(jjn1);%disp(j); %-修改节点电压- for k=3:2:n0-1 l=(k+1)./2; e(l)=e(l)-j(k,n1); %修改节点电压实部 k1=k+1; f(l)=f(l)-j(k1,n1); %修改节点电压虚部 u(l)=sqrt(e(l)2+f(l)2); end disp(各个节点电压模);disp(u); %= 结束一次迭代 =end%* 下面为迭代计算结束后的有关输出过程 *disp(迭代次数：);disp(ict1-1);disp(没有达到精度要求的个数：);disp(ict2);for k=1:n v(k)=sqrt(e(k)2+f(k)2); %计算各节点电压的模值 sida(k)=atan(f(k)./e(k)*180./pi; %计算各节点电压的角度 e(k)=e(k)+f(k)*j; %将各节点电压用复数表示end%= 计算各输出量 =disp(各节点的电压复数值e为(节点号从小到大排列)：);disp(e); %显示各节点的实际电压值e用复数表示disp(-);disp(各节点的电压模值大小v为(节点号从小到大排列)：);disp(v); %显示各节点的电压大小v的模值disp(-);disp(各节点的电压相角sida为(节点号从小到大排列)：);disp(sida); %显示各节点的电压相角for p=1:n c(p)=0; for q=1:n c(p)=c(p)+conj(y(p,q)*conj(e(q); %计算各节点注入电流的共轭值 end s(p)=e(p)*c(p); %计算各节点的功率 s = 电压 x 注入电流的共轭值enddisp(各节点的功率s为(节点号从小到大排列)：);disp(s); %显示各节点的注入功率disp(-);disp(各条支路的首端功率si为(顺序同您输入b1时一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); if b1(i,7)=0 si(p,q)=e(p)*conj(e(p)*b1(i,4)./2+(e(p)-e(q)./b1(i,3); siz(i)=si(p,q); else if b1(i,6)=0 si(p,q)=e(p)*(conj(e(p)*b1(i,4). +(e(p)*b1(i,5)-e(q)*(1./(b1(i,3)*b1(i,5); siz(i)=si(p,q); else si(p,q)=e(p)*conj(e(p)-e(q)*b1(i,5)*(1./(b1(i,3)*b1(i,5)2); siz(i)=si(p,q); end end zf=s(,num2str(p),num2str(q),)=,num2str(si(p,q); disp(zf); disp(-);enddisp(各条支路的末端功率sj为(顺序同您输入b1时一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); if b1(i,7)=0 sj(q,p)=e(q)*conj(e(q)*b1(i,4)./2+(e(q)-e(p)./b1(i,3); sjy(i)=sj(q,p); else if b1(i,6)=0 sj(q,p)=e(q)*conj(e(q)-e(p)*b1(i,5)*(1./(b1(i,3)*b1(i,5)2); sjy(i)=sj(q,p); else sj(q,p)=e(q)*(conj(e(q)*b1(i,4). +(e(q)*b1(i,5)-e(p)*(1./(b1(i,3)*b1(i,5); sjy(i)=sj(q,p); end end zf=s(,num2str(q),num2str(p),)=,num2str(sj(q,p); disp(zf); disp(-);enddisp(各条支路的功率损耗ds为(顺序同您输入b1时一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); ds(i)=si(p,q)+sj(q,p); zf=ds(,num2str(p),num2str(q),)=,num2str(ds(i); disp(zf); disp(-); end zws=0;jddy=;jdp=;jdq=;jddyjd=; for i=1:n %总网损为所有节点注入功率的代数和 zws=zws+s(i); jddyjd=strcat(jddyjd,num2str(i),(,num2str(sida(i),),); jddy=strcat(jddy,num2str(i),(,num2str(v(i),),); jdp=strcat(jdp,num2str(i),(,num2str(real(s(i),),); jdq=strcat(jdq,num2str(i),(,num2str(imag(s(i),),);endjddyjd=strcat(节点电压角度：,jddyjd);jddy=strcat(节点电压模值：,jddy);jdp=strcat(节点有功：,jdp);jdq=strcat(节点无功：,jdq); zf subplot(4,1,2);plot(sida);title(jddyjd);ylabel(节点电压角度);grid on;subplot(4,1,3);p=real(s);q=imag(s);bar(p);title(jdp);ylabel(节点注入有功);grid on;subplot(4,1,4);bar(q);title(jdq);xlabel(zf);ylabel(节点注入无功);grid on;%*figure(2);subplot(3,2,1);jdh=;jdh1=;for i=1:nl jdh=strcat(jdh,num2str(i),(,num2str(b1(i,1),num2str(b1(i,2),), ); jdh1=strcat(jdh1,num2str(i),(,num2str(b1(i,2),num2str(b1(i,1),), );endp1=real(siz);q1=imag(siz);bar(p1);title(jdh);ylabel(支路首端注入有功);%xlabel(支路号);grid on;subplot(3,2,2);bar(q1);title(jdh);ylabel(支路首端注入无功);grid on;subplot(3,2,3);p2=real(sjy);q2=imag(sjy);bar(p2);