[Débutant] Liaison function MATLAB

salut , j'ai 6 fichiers matlab function (.m) ,et ils devraient marcher ensemble ,care chaque fonction est associée a l'autre , alors je veux bien si quelqu’ un pourrai m'indiquer comment faire la liaison entre ces programmes , voci les prog :
rq: les fonctions sont correctes , ce qui manque c'est la liaison .

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
function [baseMVA,bus,gen,branch]=data_name
baseMVA = 100.0;
%bus type     Pd    Qd           Gs         Bs     area        Vm             Va         baseKV zone      Vmax       Vmin
bus = [
1    2       0.00   0.00       0.000       0.000    1        1.02000         0.000       500.00   0     1.06000     0.94000
2    2       0.00   0.00       0.000       0.000    1        1.02000         0.000       500.00   0     1.06000     0.94000
3    1       0.00   0.00       0.000       -600.000 1        1.00000         0.000       230.00   0     1.06000     0.94000
4    1       0.00   0.00       0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
5    1       200.00 100.00     0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
6    1       1000.00 800.00    0.000       200.000  1        1.00000         0.000       21.40    0     1.06000     0.94000
7    2       0.00   0.00       0.000       -300.000 1        1.05000         0.000       230.00   0     1.06000     0.94000
8    1       0.00   0.00       0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
9    1       300.00 150.00     0.000       50.000   1        1.00000         0.000       69.00    0     1.06000     0.94000
10   1       0.00   0.00       0.000       0.000    1        1.00000         0.000       69.00    0     1.06000     0.94000
11   1       1200.00 700.00    0.000       300.000  1        1.00000         0.000       69.00    0     1.06000     0.94000
12   2       0.00   0.00       0.000       0.000    1        0.98000         0.000       69.00    0     1.06000     0.94000
13   1       0.00   0.00       0.000       0.000    1        1.00000         0.000       69.00    0     1.06000     0.94000
14   1       0.00   0.00       0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
15   1       0.00   0.00       0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
16   1       100.00  50.00     0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
17   1       0.00   0.00       0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
18   1       200.00  75.00     0.000       0.000    1        1.00000         0.000       230.00   0     1.06000     0.94000
19   2       200.00  75.00     0.000       0.000    1        1.02000         0.000       230.00   0     1.06000     0.94000
20   3       0.00    0.00      0.000       0.000    1        1.05000         0.000       230.00   0     1.06000     0.94000
];
%bus Pg Qg Qmax Qmin Vsp baseMVA status Pmax Pmin
gen = [
1 750.00 0.00 600.00 -100.00 1.02000 100.00 1 1000.00 0.00
2 750.00 0.00 600.00 -100.00 1.02000 100.00 1 1000.00 0.00
7 600.00 0.00 400.00 -100.00 1.05000 100.00 1 1000.00 0.00
12 800.00 0.00 600.00 0.00 0.98000 100.00 1 1000.00 0.00
19 100.00 0.00 80.00 0.00 1.02000 100.00 1 1000.00 0.00
20 0.00 0.00 600.00 -100.00 1.05000 100.00 1 1000.00 0.00
];
%fbus tbus r x b ratea rateb ratec ratio angle status
branch = [
3 1 0.00000 0.01000 0.00000 9900.00 0.00 0.00 0.00000 0.000 1
3 2 0.00000 0.01000 0.00000 9900.00 0.00 0.00 0.00000 0.000 1
3 4 0.00260 0.04600 3.50000 9900.00 0.00 0.00 0.00000 0.000 1
3 4 0.00260 0.04600 3.50000 9900.00 0.00 0.00 0.00000 0.000 1
3 7 0.00100 0.01500 1.20000 9900.00 0.00 0.00 0.00000 0.000 1
4 5 0.00000 0.00500 0.00000 9900.00 0.00 0.00 1.01000 0.000 1
4 8 0.00080 0.01000 0.95000 9900.00 0.00 0.00 0.00000 0.000 1
4 15 0.00300 0.03000 2.50000 9900.00 0.00 0.00 0.00000 0.000 1
5 6 0.00500 0.04500 0.10000 9900.00 0.00 0.00 0.00000 0.000 1
5 6 0.00600 0.05400 0.15000 9900.00 0.00 0.00 0.00000 0.000 1
5 17 0.00100 0.01200 0.03000 9900.00 0.00 0.00 0.00000 0.000 1
6 9 0.00400 0.04000 0.10000 9900.00 0.00 0.00 0.00000 0.000 1
6 11 0.00033 0.00333 0.09000 9900.00 0.00 0.00 0.00000 0.000 1
6 19 0.00270 0.02200 0.30000 9900.00 0.00 0.00 0.00000 0.000 1
7 8 0.00200 0.02500 2.00000 9900.00 0.00 0.00 0.00000 0.000 1
7 10 0.00300 0.03000 2.50000 9900.00 0.00 0.00 0.00000 0.000 1
8 9 0.00000 0.01000 0.00000 9900.00 0.00 0.00 1.01000 0.000 1
9 11 0.00500 0.04500 0.08000 9900.00 0.00 0.00 0.00000 0.000 1
9 11 0.00500 0.04500 0.08000 9900.00 0.00 0.00 0.00000 0.000 1
10 11 0.00000 0.01000 0.00000 9900.00 0.00 0.00 1.01000 0.000 1
11 12 0.00000 0.01000 0.00000 9900.00 0.00 0.00 1.02000 0.000 1
20 13 0.00000 0.02500 0.00000 9900.00 0.00 0.00 0.00000 0.000 1
20 14 0.00000 0.00800 0.00000 9900.00 0.00 0.00 0.00000 0.000 1
13 15 0.00600 0.05400 0.09000 9900.00 0.00 0.00 0.00000 0.000 1
14 16 0.00600 0.05400 0.09000 9900.00 0.00 0.00 0.00000 0.000 1
14 16 0.00600 0.05400 0.09000 9900.00 0.00 0.00 0.00000 0.000 1
15 16 0.00000 0.01000 0.00000 9900.00 0.00 0.00 0.00000 0.000 1
16 17 0.00350 0.03000 0.07000 9900.00 0.00 0.00 0.00000 0.000 1
16 18 0.00300 0.02500 0.06000 9900.00 0.00 0.00 0.00000 0.000 1
16 19 0.00600 0.05000 0.12000 9900.00 0.00 0.00 0.00000 0.000 1
18 19 0.00300 0.02500 0.06000 9900.00 0.00 0.00 0.00000 0.000 1
];

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
function pflow(data_name,tol,a,b)
global nb nbr
[baseMVA,bus,gen,branch] = feval(data_name);
[bus,gen,branch,ext]=intnum(bus,gen,branch); 
[Y]=network_f(branch,bus,baseMVA);
% If number of inputs is less than 1, take tolerance tol = 0.001.
% Otherwise, tol is given as an input.
if nargin < 2
tol = 0.001;
end
[Vm,deltad,Pgt,Qgt,Plt,Qlt,iter,PSt,QSt,PSr,QSr,PSl,QSl,PStt,QStt,JC,z,tol]= newtraph(bus,gen,Y,baseMVA,branch,tol);
[bus,gen,branch]=extnum(bus,gen,branch,ext); 
%nb = size(bus,1);
%nbr = size(branch,1);
ngn = size(gen,1);
bus(:,8)=Vm;
bus(:,9)=deltad;
bustype=bus(:,2);
fd= fopen('output.txt', 'wt');
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n| Solution Summary|');
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n Base MVA: %4d Tolerance: %4d',baseMVA,tol);
fprintf(fd, '\n');
fprintf(fd, '\n Power-flow solution converged in %2f seconds (in %1diterations).',z,iter);
fprintf(fd, '\n');
% Bus data starts here.
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n| Bus Data|');
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n Bus Bus Voltage GenerationLoad ');
fprintf(fd, '\n No Code Mag(pu) Ang(deg) P(MW) Q(MVAR) P(MW)Q(MVAR) ');
fprintf(fd, '\n --- ---- ------- -------- -------- -------- -------- -------- ');
for i = 1:nb
fprintf(fd, '\n%4d%6d%11.3f %8.2f', bus(i, [1,2, 8, 9]));
fprintf(fd, ' %10.2f%10.2f', Pgt(i),Qgt(i));
if bus(i, 3) | bus(i, 4)
fprintf(fd, '%10.2f%10.2f', bus(i, [3, 4]));
else
fprintf(fd, ' - - ');
end
end
fprintf(fd, '\n -------- -------- -------- -------- ');
fprintf(fd,' \n'); fprintf(fd,' Total:');
fprintf(fd,' %10.2f',sum( Pgt)); fprintf(fd,' %9.2f',sum( Qgt));
fprintf(fd,' %9.2f', Plt); fprintf(fd,' %9.2f', Qlt);
fprintf(fd, '\n');
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n| Branch Data|');
fprintf(fd,'\n=======================================================================');
fprintf(fd, '\n From To From Bus Injection To Bus InjectionLoss ');
fprintf(fd, '\n Bus Bus P(MW) Q(MVAR) P(MW) Q(MVAR) P(MW)Q(MVAR) ');
fprintf(fd, '\n --- --- -------- -------- -------- -------- --------------- ');
for i=1:nbr
fprintf(fd, '\n%4d%6d%', branch(i,[1,2]));
fprintf(fd, '%10.2f%10.2f', PSr(i)*baseMVA,QSr(i)*baseMVA);
fprintf(fd, '%10.2f%10.2f', PSl(i)*baseMVA,QSl(i)*baseMVA);
fprintf(fd, '%10.2f%10.2f', PSt(i)*baseMVA,QSt(i)*baseMVA);
end
fprintf(fd, '\n ------- --------');
fprintf(fd,' \n');, fprintf(fd,'Total:');
fprintf(fd,' %9.2f',sum( PStt)*baseMVA); fprintf(fd,' %9.2f',sum(QStt)*baseMVA);
fclose(fd);
if nargin == 3
len=length(Y);
I=fopen('yb.txt','w');
for i=1:len
for k=1:len
if Y(i,k) ~=0
fprintf(I,'\n(%3d,%3d) %10.2f%10.2f j',i,k,real(Y(i,k)),imag(Y(i,k)));
end
end
end
fclose(I);
end
if nargin == 4
Jm=length(JC);
I=fopen('Jac.txt','w');
for i=1:Jm
for k=1:Jm
if JC(i,k) ~=0
fprintf(I,'\n(%3d,%3d) %10.2f%10.2f ',i,k, JC(i,k));
end
end
end
fclose(I);

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function [Y,fbus,tbus,R,X,b,Bs,fbuslen]=network_f(branch,bus,baseMV)
 
