xref: /petsc/src/snes/tutorials/network/power/PFReadData.c (revision 862e4a309d45a165aaa4da0d704ba733429d833a)
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;
59   Gen    = pf->gen;
60   Load   = pf->load;
61   Branch = pf->branch;
62 
63   /* Setting pf->sbase to 100 */
64   pf->sbase = 100.0;
65 
66   PetscCall(PetscMalloc1(maxbusnum + 1, &busext2intmap));
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     PetscCheck(fgets(line, MAXLINE, fp), 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;
80       Bus[busi].ide                  = ide;
81       Bus[busi].area                 = area;
82       Bus[busi].gl                   = gl;
83       Bus[busi].bl                   = bl;
84       Bus[busi].vm                   = vm;
85       Bus[busi].va                   = va;
86       Bus[busi].basekV               = basekV;
87       Bus[busi].internal_i           = busi;
88       busext2intmap[Bus[busi].bus_i] = busi;
89 
90       if (!((Pd == 0.0) && (Qd == 0.0))) {
91         Load[loadi].bus_i                 = Bus[busi].bus_i;
92         Load[loadi].status                = 1;
93         Load[loadi].pl                    = Pd;
94         Load[loadi].ql                    = Qd;
95         Load[loadi].area                  = Bus[busi].area;
96         Load[loadi].internal_i            = busi;
97         Bus[busi].lidx[Bus[busi].nload++] = loadi;
98         PetscCheck(Bus[busi].nload <= NLOAD_AT_BUS_MAX, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exceeded maximum number of loads allowed at bus");
99         loadi++;
100       }
101       busi++;
102     }
103 
104     /* Read generator data */
105     if (i >= gen_start_line && i < gen_end_line) {
106       double pg, qg, qt, qb, vs, mbase, pt, pb;
107       int    bus_i, status;
108       sscanf(line, "%d %lf %lf %lf %lf %lf %lf %d %lf %lf", &bus_i, &pg, &qg, &qt, &qb, &vs, &mbase, &status, &pt, &pb);
109       Gen[geni].bus_i  = bus_i;
110       Gen[geni].status = status;
111       Gen[geni].pg     = pg;
112       Gen[geni].qg     = qg;
113       Gen[geni].qt     = qt;
114       Gen[geni].qb     = qb;
115       Gen[geni].vs     = vs;
116       Gen[geni].mbase  = mbase;
117       Gen[geni].pt     = pt;
118       Gen[geni].pb     = pb;
119 
120       intbusnum                                  = busext2intmap[Gen[geni].bus_i];
121       Gen[geni].internal_i                       = intbusnum;
122       Bus[intbusnum].gidx[Bus[intbusnum].ngen++] = geni;
123 
124       Bus[intbusnum].vm = Gen[geni].vs;
125 
126       PetscCheck(Bus[intbusnum].ngen <= NGEN_AT_BUS_MAX, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exceeded maximum number of generators allowed at bus");
127       geni++;
128     }
129 
130     if (i >= br_start_line && i < br_end_line) {
131       PetscScalar R, X, Bc, B, G, Zm, tap, shift, tap2, tapr, tapi;
132       double      r, x, b, rateA, rateB, rateC, tapratio, phaseshift;
133       int         fbus, tbus, status;
134       sscanf(line, "%d %d %lf %lf %lf %lf %lf %lf %lf %lf %d", &fbus, &tbus, &r, &x, &b, &rateA, &rateB, &rateC, &tapratio, &phaseshift, &status);
135       Branch[bri].fbus       = fbus;
136       Branch[bri].tbus       = tbus;
137       Branch[bri].status     = status;
138       Branch[bri].r          = r;
139       Branch[bri].x          = x;
140       Branch[bri].b          = b;
141       Branch[bri].rateA      = rateA;
142       Branch[bri].rateB      = rateB;
143       Branch[bri].rateC      = rateC;
144       Branch[bri].tapratio   = tapratio;
145       Branch[bri].phaseshift = phaseshift;
146 
147       if (Branch[bri].tapratio == 0.0) Branch[bri].tapratio = 1.0;
148       Branch[bri].phaseshift *= PETSC_PI / 180.0;
149 
150       intbusnum              = busext2intmap[Branch[bri].fbus];
151       Branch[bri].internal_i = intbusnum;
152 
153       intbusnum              = busext2intmap[Branch[bri].tbus];
154       Branch[bri].internal_j = intbusnum;
155 
156       /* Compute self and transfer admittances */
157       R  = Branch[bri].r;
158       X  = Branch[bri].x;
159       Bc = Branch[bri].b;
160 
161       Zm = R * R + X * X;
162       G  = R / Zm;
163       B  = -X / Zm;
164 
165       tap   = Branch[bri].tapratio;
166       shift = Branch[bri].phaseshift;
167       tap2  = tap * tap;
168       tapr  = tap * PetscCosScalar(shift);
169       tapi  = tap * PetscSinScalar(shift);
170 
171       Branch[bri].yff[0] = G / tap2;
172       Branch[bri].yff[1] = (B + Bc / 2.0) / tap2;
173 
174       Branch[bri].yft[0] = -(G * tapr - B * tapi) / tap2;
175       Branch[bri].yft[1] = -(B * tapr + G * tapi) / tap2;
176 
177       Branch[bri].ytf[0] = -(G * tapr + B * tapi) / tap2;
178       Branch[bri].ytf[1] = -(B * tapr - G * tapi) / tap2;
179 
180       Branch[bri].ytt[0] = G;
181       Branch[bri].ytt[1] = B + Bc / 2.0;
182 
183       bri++;
184     }
185   }
186   fclose(fp);
187 
188   /* Reorder the generator data structure according to bus numbers */
189   genj  = 0;
190   loadj = 0;
191   PetscCall(PetscMalloc1(pf->ngen, &newgen));
192   PetscCall(PetscMalloc1(pf->nload, &newload));
193   for (i = 0; i < pf->nbus; i++) {
194     for (j = 0; j < pf->bus[i].ngen; j++) PetscCall(PetscMemcpy(&newgen[genj++], &pf->gen[pf->bus[i].gidx[j]], sizeof(struct _p_GEN)));
195     for (j = 0; j < pf->bus[i].nload; j++) PetscCall(PetscMemcpy(&newload[loadj++], &pf->load[pf->bus[i].lidx[j]], sizeof(struct _p_LOAD)));
196   }
197   PetscCall(PetscFree(pf->gen));
198   PetscCall(PetscFree(pf->load));
199   pf->gen  = newgen;
200   pf->load = newload;
201 
202   PetscCall(PetscFree(busext2intmap));
203   PetscFunctionReturn(0);
204 }
205