xref: /petsc/src/snes/tutorials/network/power/PFReadData.c (revision 38403884fb32dc03829b32f29c201f60a29d7265)
1 #include "petscmat.h"
2 #include "power.h"
3 #include <string.h>
4 #include <ctype.h>
5 
6 PetscErrorCode PFReadMatPowerData(PFDATA *pf,char *filename)
7 {
8   FILE           *fp;
9   VERTEX_Power   Bus;
10   LOAD           Load;
11   GEN            Gen;
12   EDGE_Power     Branch;
13   PetscInt       line_counter=0;
14   PetscInt       bus_start_line=-1,bus_end_line=-1; /* xx_end_line points to the next line after the record ends */
15   PetscInt       gen_start_line=-1,gen_end_line=-1;
16   PetscInt       br_start_line=-1,br_end_line=-1;
17   char           line[MAXLINE];
18   PetscInt       loadi=0,geni=0,bri=0,busi=0,i,j;
19   int            extbusnum,bustype_i;
20   double         Pd,Qd;
21   PetscInt       maxbusnum=-1,intbusnum,*busext2intmap,genj,loadj;
22   GEN            newgen;
23   LOAD           newload;
24 
25   PetscFunctionBegin;
26   fp = fopen(filename,"r");
27   /* Check for valid file */
28   PetscCheck(fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Can't open Matpower data file %s",filename);
29   pf->nload=0;
30   while (fgets(line,MAXLINE,fp)) {
31     if (strstr(line,"mpc.bus = ["))    bus_start_line = line_counter+1; /* Bus data starts from next line */
32     if (strstr(line,"mpc.gen") && gen_start_line == -1)    gen_start_line = line_counter+1; /* Generator data starts from next line */
33     if (strstr(line,"mpc.branch")) br_start_line = line_counter+1; /* Branch data starts from next line */
34     if (strstr(line,"];")) {
35       if (bus_start_line != -1 && bus_end_line == -1) bus_end_line = line_counter;
36       if (gen_start_line != -1 && gen_end_line == -1) gen_end_line = line_counter;
37       if (br_start_line  != -1 && br_end_line == -1) br_end_line = line_counter;
38     }
39 
40     /* Count the number of pq loads */
41     if (bus_start_line != -1 && line_counter >= bus_start_line && bus_end_line == -1) {
42       sscanf(line,"%d %d %lf %lf",&extbusnum,&bustype_i,&Pd,&Qd);
43       if (!((Pd == 0.0) && (Qd == 0.0))) pf->nload++;
44       if (extbusnum > maxbusnum) maxbusnum = extbusnum;
45     }
46     line_counter++;
47   }
48   fclose(fp);
49 
50   pf->nbus    = bus_end_line - bus_start_line;
51   pf->ngen    = gen_end_line - gen_start_line;
52   pf->nbranch = br_end_line  - br_start_line;
53 
54   PetscCall(PetscCalloc1(pf->nbus,&pf->bus));
55   PetscCall(PetscCalloc1(pf->ngen,&pf->gen));
56   PetscCall(PetscCalloc1(pf->nload,&pf->load));
57   PetscCall(PetscCalloc1(pf->nbranch,&pf->branch));
58   Bus = pf->bus; Gen = pf->gen; Load = pf->load; Branch = pf->branch;
59 
60   /* Setting pf->sbase to 100 */
61   pf->sbase = 100.0;
62 
63   PetscCall(PetscMalloc1(maxbusnum+1,&busext2intmap));
64   for (i=0; i < maxbusnum+1; i++) busext2intmap[i] = -1;
65 
66   fp = fopen(filename,"r");
67   /* Reading data */
68   for (i=0;i<line_counter;i++) {
69     PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_SUP,"File is incorrectly formatted");
70 
71     if ((i >= bus_start_line) && (i < bus_end_line)) {
72       double gl,bl,vm,va,basekV;
73       int    bus_i,ide,area;
74       /* Bus data */
75       sscanf(line,"%d %d %lf %lf %lf %lf %d %lf %lf %lf",&bus_i,&ide,&Pd,&Qd,&gl,&bl,&area,&vm,&va,&basekV);
76       Bus[busi].bus_i = bus_i; Bus[busi].ide = ide; Bus[busi].area = area;
77       Bus[busi].gl = gl; Bus[busi].bl = bl;
78       Bus[busi].vm = vm; Bus[busi].va = va; Bus[busi].basekV = basekV;
79       Bus[busi].internal_i = busi;
80       busext2intmap[Bus[busi].bus_i] = busi;
81 
82       if (!((Pd == 0.0) && (Qd == 0.0))) {
83         Load[loadi].bus_i = Bus[busi].bus_i;
84         Load[loadi].status = 1;
85         Load[loadi].pl = Pd;
86         Load[loadi].ql = Qd;
87         Load[loadi].area = Bus[busi].area;
88         Load[loadi].internal_i = busi;
89         Bus[busi].lidx[Bus[busi].nload++] = loadi;
90         PetscCheck(Bus[busi].nload <= NLOAD_AT_BUS_MAX,PETSC_COMM_SELF,PETSC_ERR_SUP,"Exceeded maximum number of loads allowed at bus");
91         loadi++;
92       }
93       busi++;
94     }
95 
96     /* Read generator data */
97     if (i >= gen_start_line && i < gen_end_line) {
98       double pg,qg,qt,qb,vs,mbase,pt,pb;
99       int    bus_i,status;
100       sscanf(line,"%d %lf %lf %lf %lf %lf %lf %d %lf %lf",&bus_i,&pg,&qg,&qt,&qb,&vs,&mbase,&status,&pt,&pb);
101       Gen[geni].bus_i = bus_i; Gen[geni].status = status;
102       Gen[geni].pg = pg; Gen[geni].qg = qg; Gen[geni].qt = qt; Gen[geni].qb = qb;
103       Gen[geni].vs = vs; Gen[geni].mbase = mbase; Gen[geni].pt = pt; Gen[geni].pb = pb;
104 
105       intbusnum = busext2intmap[Gen[geni].bus_i];
106       Gen[geni].internal_i = intbusnum;
107       Bus[intbusnum].gidx[Bus[intbusnum].ngen++] = geni;
108 
109       Bus[intbusnum].vm = Gen[geni].vs;
110 
111       PetscCheck(Bus[intbusnum].ngen <= NGEN_AT_BUS_MAX,PETSC_COMM_SELF,PETSC_ERR_SUP,"Exceeded maximum number of generators allowed at bus");
112       geni++;
113     }
114 
115     if (i >= br_start_line && i < br_end_line) {
116       PetscScalar R,X,Bc,B,G,Zm,tap,shift,tap2,tapr,tapi;
117       double      r,x,b,rateA,rateB,rateC,tapratio,phaseshift;
118       int         fbus,tbus,status;
119       sscanf(line,"%d %d %lf %lf %lf %lf %lf %lf %lf %lf %d",&fbus,&tbus,&r,&x,&b,&rateA,&rateB,&rateC,&tapratio,&phaseshift,&status);
120       Branch[bri].fbus = fbus; Branch[bri].tbus = tbus; Branch[bri].status = status;
121       Branch[bri].r = r; Branch[bri].x = x; Branch[bri].b = b;
122       Branch[bri].rateA = rateA; Branch[bri].rateB = rateB; Branch[bri].rateC = rateC;
123       Branch[bri].tapratio = tapratio; Branch[bri].phaseshift = phaseshift;
124 
125       if (Branch[bri].tapratio == 0.0) Branch[bri].tapratio = 1.0;
126       Branch[bri].phaseshift *= PETSC_PI/180.0;
127 
128       intbusnum = busext2intmap[Branch[bri].fbus];
129       Branch[bri].internal_i = intbusnum;
130 
131       intbusnum = busext2intmap[Branch[bri].tbus];
132       Branch[bri].internal_j = intbusnum;
133 
134       /* Compute self and transfer admittances */
135       R = Branch[bri].r;
136       X = Branch[bri].x;
137       Bc = Branch[bri].b;
138 
139       Zm = R*R + X*X;
140       G  = R/Zm;
141       B  = -X/Zm;
142 
143       tap = Branch[bri].tapratio;
144       shift = Branch[bri].phaseshift;
145       tap2 = tap*tap;
146       tapr = tap*PetscCosScalar(shift);
147       tapi = tap*PetscSinScalar(shift);
148 
149       Branch[bri].yff[0] = G/tap2;
150       Branch[bri].yff[1] = (B+Bc/2.0)/tap2;
151 
152       Branch[bri].yft[0] = -(G*tapr - B*tapi)/tap2;
153       Branch[bri].yft[1] = -(B*tapr + G*tapi)/tap2;
154 
155       Branch[bri].ytf[0] = -(G*tapr + B*tapi)/tap2;
156       Branch[bri].ytf[1] = -(B*tapr - G*tapi)/tap2;
157 
158       Branch[bri].ytt[0] = G;
159       Branch[bri].ytt[1] = B+Bc/2.0;
160 
161       bri++;
162     }
163   }
164   fclose(fp);
165 
166   /* Reorder the generator data structure according to bus numbers */
167   genj=0; loadj=0;
168   PetscCall(PetscMalloc1(pf->ngen,&newgen));
169   PetscCall(PetscMalloc1(pf->nload,&newload));
170   for (i = 0; i < pf->nbus; i++) {
171     for (j = 0; j < pf->bus[i].ngen; j++) {
172       PetscCall(PetscMemcpy(&newgen[genj++],&pf->gen[pf->bus[i].gidx[j]],sizeof(struct _p_GEN)));
173     }
174     for (j = 0; j < pf->bus[i].nload; j++) {
175       PetscCall(PetscMemcpy(&newload[loadj++],&pf->load[pf->bus[i].lidx[j]],sizeof(struct _p_LOAD)));
176     }
177   }
178   PetscCall(PetscFree(pf->gen));
179   PetscCall(PetscFree(pf->load));
180   pf->gen = newgen;
181   pf->load = newload;
182 
183   PetscCall(PetscFree(busext2intmap));
184   PetscFunctionReturn(0);
185 }
186