j=sqrt(-1); i = sqrt(-1);
fbus=branch(:,1);
tbus=branch(:,2);
R=branch(:,3);
X=branch(:,4);
b=branch(:,5)/2;
Bs=bus(:,6)/baseMVA;
fbuslen=length(branch(:,1));
maxbus=max(max(fbus),max(tbus));
Z=R+j*X;
y=ones(fbuslen,1)./Z;
Y=zeros(maxbus,maxbus);
a=branch(:,9);
for n=1:maxbus
for k=1:fbuslen
if fbus(k)==n
if a(k)==0
Y(n,n) = Y(n,n)+y(k)+j*b(k);
else
Y(n,n) = Y(n,n)+y(k)/(a(k)^2) + j*b(k);
end
elseif tbus(k)==n
Y(n,n) = Y(n,n)+y(k)+j*b(k);
end
end
end
for n=1:maxbus
Y(n,n) = Y(n,n)+j*Bs(n);
end
for i=1:fbuslen
if a(i)==0
Y(fbus(i),tbus(i))=Y(fbus(i),tbus(i))-y(i);
Y(tbus(i),fbus(i))=Y(fbus(i),tbus(i));
else
Y(fbus(i),tbus(i))=Y(fbus(i),tbus(i))-y(i)/a(i);
Y(tbus(i),fbus(i))=Y(fbus(i),tbus(i));
end
end
len=length(Y);
I=fopen('yb.txt','w');
for i=1:len
for k=1:len
if Y(i,k) ~=0
fprintf(I,'\n(%3d,%3d) %10.2f%10.2f j',i,k,real(Y(i,k)),imag(Y(i,k)));
end
end
end
fclose(I);

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
function [Vm,deltad,Pgt,Qgt,Plt,Qlt,iter,PSt,QSt,PSr,QSr,PSl,QSl,PStt,QStt,JC,z,tol] = newtraph(bus,gen,Y,baseMVA,branch,tol);
% NEWTRAPH Calculates power flow using Newton-Raphson method.
 
