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