dr69.m

这个m文件实现了Newton raphsom算法

clear
clc
tic;
basemva=10;
basekv=12.66;
% 1.bus no 2.realpower(KW) 3. Reactive Power(KVAR)
busdata=[1	0	0
2	0	0
3	0	0
4	0	0
5	0	0
6	2.6	2.2
7	40.4	30
8	75	54
9	30	22
10	28	19
11	145	104
12	145	104
13	8	5.5
14	8	5.5
15	0	0
16	45.5	30
17	60	35
18	60	35
19	0	0
20	1	0.6
21	114	81
22	5.3	3.5
23	0	0
24	28	20
25	0	0
26	14	10
27	14	10
28	26	18.6
29	26	18.6
30	0	0
31	0	0
32	0	0
33	14	10
34	19.5	14
35	6	4
36	26	18.55
37	26	18.55
38	0	0
39	24	17
40	24	17
41	1.2	1
42	0	0
43	6	4.3
44	0	0
45	39.2	26.3
46	39.2	26.3
47	0	0
48	79	56.4
49	384.7	274.5
50	384.7	274.5
51	40.5	28.3
52	3.6	2.7
53	4.3	3.5
54	26.4	19
55	24	17.2
56	0	0
57	0	0
58	0	0
59	100	72
60	0	0
61	1244	888
62	32	23
63	0	0
64	227	162
65	59	42
66	18	13
67	18	13
68	28	20
69	28	20];
% starting node end node Bus Reactive power(MW) resitance(ohms)  reactance(ohms)

linedata=[1	2	0.0005	0.0012	0	1
2	3	0.0005	0.0012	0	1
3	4	0.0015	0.0036	0	1
4	5	0.0251	0.0294	0	1
5	6	0.366	0.1864	0	1
6	7	0.3811	0.1941	0	1
7	8	0.0922	0.047	0	1
8	9	0.0493	0.0251	0	1
9	10	0.819	0.2707	0	1
10	11	0.1872	0.0619	0	1
11	12	0.7114	0.2351	0	1
12	13	1.03	0.34	0	1
13	14	1.044	0.345	0	1
14	15	1.058	0.3496	0	1
15	16	0.1966	0.065	0	1
16	17	0.3744	0.1238	0	1
17	18	0.0047	0.0016	0	1
18	19	0.3276	0.1083	0	1
19	20	0.2106	0.0696	0	1
20	21	0.3416	0.1129	0	1
21	22	0.014	0.0046	0	1
22	23	0.1591	0.0526	0	1
23	24	0.3463	0.1145	0	1
24	25	0.7488	0.2475	0	1
25	26	0.3089	0.1021	0	1
26	27	0.1732	0.0572	0	1
3	28	0.0044	0.0108	0	1
28	29	0.064	0.1565	0	1
29	30	0.3978	0.1315	0	1
30	31	0.0702	0.0232	0	1
31	32	0.351	0.116	0	1
32	33	0.839	0.2816	0	1
33	34	1.708	0.5646	0	1
34	35	1.474	0.4873	0	1
3	36	0.0044	0.0108	0	1
36	37	0.064	0.1565	0	1
37	38	0.1053	0.123	0	1
38	39	0.0304	0.0355	0	1
39	40	0.0018	0.0021	0	1
40	41	0.7283	0.8509	0	1
41	42	0.31	0.3623	0	1
42	43	0.041	0.0478	0	1
43	44	0.0092	0.0116	0	1
44	45	0.1089	0.1373	0	1
45	46	0.0009	0.0012	0	1
4	47	0.0034	0.0084	0	1
47	48	0.0851	0.2083	0	1
48	49	0.2898	0.7091	0	1
49	50	0.0822	0.2011	0	1
8	51	0.0928	0.0473	0	1
51	52	0.3319	0.1114	0	1
9	53	0.174	0.0886	0	1
53	54	0.203	0.1034	0	1
54	55	0.2842	0.1447	0	1
55	56	0.2813	0.1433	0	1
56	57	1.59	0.5337	0	1
57	58	0.7837	0.263	0	1
58	59	0.3042	0.1006	0	1
59	60	0.3861	0.1172	0	1
60	61	0.5075	0.2585	0	1
61	62	0.0974	0.0496	0	1
62	63	0.145	0.0738	0	1
63	64	0.7105	0.3619	0	1
64	65	1.041	0.5302	0	1
11	66	0.2012	0.0611	0	1
66	67	0.0047	0.0014	0	1
12	68	0.7394	0.2444	0	1
68	69	0.0047	0.0016	0	1
];
BZ=(basekv^2)/basemva;
linedata(:,3)=(1/BZ)*linedata(:,3);
linedata(:,4)=(1/BZ)*linedata(:,4);
P=-(.001/basemva)*busdata(:,2)';
Q=-(.001/basemva)*busdata(:,3)';
n=length(P);
% building admittance matrix
nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);
X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);
nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));
Z = R + j*X; y= ones(nbr,1)./Z;        %branch admittance
for n = 1:nbr
Ybus=zeros(nbus,nbus);     % initialize Ybus to zero
               % formation of the off diagonal elements
for k=1:nbr;
       Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);
       Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));
    end
end
              % formation of the diagonal elements
for  n=1:nbus
     for k=1:nbr
         if nl(k)==n
         Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);
         elseif nr(k)==n
         Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);
         else, end
     end
end
PQ=P+i*Q;
V=ones(n,1);
chV=zeros(n-1,1);
accuracy=1;
iter=0;
while accuracy >1e-7;
    iter=iter+1;
    V1=V(2:n);
    V=[1;V1]+conj([0; chV]);
    cPQ=diag(V)*conj(Ybus)*conj(V);
    chPQ=PQ.'-cPQ;
    chPQ1=chPQ(2:n);
    J=diag(V)*conj(Ybus)+diag(conj(Ybus*V));
    J1=J(2:n,2:n);
    chV=(inv(J1)*chPQ1);
    accuracy=max(abs(chV));
end
   toc;
    s=abs(V);
    Tloss=basemva*sum(cPQ)
    plot(s)