global nb nbr tic;
nb = length(bus(:,1));
nbr = size(branch,1);
fbus = branch(:,1);
tbus = branch(:,2);
% Reading the system data and initial values from the file-------
bustype = bus(:,2);
V = bus(:,8);
delta = bus(:,9);
PL = bus(:,3)/baseMVA;
QL = bus(:,4)/baseMVA;
QBs = bus(:,6);
PG = zeros(nb,1);
Qmax = zeros(nb,1);
Qmin = zeros(nb,1);
gindx = find(bustype == 2);
ln = length(gindx);
for i = 1:ln
ing(i) = find(gen(:,1) == gindx(i));
end
PG(gindx) = gen(ing,2)/baseMVA;
Qmax(gindx) = gen(ing,4)/baseMVA;
Qmin(gindx) = gen(ing,5)/baseMVA;
j=sqrt(-1);
for n=1:nb
if V(n)<=0 V(n)=1.0; V(n)=1+j*0;
else delta(n)=pi/180*delta(n);
Vm(n)=abs(V(n)*(cos(delta(n))+j*sin(delta(n))));
end
end
% Computing the magnitudes and angles of the Ybus matrix elements--
Ym=abs(Y);
Yang = angle(Y);
m=2*nb-2;
JC=zeros(m,m);
% Computing the elements of the jacobian matrix--
iter = 0;
nc=nb-1;
for i=1:nc
DP(i)=100;
DQ(i)=100;
end
% Checking the convergence--
while (any(abs(DP)> tol) | any(abs(DQ) > tol))
iter=iter+1;
for i=1:nb
Vm2(i)=Vm(i);
end
for i=1:nc
if bustype(i) == 2
Vm2(i)=abs(V(i));
end
PI=PG(i)-PL(i);
DP(i)=0;
DQ(i)=0;
JC(i,i)=0;
JC(i,i+nc)=0;
JC(i+nc,i)=0;
JC(i+nc,i+nc)=0;
for j=1:nb
DP(i)=DP(i)-Ym(i,j)*Vm(i)*Vm(j)*cos(delta(i)-delta(j)-Yang(i,j));
DQ(i)=DQ(i)-Ym(i,j)*Vm(i)*Vm(j)*sin(delta(i)-delta(j)-Yang(i,j));
if j==i
continue
else
JC(i,i)=JC(i,i)+Vm(i)*Vm(j)*Ym(i,j)*sin(delta(i)-delta(j)-Yang(i,j));
JC(i,i+nc)=JC(i,i+nc)-Vm(j)*Ym(i,j)*cos(delta(i)-delta(j)-Yang(i,j));
JC(i+nc,i)=JC(i+nc,i)-Vm(i)*Vm(j)*Ym(i,j)*cos(delta(i)-delta(j)-Yang(i,j));
JC(i+nc,i+nc)=JC(i+nc,i+nc)-Vm(j)*Ym(i,j)*sin(delta(i)-delta(j)-Yang(i,j));
end
end
DP(i)=DP(i)+PI;
if bustype(i)==2
Q(i)=0;
for j=1:nb
Q(i)=Q(i)+Ym(i,j)*Vm2(i)*Vm2(j)*sin(delta(i)-delta(j)-Yang(i,j));
end
Qmax1=Qmax(i)-QL(i);
Qmin1=Qmin(i)-QL(i);
if Q(i) > Qmax1
DQ(i)=DQ(i)+Qmax1;
elseif Q(i) < Qmin1
DQ(i)=DQ(i)+Qmin1;
else
DQ(i)=0;
end
else
DQ(i)=DQ(i)-QL(i);
end
JC(i,i+nc)=JC(i,i+nc)-2*Vm(i)*Ym(i,i)*cos(-Yang(i,i));
JC(i+nc,i+nc)=JC(i+nc,i+nc)-2*Vm(i)*Ym(i,i)*sin(-Yang(i,i));
for k=1:nc
if k ==i
continue
else
JC(i,k)=-Vm(i)*Vm(k)*Ym(i,k)*sin(delta(i)-delta(k)-Yang(i,k));
JC(i,k+nc)=-Vm(i)*Ym(i,k)*cos(delta(i)-delta(k)-Yang(i,k));
JC(i+nc,k)=Vm(i)*Vm(k)*Ym(i,k)*cos(delta(i)-delta(k)-Yang(i,k));
JC(i+nc,k+nc)=-Vm(i)*Ym(i,k)*sin(delta(i)-delta(k)-Yang(i,k));
end
end
end
% Modifying the jacobian matrix for non- violating P-V buses
for i=1:nb-1
for j=1:nb-1
if bustype(i)==2 & DQ(i)==0
JC(i+nb-1,j)=0;
JC(j,i+nb-1)=0;
elseif i ~= j
JC(i+nb-1,j+nb-1)=0;
JC(j+nb-1,i+nb-1)=0;
end
 
