1 /* function subroutines used by power.c */ 2 3 #include "power.h" 4 5 PetscErrorCode GetListofEdges_Power(PFDATA *pfdata,PetscInt *edgelist) 6 { 7 PetscErrorCode ierr; 8 PetscInt i,fbus,tbus,nbranches=pfdata->nbranch; 9 EDGE_Power branch=pfdata->branch; 10 PetscBool netview=PETSC_FALSE; 11 12 PetscFunctionBegin; 13 ierr = PetscOptionsHasName(NULL,NULL, "-powernet_view",&netview);CHKERRQ(ierr); 14 for (i=0; i<nbranches; i++) { 15 fbus = branch[i].internal_i; 16 tbus = branch[i].internal_j; 17 edgelist[2*i] = fbus; 18 edgelist[2*i+1] = tbus; 19 if (netview) { 20 ierr = PetscPrintf(PETSC_COMM_SELF,"branch %d, bus[%d] -> bus[%d]\n",i,fbus,tbus);CHKERRQ(ierr); 21 } 22 } 23 if (netview) { 24 for (i=0; i<pfdata->nbus; i++) { 25 if (pfdata->bus[i].ngen) { 26 ierr = PetscPrintf(PETSC_COMM_SELF," bus %D: gen\n",i);CHKERRQ(ierr); 27 } else if (pfdata->bus[i].nload) { 28 ierr = PetscPrintf(PETSC_COMM_SELF," bus %D: load\n",i);CHKERRQ(ierr); 29 } 30 } 31 } 32 PetscFunctionReturn(0); 33 } 34 35 PetscErrorCode FormJacobian_Power_private(DM networkdm,Vec localX,Mat J,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx) 36 { 37 PetscErrorCode ierr; 38 const PetscScalar *xarr; 39 PetscInt i,v,row[2],col[8],e,vfrom,vto; 40 PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto,numComps; 41 PetscScalar values[8]; 42 PetscInt j,key,offset,goffset; 43 PetscScalar Vm; 44 UserCtx_Power *user_power=(UserCtx_Power*)appctx; 45 PetscScalar Sbase=user_power->Sbase; 46 VERTEX_Power bus; 47 PetscBool ghostvtex; 48 void* component; 49 50 PetscFunctionBegin; 51 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 52 53 for (v=0; v<nv; v++) { 54 ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr); 55 56 ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr); 57 for (j = 0; j < numComps; j++) { 58 ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr); 59 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vtx[v],&goffset);CHKERRQ(ierr); 60 ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr); 61 62 if (key == user_power->compkey_bus) { 63 PetscInt nconnedges; 64 const PetscInt *connedges; 65 66 bus = (VERTEX_Power)(component); 67 if (!ghostvtex) { 68 /* Handle reference bus constrained dofs */ 69 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 70 row[0] = goffset; row[1] = goffset+1; 71 col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1; 72 values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0; 73 ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr); 74 break; 75 } 76 77 Vm = xarr[offset+1]; 78 79 /* Shunt injections */ 80 row[0] = goffset; row[1] = goffset+1; 81 col[0] = goffset; col[1] = goffset+1; 82 values[0] = values[1] = values[2] = values[3] = 0.0; 83 if (bus->ide != PV_BUS) { 84 values[1] = 2.0*Vm*bus->gl/Sbase; 85 values[3] = -2.0*Vm*bus->bl/Sbase; 86 } 87 ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr); 88 } 89 90 ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr); 91 92 for (i=0; i < nconnedges; i++) { 93 EDGE_Power branch; 94 VERTEX_Power busf,bust; 95 PetscInt keyf,keyt; 96 PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt; 97 const PetscInt *cone; 98 PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf; 99 100 e = connedges[i]; 101 ierr = DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch);CHKERRQ(ierr); 102 if (!branch->status) continue; 103 104 Gff = branch->yff[0]; 105 Bff = branch->yff[1]; 106 Gft = branch->yft[0]; 107 Bft = branch->yft[1]; 108 Gtf = branch->ytf[0]; 109 Btf = branch->ytf[1]; 110 Gtt = branch->ytt[0]; 111 Btt = branch->ytt[1]; 112 113 ierr = DMNetworkGetConnectedVertices(networkdm,edges[e],&cone);CHKERRQ(ierr); 114 vfrom = cone[0]; 115 vto = cone[1]; 116 117 ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr); 118 ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr); 119 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);CHKERRQ(ierr); 120 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);CHKERRQ(ierr); 121 122 if (goffsetto < 0) goffsetto = -goffsetto - 1; 123 124 thetaf = xarr[offsetfrom]; 125 Vmf = xarr[offsetfrom+1]; 126 thetat = xarr[offsetto]; 127 Vmt = xarr[offsetto+1]; 128 thetaft = thetaf - thetat; 129 thetatf = thetat - thetaf; 130 131 ierr = DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf);CHKERRQ(ierr); 132 ierr = DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust);CHKERRQ(ierr); 133 134 if (vfrom == vtx[v]) { 135 if (busf->ide != REF_BUS) { 136 /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */ 137 row[0] = goffsetfrom; 138 col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1; 139 values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */ 140 values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */ 141 values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */ 142 values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */ 143 144 ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr); 145 } 146 if (busf->ide != PV_BUS && busf->ide != REF_BUS) { 147 row[0] = goffsetfrom+1; 148 col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1; 149 /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */ 150 values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft)); 151 values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); 152 values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft)); 153 values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); 154 155 ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr); 156 } 157 } else { 158 if (bust->ide != REF_BUS) { 159 row[0] = goffsetto; 160 col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1; 161 /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */ 162 values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */ 163 values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */ 164 values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */ 165 values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */ 166 167 ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr); 168 } 169 if (bust->ide != PV_BUS && bust->ide != REF_BUS) { 170 row[0] = goffsetto+1; 171 col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1; 172 /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */ 173 values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf)); 174 values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); 175 values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf)); 176 values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); 177 178 ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr); 179 } 180 } 181 } 182 if (!ghostvtex && bus->ide == PV_BUS) { 183 row[0] = goffset+1; col[0] = goffset+1; values[0] = 1.0; 184 if (user_power->jac_error) values[0] = 50.0; 185 ierr = MatSetValues(J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); 186 } 187 } 188 } 189 } 190 191 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode FormJacobian_Power(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx) 196 { 197 PetscErrorCode ierr; 198 DM networkdm; 199 Vec localX; 200 PetscInt nv,ne; 201 const PetscInt *vtx,*edges; 202 203 PetscFunctionBegin; 204 ierr = MatZeroEntries(J);CHKERRQ(ierr); 205 206 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 207 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 208 209 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 210 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 211 212 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 213 ierr = FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);CHKERRQ(ierr); 214 215 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 216 217 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 218 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 PetscErrorCode FormFunction_Power(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx) 223 { 224 PetscErrorCode ierr; 225 UserCtx_Power *User=(UserCtx_Power*)appctx; 226 PetscInt e,v,vfrom,vto; 227 const PetscScalar *xarr; 228 PetscScalar *farr; 229 PetscInt offsetfrom,offsetto,offset,i,j,key,numComps; 230 PetscScalar Vm; 231 PetscScalar Sbase=User->Sbase; 232 VERTEX_Power bus=NULL; 233 GEN gen; 234 LOAD load; 235 PetscBool ghostvtex; 236 void* component; 237 238 PetscFunctionBegin; 239 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 240 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 241 242 for (v=0; v<nv; v++) { 243 ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr); 244 ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr); 245 ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr); 246 247 for (j = 0; j < numComps; j++) { 248 ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr); 249 if (key == User->compkey_bus) { 250 PetscInt nconnedges; 251 const PetscInt *connedges; 252 253 bus = (VERTEX_Power)(component); 254 /* Handle reference bus constrained dofs */ 255 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 256 farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0; 257 farr[offset+1] = xarr[offset+1] - bus->vm; 258 break; 259 } 260 261 if (!ghostvtex) { 262 Vm = xarr[offset+1]; 263 264 /* Shunt injections */ 265 farr[offset] += Vm*Vm*bus->gl/Sbase; 266 if (bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase; 267 } 268 269 ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr); 270 for (i=0; i < nconnedges; i++) { 271 EDGE_Power branch; 272 PetscInt keye; 273 PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt; 274 const PetscInt *cone; 275 PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf; 276 277 e = connedges[i]; 278 ierr = DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);CHKERRQ(ierr); 279 if (!branch->status) continue; 280 Gff = branch->yff[0]; 281 Bff = branch->yff[1]; 282 Gft = branch->yft[0]; 283 Bft = branch->yft[1]; 284 Gtf = branch->ytf[0]; 285 Btf = branch->ytf[1]; 286 Gtt = branch->ytt[0]; 287 Btt = branch->ytt[1]; 288 289 ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr); 290 vfrom = cone[0]; 291 vto = cone[1]; 292 293 ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr); 294 ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr); 295 296 thetaf = xarr[offsetfrom]; 297 Vmf = xarr[offsetfrom+1]; 298 thetat = xarr[offsetto]; 299 Vmt = xarr[offsetto+1]; 300 thetaft = thetaf - thetat; 301 thetatf = thetat - thetaf; 302 303 if (vfrom == vtx[v]) { 304 farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); 305 farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); 306 } else { 307 farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); 308 farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); 309 } 310 } 311 } else if (key == User->compkey_gen) { 312 if (!ghostvtex) { 313 gen = (GEN)(component); 314 if (!gen->status) continue; 315 farr[offset] += -gen->pg/Sbase; 316 farr[offset+1] += -gen->qg/Sbase; 317 } 318 } else if (key == User->compkey_load) { 319 if (!ghostvtex) { 320 load = (LOAD)(component); 321 farr[offset] += load->pl/Sbase; 322 farr[offset+1] += load->ql/Sbase; 323 } 324 } 325 } 326 if (bus && bus->ide == PV_BUS) { 327 farr[offset+1] = xarr[offset+1] - bus->vm; 328 } 329 } 330 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 331 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 PetscErrorCode SetInitialGuess_Power(DM networkdm,Vec localX,PetscInt nv,PetscInt ne,const PetscInt *vtx,const PetscInt *edges,void* appctx) 336 { 337 PetscErrorCode ierr; 338 VERTEX_Power bus; 339 PetscInt i; 340 GEN gen; 341 PetscBool ghostvtex; 342 PetscScalar *xarr; 343 PetscInt key,numComps,j,offset; 344 void* component; 345 PetscMPIInt rank; 346 MPI_Comm comm; 347 UserCtx_Power *User=(UserCtx_Power*)appctx; 348 349 PetscFunctionBegin; 350 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 351 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 352 ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); 353 for (i = 0; i < nv; i++) { 354 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 355 if (ghostvtex) continue; 356 357 ierr = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr); 358 ierr = DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);CHKERRQ(ierr); 359 for (j=0; j < numComps; j++) { 360 ierr = DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component);CHKERRQ(ierr); 361 if (key == User->compkey_bus) { 362 bus = (VERTEX_Power)(component); 363 xarr[offset] = bus->va*PETSC_PI/180.0; 364 xarr[offset+1] = bus->vm; 365 } else if (key == User->compkey_gen) { 366 gen = (GEN)(component); 367 if (!gen->status) continue; 368 xarr[offset+1] = gen->vs; 369 break; 370 } 371 } 372 } 373 ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); 374 PetscFunctionReturn(0); 375 } 376