end
end
% Calculating bus voltage angles and magnitudes
J=inv(JC);
for i=1:nb-1
for j=1:nb-1
delta(i)=delta(i)-J(i,j)*DP(j)-J(i,j+nb-1)*DQ(j);
Vm(i)=Vm(i)-J(i+nb-1,j)*DP(j)-J(i+nb-1,j+nb-1)*DQ(j);
end
if DQ(i) == 0
Vm(i)=abs(V(i));
end
end
end
% Initialiazing the real and reactive power flow vectors
for j=1:nb
for i=1:nb
PS(i,j)=0;
QS(i,j)=0;
end
end
% Calculating the real and reactive bus power
for i=1:nb
for j=1:nb
PS(i,i)=PS(i,i)+Ym(i,j)*Vm(i)*Vm(j)*cos(delta(i)-delta(j)-Yang(i,j));
QS(i,i)=QS(i,i)+Ym(i,j)*Vm(i)*Vm(j)*sin(delta(i)-delta(j)-Yang(i,j));
end
end
% Calculating the real and reactive power flow
for j=1:nb
for i=j+1:nb
if Ym(i,j)~=0
PS(i,j)=Ym(i,j)*Vm(i)*Vm(j)*cos(delta(i)-delta(j)-Yang(i,j))-Ym(i,j)*Vm(i)^2*cos(Yang(i,j));
PS(j,i)=Ym(j,i)*Vm(j)*Vm(i)*cos(delta(j)-delta(i)-Yang(j,i))-Ym(j,i)*Vm(j)^2*cos(Yang(j,i));
QS(i,j)=Ym(i,j)*Vm(i)*Vm(j)*sin(delta(i)-delta(j)-Yang(i,j))+Ym(i,j)*Vm(i)^2*sin(Yang(i,j));
QS(j,i)=Ym(j,i)*Vm(j)*Vm(i)*sin(delta(j)-delta(i)-Yang(j,i))+Ym(j,i)*Vm(j)^2*sin(Yang(j,i));
end
end
end
for i=1:nb
Ps(i,i)=0; Qs(i,i)=0;
if bustype(i) == 3
Ps(i,i)= PS(i,i);
Qs(i,i)= QS(i,i);
else
if bustype(i) == 2
Qs(i,i)=QS(i,i);
Ps(i,i)=PG(i);
end
end
end
for i=1:nbr
PSr(i)=PS(fbus(i),tbus(i));
QSr(i)=QS(fbus(i),tbus(i));
PSl(i)=PS(tbus(i),fbus(i));
QSl(i)=QS(tbus(i),fbus(i));
PSt(i)=PS(fbus(i),tbus(i))+PS(tbus(i),fbus(i));
QSt(i)=QS(fbus(i),tbus(i))+QS(tbus(i),fbus(i));
end
deltad=180/pi*delta;
bus(:,8)=Vm; cd=bus(:,2);
 
Pgt= sum(Ps)*baseMVA; Qgt = sum(Qs)*baseMVA;
Plt=sum(PL)*baseMVA;Qlt = sum(QL)*baseMVA;
PStt=sum(PSt);
QStt=sum(QSt);
toc;
z=toc;

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
function plotcurv(data_name,tol)
if nargin < 2
tol = 0.001;
end
[baseMVA,bus,gen,branch,Vm]=pcal(data_name,tol);
busno =input('Bus no = ');
per=input('Period = ');
limit=input('Limit = ');
m=input('For active power = 3 For reactive power = 4 = ');
x=bus(busno,m);
i=(limit-x)/per;
i=i+1;
V(1)=Vm(:,busno);
z=(limit-x)/baseMVA;
for s=2:i
bus(busno,m)=bus(busno,m)+per;
[bus,gen,branch,ext]=intnum(bus,gen,branch);
[Y]=network_f(branch,bus,baseMVA);
[Vm]=newtraph(bus,gen,Y,baseMVA,branch,tol);
[bus,gen,branch]=extnum(bus,gen,branch,ext);
V(s)=Vm(:,busno);
end
a = (x:per:(limit))/baseMVA;
t=1:i;
c = V(t);
plot(a,c,'-rs','LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',5)
if m == 3
title(['P-V Curve for Real Power Injection of ',num2str(z),' p.u. at Bus',num2str(busno)])
else
title(['Q-V Curve for Reactive Power Injection of ',num2str(z),' p.u. at Bus ',num2str(busno)])
end
if m == 3
xlabel('P (p.u.)')
else
xlabel('Q (p.u.)')
end
ylabel('V (p.u.)')
grid on

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
function [baseMVA,bus,gen,branch,Vm] = pcal(data_name,tol)
 
global nb nbr
[baseMVA,bus,gen,branch] = feval(data_name);
[bus,gen,branch,ext]=intnum(bus,gen,branch); 
[Y]=network_f(branch,bus,baseMVA);
% If number of inputs is less than 1, take tolerance tol = 0.001.
% Otherwise, tol is given as an input.
if nargin < 2
tol = 0.001;
end
Vm = newtraph(bus,gen,Y,baseMVA,branch,tol);
[bus,gen,branch]=extnum(bus,gen,branch,ext);