1 static char help[] = "Application of adjoint sensitivity analysis for power grid stability analysis of WECC 9 bus system.\n\ 2 This example is based on the 9-bus (node) example given in the book Power\n\ 3 Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\ 4 The power grid in this example consists of 9 buses (nodes), 3 generators,\n\ 5 3 loads, and 9 transmission lines. The network equations are written\n\ 6 in current balance form using rectangular coordinates.\n\n"; 7 8 /* 9 This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS. 10 The objectivie is to find optimal parameter PG for each generator to minizie the frequency violations due to faults. 11 The problem features discontinuities and a cost function in integral form. 12 The gradient is computed with the discrete adjoint of an implicit theta method, see ex9busadj.c for details. 13 */ 14 15 #include <petsctao.h> 16 #include <petscts.h> 17 #include <petscdm.h> 18 #include <petscdmda.h> 19 #include <petscdmcomposite.h> 20 #include <petsctime.h> 21 22 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 23 24 #define freq 60 25 #define w_s (2*PETSC_PI*freq) 26 27 /* Sizes and indices */ 28 const PetscInt nbus = 9; /* Number of network buses */ 29 const PetscInt ngen = 3; /* Number of generators */ 30 const PetscInt nload = 3; /* Number of loads */ 31 const PetscInt gbus[3] = {0,1,2}; /* Buses at which generators are incident */ 32 const PetscInt lbus[3] = {4,5,7}; /* Buses at which loads are incident */ 33 34 /* Generator real and reactive powers (found via loadflow) */ 35 PetscScalar PG[3] = { 0.69,1.59,0.69}; 36 /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/ 37 38 const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588}; 39 /* Generator constants */ 40 const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */ 41 const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */ 42 const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */ 43 const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */ 44 const PetscScalar Xq[3] = {0.4360,0.8645,1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */ 45 const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */ 46 const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */ 47 const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */ 48 PetscScalar M[3]; /* M = 2*H/w_s */ 49 PetscScalar D[3]; /* D = 0.1*M */ 50 51 PetscScalar TM[3]; /* Mechanical Torque */ 52 /* Exciter system constants */ 53 const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */ 54 const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */ 55 const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */ 56 const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */ 57 const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */ 58 const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */ 59 const PetscScalar k1[3] = {0.0039,0.0039,0.0039}; 60 const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */ 61 62 PetscScalar Vref[3]; 63 /* Load constants 64 We use a composite load model that describes the load and reactive powers at each time instant as follows 65 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i 66 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i 67 where 68 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads 69 ld_alphap,ld_alphap - Percentage contribution (weights) or loads 70 P_D0 - Real power load 71 Q_D0 - Reactive power load 72 V_m(t) - Voltage magnitude at time t 73 V_m0 - Voltage magnitude at t = 0 74 ld_betap, ld_betaq - exponents describing the load model for real and reactive part 75 76 Note: All loads have the same characteristic currently. 77 */ 78 const PetscScalar PD0[3] = {1.25,0.9,1.0}; 79 const PetscScalar QD0[3] = {0.5,0.3,0.35}; 80 const PetscInt ld_nsegsp[3] = {3,3,3}; 81 const PetscScalar ld_alphap[3] = {1.0,0.0,0.0}; 82 const PetscScalar ld_betap[3] = {2.0,1.0,0.0}; 83 const PetscInt ld_nsegsq[3] = {3,3,3}; 84 const PetscScalar ld_alphaq[3] = {1.0,0.0,0.0}; 85 const PetscScalar ld_betaq[3] = {2.0,1.0,0.0}; 86 87 typedef struct { 88 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */ 89 DM dmpgrid; /* Composite DM to manage the entire power grid */ 90 Mat Ybus; /* Network admittance matrix */ 91 Vec V0; /* Initial voltage vector (Power flow solution) */ 92 PetscReal tfaulton,tfaultoff; /* Fault on and off times */ 93 PetscInt faultbus; /* Fault bus */ 94 PetscScalar Rfault; 95 PetscReal t0,tmax; 96 PetscInt neqs_gen,neqs_net,neqs_pgrid; 97 Mat Sol; /* Matrix to save solution at each time step */ 98 PetscInt stepnum; 99 PetscBool alg_flg; 100 PetscReal t; 101 IS is_diff; /* indices for differential equations */ 102 IS is_alg; /* indices for algebraic equations */ 103 PetscReal freq_u,freq_l; /* upper and lower frequency limit */ 104 PetscInt pow; /* power coefficient used in the cost function */ 105 PetscBool jacp_flg; 106 Mat J,Jacp; 107 Mat DRDU,DRDP; 108 } Userctx; 109 110 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 111 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr, PetscScalar *Fi) 112 { 113 PetscFunctionBegin; 114 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 115 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 116 PetscFunctionReturn(0); 117 } 118 119 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 120 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd, PetscScalar *Fq) 121 { 122 PetscFunctionBegin; 123 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 124 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 125 PetscFunctionReturn(0); 126 } 127 128 /* Saves the solution at each time to a matrix */ 129 PetscErrorCode SaveSolution(TS ts) 130 { 131 Userctx *user; 132 Vec X; 133 PetscScalar *mat; 134 const PetscScalar *x; 135 PetscInt idx; 136 PetscReal t; 137 138 PetscFunctionBegin; 139 CHKERRQ(TSGetApplicationContext(ts,&user)); 140 CHKERRQ(TSGetTime(ts,&t)); 141 CHKERRQ(TSGetSolution(ts,&X)); 142 idx = user->stepnum*(user->neqs_pgrid+1); 143 CHKERRQ(MatDenseGetArray(user->Sol,&mat)); 144 CHKERRQ(VecGetArrayRead(X,&x)); 145 mat[idx] = t; 146 CHKERRQ(PetscArraycpy(mat+idx+1,x,user->neqs_pgrid)); 147 CHKERRQ(MatDenseRestoreArray(user->Sol,&mat)); 148 CHKERRQ(VecRestoreArrayRead(X,&x)); 149 user->stepnum++; 150 PetscFunctionReturn(0); 151 } 152 153 PetscErrorCode SetInitialGuess(Vec X,Userctx *user) 154 { 155 Vec Xgen,Xnet; 156 PetscScalar *xgen,*xnet; 157 PetscInt i,idx=0; 158 PetscScalar Vr,Vi,IGr,IGi,Vm,Vm2; 159 PetscScalar Eqp,Edp,delta; 160 PetscScalar Efd,RF,VR; /* Exciter variables */ 161 PetscScalar Id,Iq; /* Generator dq axis currents */ 162 PetscScalar theta,Vd,Vq,SE; 163 164 PetscFunctionBegin; 165 M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s; 166 D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2]; 167 168 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 169 170 /* Network subsystem initialization */ 171 CHKERRQ(VecCopy(user->V0,Xnet)); 172 173 /* Generator subsystem initialization */ 174 CHKERRQ(VecGetArray(Xgen,&xgen)); 175 CHKERRQ(VecGetArray(Xnet,&xnet)); 176 177 for (i=0; i < ngen; i++) { 178 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 179 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 180 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 181 IGr = (Vr*PG[i] + Vi*QG[i])/Vm2; 182 IGi = (Vi*PG[i] - Vr*QG[i])/Vm2; 183 184 delta = PetscAtan2Real(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */ 185 186 theta = PETSC_PI/2.0 - delta; 187 188 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */ 189 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */ 190 191 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 192 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 193 194 Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */ 195 Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */ 196 197 TM[i] = PG[i]; 198 199 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */ 200 xgen[idx] = Eqp; 201 xgen[idx+1] = Edp; 202 xgen[idx+2] = delta; 203 xgen[idx+3] = w_s; 204 205 idx = idx + 4; 206 207 xgen[idx] = Id; 208 xgen[idx+1] = Iq; 209 210 idx = idx + 2; 211 212 /* Exciter */ 213 Efd = Eqp + (Xd[i] - Xdp[i])*Id; 214 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 215 VR = KE[i]*Efd + SE; 216 RF = KF[i]*Efd/TF[i]; 217 218 xgen[idx] = Efd; 219 xgen[idx+1] = RF; 220 xgen[idx+2] = VR; 221 222 Vref[i] = Vm + (VR/KA[i]); 223 224 idx = idx + 3; 225 } 226 227 CHKERRQ(VecRestoreArray(Xgen,&xgen)); 228 CHKERRQ(VecRestoreArray(Xnet,&xnet)); 229 230 /* CHKERRQ(VecView(Xgen,0)); */ 231 CHKERRQ(DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet)); 232 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 233 PetscFunctionReturn(0); 234 } 235 236 PetscErrorCode InitialGuess(Vec X,Userctx *user, const PetscScalar PGv[]) 237 { 238 Vec Xgen,Xnet; 239 PetscScalar *xgen,*xnet; 240 PetscInt i,idx=0; 241 PetscScalar Vr,Vi,IGr,IGi,Vm,Vm2; 242 PetscScalar Eqp,Edp,delta; 243 PetscScalar Efd,RF,VR; /* Exciter variables */ 244 PetscScalar Id,Iq; /* Generator dq axis currents */ 245 PetscScalar theta,Vd,Vq,SE; 246 247 PetscFunctionBegin; 248 M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s; 249 D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2]; 250 251 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 252 253 /* Network subsystem initialization */ 254 CHKERRQ(VecCopy(user->V0,Xnet)); 255 256 /* Generator subsystem initialization */ 257 CHKERRQ(VecGetArray(Xgen,&xgen)); 258 CHKERRQ(VecGetArray(Xnet,&xnet)); 259 260 for (i=0; i < ngen; i++) { 261 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 262 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 263 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 264 IGr = (Vr*PGv[i] + Vi*QG[i])/Vm2; 265 IGi = (Vi*PGv[i] - Vr*QG[i])/Vm2; 266 267 delta = PetscAtan2Real(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */ 268 269 theta = PETSC_PI/2.0 - delta; 270 271 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */ 272 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */ 273 274 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 275 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 276 277 Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */ 278 Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */ 279 280 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */ 281 xgen[idx] = Eqp; 282 xgen[idx+1] = Edp; 283 xgen[idx+2] = delta; 284 xgen[idx+3] = w_s; 285 286 idx = idx + 4; 287 288 xgen[idx] = Id; 289 xgen[idx+1] = Iq; 290 291 idx = idx + 2; 292 293 /* Exciter */ 294 Efd = Eqp + (Xd[i] - Xdp[i])*Id; 295 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 296 VR = KE[i]*Efd + SE; 297 RF = KF[i]*Efd/TF[i]; 298 299 xgen[idx] = Efd; 300 xgen[idx+1] = RF; 301 xgen[idx+2] = VR; 302 303 idx = idx + 3; 304 } 305 306 CHKERRQ(VecRestoreArray(Xgen,&xgen)); 307 CHKERRQ(VecRestoreArray(Xnet,&xnet)); 308 309 /* CHKERRQ(VecView(Xgen,0)); */ 310 CHKERRQ(DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet)); 311 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 312 PetscFunctionReturn(0); 313 } 314 315 PetscErrorCode DICDPFiniteDifference(Vec X,Vec *DICDP, Userctx *user) 316 { 317 Vec Y; 318 PetscScalar PGv[3],eps; 319 PetscInt i,j; 320 321 PetscFunctionBegin; 322 eps = 1.e-7; 323 CHKERRQ(VecDuplicate(X,&Y)); 324 325 for (i=0;i<ngen;i++) { 326 for (j=0;j<3;j++) PGv[j] = PG[j]; 327 PGv[i] = PG[i]+eps; 328 CHKERRQ(InitialGuess(Y,user,PGv)); 329 CHKERRQ(InitialGuess(X,user,PG)); 330 331 CHKERRQ(VecAXPY(Y,-1.0,X)); 332 CHKERRQ(VecScale(Y,1./eps)); 333 CHKERRQ(VecCopy(Y,DICDP[i])); 334 } 335 CHKERRQ(VecDestroy(&Y)); 336 PetscFunctionReturn(0); 337 } 338 339 /* Computes F = [-f(x,y);g(x,y)] */ 340 PetscErrorCode ResidualFunction(SNES snes,Vec X, Vec F, Userctx *user) 341 { 342 Vec Xgen,Xnet,Fgen,Fnet; 343 PetscScalar *xgen,*xnet,*fgen,*fnet; 344 PetscInt i,idx=0; 345 PetscScalar Vr,Vi,Vm,Vm2; 346 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 347 PetscScalar Efd,RF,VR; /* Exciter variables */ 348 PetscScalar Id,Iq; /* Generator dq axis currents */ 349 PetscScalar Vd,Vq,SE; 350 PetscScalar IGr,IGi,IDr,IDi; 351 PetscScalar Zdq_inv[4],det; 352 PetscScalar PD,QD,Vm0,*v0; 353 PetscInt k; 354 355 PetscFunctionBegin; 356 CHKERRQ(VecZeroEntries(F)); 357 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 358 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet)); 359 CHKERRQ(DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet)); 360 CHKERRQ(DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet)); 361 362 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 363 The generator current injection, IG, and load current injection, ID are added later 364 */ 365 /* Note that the values in Ybus are stored assuming the imaginary current balance 366 equation is ordered first followed by real current balance equation for each bus. 367 Thus imaginary current contribution goes in location 2*i, and 368 real current contribution in 2*i+1 369 */ 370 CHKERRQ(MatMult(user->Ybus,Xnet,Fnet)); 371 372 CHKERRQ(VecGetArray(Xgen,&xgen)); 373 CHKERRQ(VecGetArray(Xnet,&xnet)); 374 CHKERRQ(VecGetArray(Fgen,&fgen)); 375 CHKERRQ(VecGetArray(Fnet,&fnet)); 376 377 /* Generator subsystem */ 378 for (i=0; i < ngen; i++) { 379 Eqp = xgen[idx]; 380 Edp = xgen[idx+1]; 381 delta = xgen[idx+2]; 382 w = xgen[idx+3]; 383 Id = xgen[idx+4]; 384 Iq = xgen[idx+5]; 385 Efd = xgen[idx+6]; 386 RF = xgen[idx+7]; 387 VR = xgen[idx+8]; 388 389 /* Generator differential equations */ 390 fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; 391 fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; 392 fgen[idx+2] = -w + w_s; 393 fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; 394 395 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 396 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 397 398 CHKERRQ(ri2dq(Vr,Vi,delta,&Vd,&Vq)); 399 /* Algebraic equations for stator currents */ 400 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 401 402 Zdq_inv[0] = Rs[i]/det; 403 Zdq_inv[1] = Xqp[i]/det; 404 Zdq_inv[2] = -Xdp[i]/det; 405 Zdq_inv[3] = Rs[i]/det; 406 407 fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 408 fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 409 410 /* Add generator current injection to network */ 411 CHKERRQ(dq2ri(Id,Iq,delta,&IGr,&IGi)); 412 413 fnet[2*gbus[i]] -= IGi; 414 fnet[2*gbus[i]+1] -= IGr; 415 416 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 417 418 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 419 420 /* Exciter differential equations */ 421 fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; 422 fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; 423 fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 424 425 idx = idx + 9; 426 } 427 428 CHKERRQ(VecGetArray(user->V0,&v0)); 429 for (i=0; i < nload; i++) { 430 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 431 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 432 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 433 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 434 PD = QD = 0.0; 435 for (k=0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 436 for (k=0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 437 438 /* Load currents */ 439 IDr = (PD*Vr + QD*Vi)/Vm2; 440 IDi = (-QD*Vr + PD*Vi)/Vm2; 441 442 fnet[2*lbus[i]] += IDi; 443 fnet[2*lbus[i]+1] += IDr; 444 } 445 CHKERRQ(VecRestoreArray(user->V0,&v0)); 446 447 CHKERRQ(VecRestoreArray(Xgen,&xgen)); 448 CHKERRQ(VecRestoreArray(Xnet,&xnet)); 449 CHKERRQ(VecRestoreArray(Fgen,&fgen)); 450 CHKERRQ(VecRestoreArray(Fnet,&fnet)); 451 452 CHKERRQ(DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet)); 453 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 454 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Fgen,&Fnet)); 455 PetscFunctionReturn(0); 456 } 457 458 /* \dot{x} - f(x,y) 459 g(x,y) = 0 460 */ 461 PetscErrorCode IFunction(TS ts,PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user) 462 { 463 SNES snes; 464 PetscScalar *f; 465 const PetscScalar *xdot; 466 PetscInt i; 467 468 PetscFunctionBegin; 469 user->t = t; 470 471 CHKERRQ(TSGetSNES(ts,&snes)); 472 CHKERRQ(ResidualFunction(snes,X,F,user)); 473 CHKERRQ(VecGetArray(F,&f)); 474 CHKERRQ(VecGetArrayRead(Xdot,&xdot)); 475 for (i=0;i < ngen;i++) { 476 f[9*i] += xdot[9*i]; 477 f[9*i+1] += xdot[9*i+1]; 478 f[9*i+2] += xdot[9*i+2]; 479 f[9*i+3] += xdot[9*i+3]; 480 f[9*i+6] += xdot[9*i+6]; 481 f[9*i+7] += xdot[9*i+7]; 482 f[9*i+8] += xdot[9*i+8]; 483 } 484 CHKERRQ(VecRestoreArray(F,&f)); 485 CHKERRQ(VecRestoreArrayRead(Xdot,&xdot)); 486 PetscFunctionReturn(0); 487 } 488 489 /* This function is used for solving the algebraic system only during fault on and 490 off times. It computes the entire F and then zeros out the part corresponding to 491 differential equations 492 F = [0;g(y)]; 493 */ 494 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx) 495 { 496 Userctx *user=(Userctx*)ctx; 497 PetscScalar *f; 498 PetscInt i; 499 500 PetscFunctionBegin; 501 CHKERRQ(ResidualFunction(snes,X,F,user)); 502 CHKERRQ(VecGetArray(F,&f)); 503 for (i=0; i < ngen; i++) { 504 f[9*i] = 0; 505 f[9*i+1] = 0; 506 f[9*i+2] = 0; 507 f[9*i+3] = 0; 508 f[9*i+6] = 0; 509 f[9*i+7] = 0; 510 f[9*i+8] = 0; 511 } 512 CHKERRQ(VecRestoreArray(F,&f)); 513 PetscFunctionReturn(0); 514 } 515 516 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user) 517 { 518 PetscInt *d_nnz; 519 PetscInt i,idx=0,start=0; 520 PetscInt ncols; 521 522 PetscFunctionBegin; 523 CHKERRQ(PetscMalloc1(user->neqs_pgrid,&d_nnz)); 524 for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0; 525 /* Generator subsystem */ 526 for (i=0; i < ngen; i++) { 527 528 d_nnz[idx] += 3; 529 d_nnz[idx+1] += 2; 530 d_nnz[idx+2] += 2; 531 d_nnz[idx+3] += 5; 532 d_nnz[idx+4] += 6; 533 d_nnz[idx+5] += 6; 534 535 d_nnz[user->neqs_gen+2*gbus[i]] += 3; 536 d_nnz[user->neqs_gen+2*gbus[i]+1] += 3; 537 538 d_nnz[idx+6] += 2; 539 d_nnz[idx+7] += 2; 540 d_nnz[idx+8] += 5; 541 542 idx = idx + 9; 543 } 544 545 start = user->neqs_gen; 546 for (i=0; i < nbus; i++) { 547 CHKERRQ(MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL)); 548 d_nnz[start+2*i] += ncols; 549 d_nnz[start+2*i+1] += ncols; 550 CHKERRQ(MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL)); 551 } 552 553 CHKERRQ(MatSeqAIJSetPreallocation(J,0,d_nnz)); 554 CHKERRQ(PetscFree(d_nnz)); 555 PetscFunctionReturn(0); 556 } 557 558 /* 559 J = [-df_dx, -df_dy 560 dg_dx, dg_dy] 561 */ 562 PetscErrorCode ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx) 563 { 564 Userctx *user=(Userctx*)ctx; 565 Vec Xgen,Xnet; 566 PetscScalar *xgen,*xnet; 567 PetscInt i,idx=0; 568 PetscScalar Vr,Vi,Vm,Vm2; 569 PetscScalar Eqp,Edp,delta; /* Generator variables */ 570 PetscScalar Efd; /* Exciter variables */ 571 PetscScalar Id,Iq; /* Generator dq axis currents */ 572 PetscScalar Vd,Vq; 573 PetscScalar val[10]; 574 PetscInt row[2],col[10]; 575 PetscInt net_start=user->neqs_gen; 576 PetscInt ncols; 577 const PetscInt *cols; 578 const PetscScalar *yvals; 579 PetscInt k; 580 PetscScalar Zdq_inv[4],det; 581 PetscScalar dVd_dVr,dVd_dVi,dVq_dVr,dVq_dVi,dVd_ddelta,dVq_ddelta; 582 PetscScalar dIGr_ddelta,dIGi_ddelta,dIGr_dId,dIGr_dIq,dIGi_dId,dIGi_dIq; 583 PetscScalar dSE_dEfd; 584 PetscScalar dVm_dVd,dVm_dVq,dVm_dVr,dVm_dVi; 585 PetscScalar PD,QD,Vm0,*v0,Vm4; 586 PetscScalar dPD_dVr,dPD_dVi,dQD_dVr,dQD_dVi; 587 PetscScalar dIDr_dVr,dIDr_dVi,dIDi_dVr,dIDi_dVi; 588 589 PetscFunctionBegin; 590 CHKERRQ(MatZeroEntries(B)); 591 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 592 CHKERRQ(DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet)); 593 594 CHKERRQ(VecGetArray(Xgen,&xgen)); 595 CHKERRQ(VecGetArray(Xnet,&xnet)); 596 597 /* Generator subsystem */ 598 for (i=0; i < ngen; i++) { 599 Eqp = xgen[idx]; 600 Edp = xgen[idx+1]; 601 delta = xgen[idx+2]; 602 Id = xgen[idx+4]; 603 Iq = xgen[idx+5]; 604 Efd = xgen[idx+6]; 605 606 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */ 607 row[0] = idx; 608 col[0] = idx; col[1] = idx+4; col[2] = idx+6; 609 val[0] = 1/ Td0p[i]; val[1] = (Xd[i] - Xdp[i])/ Td0p[i]; val[2] = -1/Td0p[i]; 610 611 CHKERRQ(MatSetValues(J,1,row,3,col,val,INSERT_VALUES)); 612 613 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */ 614 row[0] = idx + 1; 615 col[0] = idx + 1; col[1] = idx+5; 616 val[0] = 1/Tq0p[i]; val[1] = -(Xq[i] - Xqp[i])/Tq0p[i]; 617 CHKERRQ(MatSetValues(J,1,row,2,col,val,INSERT_VALUES)); 618 619 /* fgen[idx+2] = - w + w_s; */ 620 row[0] = idx + 2; 621 col[0] = idx + 2; col[1] = idx + 3; 622 val[0] = 0; val[1] = -1; 623 CHKERRQ(MatSetValues(J,1,row,2,col,val,INSERT_VALUES)); 624 625 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */ 626 row[0] = idx + 3; 627 col[0] = idx; col[1] = idx + 1; col[2] = idx + 3; col[3] = idx + 4; col[4] = idx + 5; 628 val[0] = Iq/M[i]; val[1] = Id/M[i]; val[2] = D[i]/M[i]; val[3] = (Edp + (Xqp[i]-Xdp[i])*Iq)/M[i]; val[4] = (Eqp + (Xqp[i] - Xdp[i])*Id)/M[i]; 629 CHKERRQ(MatSetValues(J,1,row,5,col,val,INSERT_VALUES)); 630 631 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 632 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 633 CHKERRQ(ri2dq(Vr,Vi,delta,&Vd,&Vq)); 634 635 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 636 637 Zdq_inv[0] = Rs[i]/det; 638 Zdq_inv[1] = Xqp[i]/det; 639 Zdq_inv[2] = -Xdp[i]/det; 640 Zdq_inv[3] = Rs[i]/det; 641 642 dVd_dVr = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta); 643 dVq_dVr = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta); 644 dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta); 645 dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta); 646 647 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */ 648 row[0] = idx+4; 649 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 650 val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0]; val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta; 651 col[3] = idx + 4; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 652 val[3] = 1; val[4] = Zdq_inv[0]*dVd_dVr + Zdq_inv[1]*dVq_dVr; val[5] = Zdq_inv[0]*dVd_dVi + Zdq_inv[1]*dVq_dVi; 653 CHKERRQ(MatSetValues(J,1,row,6,col,val,INSERT_VALUES)); 654 655 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */ 656 row[0] = idx+5; 657 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 658 val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2]; val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta; 659 col[3] = idx + 5; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 660 val[3] = 1; val[4] = Zdq_inv[2]*dVd_dVr + Zdq_inv[3]*dVq_dVr; val[5] = Zdq_inv[2]*dVd_dVi + Zdq_inv[3]*dVq_dVi; 661 CHKERRQ(MatSetValues(J,1,row,6,col,val,INSERT_VALUES)); 662 663 dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta); 664 dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta); 665 dIGr_dId = PetscSinScalar(delta); dIGr_dIq = PetscCosScalar(delta); 666 dIGi_dId = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta); 667 668 /* fnet[2*gbus[i]] -= IGi; */ 669 row[0] = net_start + 2*gbus[i]; 670 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 671 val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq; 672 CHKERRQ(MatSetValues(J,1,row,3,col,val,INSERT_VALUES)); 673 674 /* fnet[2*gbus[i]+1] -= IGr; */ 675 row[0] = net_start + 2*gbus[i]+1; 676 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 677 val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq; 678 CHKERRQ(MatSetValues(J,1,row,3,col,val,INSERT_VALUES)); 679 680 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 681 682 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */ 683 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */ 684 dSE_dEfd = k1[i]*k2[i]*PetscExpScalar(k2[i]*Efd); 685 686 row[0] = idx + 6; 687 col[0] = idx + 6; col[1] = idx + 8; 688 val[0] = (KE[i] + dSE_dEfd)/TE[i]; val[1] = -1/TE[i]; 689 CHKERRQ(MatSetValues(J,1,row,2,col,val,INSERT_VALUES)); 690 691 /* Exciter differential equations */ 692 693 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */ 694 row[0] = idx + 7; 695 col[0] = idx + 6; col[1] = idx + 7; 696 val[0] = (-KF[i]/TF[i])/TF[i]; val[1] = 1/TF[i]; 697 CHKERRQ(MatSetValues(J,1,row,2,col,val,INSERT_VALUES)); 698 699 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */ 700 /* Vm = (Vd^2 + Vq^2)^0.5; */ 701 dVm_dVd = Vd/Vm; dVm_dVq = Vq/Vm; 702 dVm_dVr = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr; 703 dVm_dVi = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi; 704 row[0] = idx + 8; 705 col[0] = idx + 6; col[1] = idx + 7; col[2] = idx + 8; 706 val[0] = (KA[i]*KF[i]/TF[i])/TA[i]; val[1] = -KA[i]/TA[i]; val[2] = 1/TA[i]; 707 col[3] = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1; 708 val[3] = KA[i]*dVm_dVr/TA[i]; val[4] = KA[i]*dVm_dVi/TA[i]; 709 CHKERRQ(MatSetValues(J,1,row,5,col,val,INSERT_VALUES)); 710 idx = idx + 9; 711 } 712 713 for (i=0; i<nbus; i++) { 714 CHKERRQ(MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals)); 715 row[0] = net_start + 2*i; 716 for (k=0; k<ncols; k++) { 717 col[k] = net_start + cols[k]; 718 val[k] = yvals[k]; 719 } 720 CHKERRQ(MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES)); 721 CHKERRQ(MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals)); 722 723 CHKERRQ(MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals)); 724 row[0] = net_start + 2*i+1; 725 for (k=0; k<ncols; k++) { 726 col[k] = net_start + cols[k]; 727 val[k] = yvals[k]; 728 } 729 CHKERRQ(MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES)); 730 CHKERRQ(MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals)); 731 } 732 733 CHKERRQ(MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY)); 734 CHKERRQ(MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY)); 735 736 CHKERRQ(VecGetArray(user->V0,&v0)); 737 for (i=0; i < nload; i++) { 738 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 739 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 740 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2= Vm*Vm; Vm4 = Vm2*Vm2; 741 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 742 PD = QD = 0.0; 743 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0; 744 for (k=0; k < ld_nsegsp[i]; k++) { 745 PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 746 dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2)); 747 dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2)); 748 } 749 for (k=0; k < ld_nsegsq[i]; k++) { 750 QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 751 dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2)); 752 dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2)); 753 } 754 755 /* IDr = (PD*Vr + QD*Vi)/Vm2; */ 756 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */ 757 758 dIDr_dVr = (dPD_dVr*Vr + dQD_dVr*Vi + PD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vr)/Vm4; 759 dIDr_dVi = (dPD_dVi*Vr + dQD_dVi*Vi + QD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vi)/Vm4; 760 761 dIDi_dVr = (-dQD_dVr*Vr + dPD_dVr*Vi - QD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vr)/Vm4; 762 dIDi_dVi = (-dQD_dVi*Vr + dPD_dVi*Vi + PD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vi)/Vm4; 763 764 /* fnet[2*lbus[i]] += IDi; */ 765 row[0] = net_start + 2*lbus[i]; 766 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 767 val[0] = dIDi_dVr; val[1] = dIDi_dVi; 768 CHKERRQ(MatSetValues(J,1,row,2,col,val,ADD_VALUES)); 769 /* fnet[2*lbus[i]+1] += IDr; */ 770 row[0] = net_start + 2*lbus[i]+1; 771 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 772 val[0] = dIDr_dVr; val[1] = dIDr_dVi; 773 CHKERRQ(MatSetValues(J,1,row,2,col,val,ADD_VALUES)); 774 } 775 CHKERRQ(VecRestoreArray(user->V0,&v0)); 776 777 CHKERRQ(VecRestoreArray(Xgen,&xgen)); 778 CHKERRQ(VecRestoreArray(Xnet,&xnet)); 779 780 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 781 782 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 783 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 784 PetscFunctionReturn(0); 785 } 786 787 /* 788 J = [I, 0 789 dg_dx, dg_dy] 790 */ 791 PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx) 792 { 793 Userctx *user=(Userctx*)ctx; 794 795 PetscFunctionBegin; 796 CHKERRQ(ResidualJacobian(snes,X,A,B,ctx)); 797 CHKERRQ(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE)); 798 CHKERRQ(MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL)); 799 PetscFunctionReturn(0); 800 } 801 802 /* 803 J = [a*I-df_dx, -df_dy 804 dg_dx, dg_dy] 805 */ 806 807 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user) 808 { 809 SNES snes; 810 PetscScalar atmp = (PetscScalar) a; 811 PetscInt i,row; 812 813 PetscFunctionBegin; 814 user->t = t; 815 816 CHKERRQ(TSGetSNES(ts,&snes)); 817 CHKERRQ(ResidualJacobian(snes,X,A,B,user)); 818 for (i=0;i < ngen;i++) { 819 row = 9*i; 820 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 821 row = 9*i+1; 822 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 823 row = 9*i+2; 824 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 825 row = 9*i+3; 826 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 827 row = 9*i+6; 828 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 829 row = 9*i+7; 830 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 831 row = 9*i+8; 832 CHKERRQ(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 833 } 834 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 835 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 836 PetscFunctionReturn(0); 837 } 838 839 /* Matrix JacobianP is constant so that it only needs to be evaluated once */ 840 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A, void *ctx0) 841 { 842 PetscScalar a; 843 PetscInt row,col; 844 Userctx *ctx=(Userctx*)ctx0; 845 846 PetscFunctionBeginUser; 847 848 if (ctx->jacp_flg) { 849 CHKERRQ(MatZeroEntries(A)); 850 851 for (col=0;col<3;col++) { 852 a = 1.0/M[col]; 853 row = 9*col+3; 854 CHKERRQ(MatSetValues(A,1,&row,1,&col,&a,INSERT_VALUES)); 855 } 856 857 ctx->jacp_flg = PETSC_FALSE; 858 859 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 860 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 861 } 862 PetscFunctionReturn(0); 863 } 864 865 static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,Userctx *user) 866 { 867 const PetscScalar *u; 868 PetscInt idx; 869 Vec Xgen,Xnet; 870 PetscScalar *r,*xgen; 871 PetscInt i; 872 873 PetscFunctionBegin; 874 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 875 CHKERRQ(DMCompositeScatter(user->dmpgrid,U,Xgen,Xnet)); 876 877 CHKERRQ(VecGetArray(Xgen,&xgen)); 878 879 CHKERRQ(VecGetArrayRead(U,&u)); 880 CHKERRQ(VecGetArray(R,&r)); 881 r[0] = 0.; 882 idx = 0; 883 for (i=0;i<ngen;i++) { 884 r[0] += PetscPowScalarInt(PetscMax(0.,PetscMax(xgen[idx+3]/(2.*PETSC_PI)-user->freq_u,user->freq_l-xgen[idx+3]/(2.*PETSC_PI))),user->pow); 885 idx += 9; 886 } 887 CHKERRQ(VecRestoreArrayRead(U,&u)); 888 CHKERRQ(VecRestoreArray(R,&r)); 889 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 890 PetscFunctionReturn(0); 891 } 892 893 static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,Userctx *user) 894 { 895 Vec Xgen,Xnet,Dgen,Dnet; 896 PetscScalar *xgen,*dgen; 897 PetscInt i; 898 PetscInt idx; 899 Vec drdu_col; 900 PetscScalar *xarr; 901 902 PetscFunctionBegin; 903 CHKERRQ(VecDuplicate(U,&drdu_col)); 904 CHKERRQ(MatDenseGetColumn(DRDU,0,&xarr)); 905 CHKERRQ(VecPlaceArray(drdu_col,xarr)); 906 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 907 CHKERRQ(DMCompositeGetLocalVectors(user->dmpgrid,&Dgen,&Dnet)); 908 CHKERRQ(DMCompositeScatter(user->dmpgrid,U,Xgen,Xnet)); 909 CHKERRQ(DMCompositeScatter(user->dmpgrid,drdu_col,Dgen,Dnet)); 910 911 CHKERRQ(VecGetArray(Xgen,&xgen)); 912 CHKERRQ(VecGetArray(Dgen,&dgen)); 913 914 idx = 0; 915 for (i=0;i<ngen;i++) { 916 dgen[idx+3] = 0.; 917 if (xgen[idx+3]/(2.*PETSC_PI) > user->freq_u) dgen[idx+3] = user->pow*PetscPowScalarInt(xgen[idx+3]/(2.*PETSC_PI)-user->freq_u,user->pow-1)/(2.*PETSC_PI); 918 if (xgen[idx+3]/(2.*PETSC_PI) < user->freq_l) dgen[idx+3] = user->pow*PetscPowScalarInt(user->freq_l-xgen[idx+3]/(2.*PETSC_PI),user->pow-1)/(-2.*PETSC_PI); 919 idx += 9; 920 } 921 922 CHKERRQ(VecRestoreArray(Dgen,&dgen)); 923 CHKERRQ(VecRestoreArray(Xgen,&xgen)); 924 CHKERRQ(DMCompositeGather(user->dmpgrid,INSERT_VALUES,drdu_col,Dgen,Dnet)); 925 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Dgen,&Dnet)); 926 CHKERRQ(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 927 CHKERRQ(VecResetArray(drdu_col)); 928 CHKERRQ(MatDenseRestoreColumn(DRDU,&xarr)); 929 CHKERRQ(VecDestroy(&drdu_col)); 930 PetscFunctionReturn(0); 931 } 932 933 static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat drdp,Userctx *user) 934 { 935 PetscFunctionBegin; 936 PetscFunctionReturn(0); 937 } 938 939 PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,Vec *DICDP,Userctx *user) 940 { 941 PetscScalar *x,*y,sensip; 942 PetscInt i; 943 944 PetscFunctionBegin; 945 CHKERRQ(VecGetArray(lambda,&x)); 946 CHKERRQ(VecGetArray(mu,&y)); 947 948 for (i=0;i<3;i++) { 949 CHKERRQ(VecDot(lambda,DICDP[i],&sensip)); 950 sensip = sensip+y[i]; 951 /* CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %D th parameter: %g \n",i,(double)sensip)); */ 952 y[i] = sensip; 953 } 954 CHKERRQ(VecRestoreArray(mu,&y)); 955 PetscFunctionReturn(0); 956 } 957 958 int main(int argc,char **argv) 959 { 960 Userctx user; 961 Vec p; 962 PetscScalar *x_ptr; 963 PetscErrorCode ierr; 964 PetscMPIInt size; 965 PetscInt i; 966 PetscViewer Xview,Ybusview; 967 PetscInt *idx2; 968 Tao tao; 969 KSP ksp; 970 PC pc; 971 Vec lowerb,upperb; 972 973 ierr = PetscInitialize(&argc,&argv,"petscoptions",help);if (ierr) return ierr; 974 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 975 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 976 977 user.jacp_flg = PETSC_TRUE; 978 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 979 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 980 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 981 982 /* Create indices for differential and algebraic equations */ 983 CHKERRQ(PetscMalloc1(7*ngen,&idx2)); 984 for (i=0; i<ngen; i++) { 985 idx2[7*i] = 9*i; idx2[7*i+1] = 9*i+1; idx2[7*i+2] = 9*i+2; idx2[7*i+3] = 9*i+3; 986 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 987 } 988 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff)); 989 CHKERRQ(ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg)); 990 CHKERRQ(PetscFree(idx2)); 991 992 /* Set run time options */ 993 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 994 { 995 user.tfaulton = 1.0; 996 user.tfaultoff = 1.2; 997 user.Rfault = 0.0001; 998 user.faultbus = 8; 999 CHKERRQ(PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL)); 1000 CHKERRQ(PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL)); 1001 CHKERRQ(PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL)); 1002 user.t0 = 0.0; 1003 user.tmax = 1.3; 1004 CHKERRQ(PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL)); 1005 CHKERRQ(PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL)); 1006 user.freq_u = 61.0; 1007 user.freq_l = 59.0; 1008 user.pow = 2; 1009 CHKERRQ(PetscOptionsReal("-frequ","","",user.freq_u,&user.freq_u,NULL)); 1010 CHKERRQ(PetscOptionsReal("-freql","","",user.freq_l,&user.freq_l,NULL)); 1011 CHKERRQ(PetscOptionsInt("-pow","","",user.pow,&user.pow,NULL)); 1012 1013 } 1014 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1015 1016 /* Create DMs for generator and network subsystems */ 1017 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen)); 1018 CHKERRQ(DMSetOptionsPrefix(user.dmgen,"dmgen_")); 1019 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet)); 1020 CHKERRQ(DMSetOptionsPrefix(user.dmnet,"dmnet_")); 1021 CHKERRQ(DMSetFromOptions(user.dmnet)); 1022 CHKERRQ(DMSetUp(user.dmnet)); 1023 /* Create a composite DM packer and add the two DMs */ 1024 CHKERRQ(DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid)); 1025 CHKERRQ(DMSetOptionsPrefix(user.dmpgrid,"pgrid_")); 1026 CHKERRQ(DMSetFromOptions(user.dmgen)); 1027 CHKERRQ(DMSetUp(user.dmgen)); 1028 CHKERRQ(DMCompositeAddDM(user.dmpgrid,user.dmgen)); 1029 CHKERRQ(DMCompositeAddDM(user.dmpgrid,user.dmnet)); 1030 1031 /* Read initial voltage vector and Ybus */ 1032 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview)); 1033 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview)); 1034 1035 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user.V0)); 1036 CHKERRQ(VecSetSizes(user.V0,PETSC_DECIDE,user.neqs_net)); 1037 CHKERRQ(VecLoad(user.V0,Xview)); 1038 1039 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Ybus)); 1040 CHKERRQ(MatSetSizes(user.Ybus,PETSC_DECIDE,PETSC_DECIDE,user.neqs_net,user.neqs_net)); 1041 CHKERRQ(MatSetType(user.Ybus,MATBAIJ)); 1042 /* CHKERRQ(MatSetBlockSize(ctx->Ybus,2)); */ 1043 CHKERRQ(MatLoad(user.Ybus,Ybusview)); 1044 1045 CHKERRQ(PetscViewerDestroy(&Xview)); 1046 CHKERRQ(PetscViewerDestroy(&Ybusview)); 1047 1048 /* Allocate space for Jacobians */ 1049 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.J)); 1050 CHKERRQ(MatSetSizes(user.J,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,user.neqs_pgrid)); 1051 CHKERRQ(MatSetFromOptions(user.J)); 1052 CHKERRQ(PreallocateJacobian(user.J,&user)); 1053 1054 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 1055 CHKERRQ(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,3)); 1056 CHKERRQ(MatSetFromOptions(user.Jacp)); 1057 CHKERRQ(MatSetUp(user.Jacp)); 1058 CHKERRQ(MatZeroEntries(user.Jacp)); /* initialize to zeros */ 1059 1060 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,&user.DRDP)); 1061 CHKERRQ(MatSetUp(user.DRDP)); 1062 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,1,NULL,&user.DRDU)); 1063 CHKERRQ(MatSetUp(user.DRDU)); 1064 1065 /* Create TAO solver and set desired solution method */ 1066 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 1067 CHKERRQ(TaoSetType(tao,TAOBLMVM)); 1068 /* 1069 Optimization starts 1070 */ 1071 /* Set initial solution guess */ 1072 CHKERRQ(VecCreateSeq(PETSC_COMM_WORLD,3,&p)); 1073 CHKERRQ(VecGetArray(p,&x_ptr)); 1074 x_ptr[0] = PG[0]; x_ptr[1] = PG[1]; x_ptr[2] = PG[2]; 1075 CHKERRQ(VecRestoreArray(p,&x_ptr)); 1076 1077 CHKERRQ(TaoSetSolution(tao,p)); 1078 /* Set routine for function and gradient evaluation */ 1079 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,&user)); 1080 1081 /* Set bounds for the optimization */ 1082 CHKERRQ(VecDuplicate(p,&lowerb)); 1083 CHKERRQ(VecDuplicate(p,&upperb)); 1084 CHKERRQ(VecGetArray(lowerb,&x_ptr)); 1085 x_ptr[0] = 0.5; x_ptr[1] = 0.5; x_ptr[2] = 0.5; 1086 CHKERRQ(VecRestoreArray(lowerb,&x_ptr)); 1087 CHKERRQ(VecGetArray(upperb,&x_ptr)); 1088 x_ptr[0] = 2.0; x_ptr[1] = 2.0; x_ptr[2] = 2.0; 1089 CHKERRQ(VecRestoreArray(upperb,&x_ptr)); 1090 CHKERRQ(TaoSetVariableBounds(tao,lowerb,upperb)); 1091 1092 /* Check for any TAO command line options */ 1093 CHKERRQ(TaoSetFromOptions(tao)); 1094 CHKERRQ(TaoGetKSP(tao,&ksp)); 1095 if (ksp) { 1096 CHKERRQ(KSPGetPC(ksp,&pc)); 1097 CHKERRQ(PCSetType(pc,PCNONE)); 1098 } 1099 1100 /* SOLVE THE APPLICATION */ 1101 CHKERRQ(TaoSolve(tao)); 1102 1103 CHKERRQ(VecView(p,PETSC_VIEWER_STDOUT_WORLD)); 1104 /* Free TAO data structures */ 1105 CHKERRQ(TaoDestroy(&tao)); 1106 1107 CHKERRQ(DMDestroy(&user.dmgen)); 1108 CHKERRQ(DMDestroy(&user.dmnet)); 1109 CHKERRQ(DMDestroy(&user.dmpgrid)); 1110 CHKERRQ(ISDestroy(&user.is_diff)); 1111 CHKERRQ(ISDestroy(&user.is_alg)); 1112 1113 CHKERRQ(MatDestroy(&user.J)); 1114 CHKERRQ(MatDestroy(&user.Jacp)); 1115 CHKERRQ(MatDestroy(&user.Ybus)); 1116 /* CHKERRQ(MatDestroy(&user.Sol)); */ 1117 CHKERRQ(VecDestroy(&user.V0)); 1118 CHKERRQ(VecDestroy(&p)); 1119 CHKERRQ(VecDestroy(&lowerb)); 1120 CHKERRQ(VecDestroy(&upperb)); 1121 CHKERRQ(MatDestroy(&user.DRDU)); 1122 CHKERRQ(MatDestroy(&user.DRDP)); 1123 ierr = PetscFinalize(); 1124 return ierr; 1125 } 1126 1127 /* ------------------------------------------------------------------ */ 1128 /* 1129 FormFunction - Evaluates the function and corresponding gradient. 1130 1131 Input Parameters: 1132 tao - the Tao context 1133 X - the input vector 1134 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 1135 1136 Output Parameters: 1137 f - the newly evaluated function 1138 G - the newly evaluated gradient 1139 */ 1140 PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx0) 1141 { 1142 TS ts,quadts; 1143 SNES snes_alg; 1144 Userctx *ctx = (Userctx*)ctx0; 1145 Vec X; 1146 PetscInt i; 1147 /* sensitivity context */ 1148 PetscScalar *x_ptr; 1149 Vec lambda[1],q; 1150 Vec mu[1]; 1151 PetscInt steps1,steps2,steps3; 1152 Vec DICDP[3]; 1153 Vec F_alg; 1154 PetscInt row_loc,col_loc; 1155 PetscScalar val; 1156 Vec Xdot; 1157 1158 PetscFunctionBegin; 1159 CHKERRQ(VecGetArrayRead(P,(const PetscScalar**)&x_ptr)); 1160 PG[0] = x_ptr[0]; 1161 PG[1] = x_ptr[1]; 1162 PG[2] = x_ptr[2]; 1163 CHKERRQ(VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr)); 1164 1165 ctx->stepnum = 0; 1166 1167 CHKERRQ(DMCreateGlobalVector(ctx->dmpgrid,&X)); 1168 1169 /* Create matrix to save solutions at each time step */ 1170 /* CHKERRQ(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */ 1171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1172 Create timestepping solver context 1173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1174 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 1175 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 1176 CHKERRQ(TSSetType(ts,TSCN)); 1177 CHKERRQ(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx)); 1178 CHKERRQ(TSSetIJacobian(ts,ctx->J,ctx->J,(TSIJacobian)IJacobian,ctx)); 1179 CHKERRQ(TSSetApplicationContext(ts,ctx)); 1180 /* Set RHS JacobianP */ 1181 CHKERRQ(TSSetRHSJacobianP(ts,ctx->Jacp,RHSJacobianP,ctx)); 1182 1183 CHKERRQ(TSCreateQuadratureTS(ts,PETSC_FALSE,&quadts)); 1184 CHKERRQ(TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,ctx)); 1185 CHKERRQ(TSSetRHSJacobian(quadts,ctx->DRDU,ctx->DRDU,(TSRHSJacobian)DRDUJacobianTranspose,ctx)); 1186 CHKERRQ(TSSetRHSJacobianP(quadts,ctx->DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,ctx)); 1187 1188 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1189 Set initial conditions 1190 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1191 CHKERRQ(SetInitialGuess(X,ctx)); 1192 1193 /* Approximate DICDP with finite difference, we want to zero out network variables */ 1194 for (i=0;i<3;i++) { 1195 CHKERRQ(VecDuplicate(X,&DICDP[i])); 1196 } 1197 CHKERRQ(DICDPFiniteDifference(X,DICDP,ctx)); 1198 1199 CHKERRQ(VecDuplicate(X,&F_alg)); 1200 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes_alg)); 1201 CHKERRQ(SNESSetFunction(snes_alg,F_alg,AlgFunction,ctx)); 1202 CHKERRQ(MatZeroEntries(ctx->J)); 1203 CHKERRQ(SNESSetJacobian(snes_alg,ctx->J,ctx->J,AlgJacobian,ctx)); 1204 CHKERRQ(SNESSetOptionsPrefix(snes_alg,"alg_")); 1205 CHKERRQ(SNESSetFromOptions(snes_alg)); 1206 ctx->alg_flg = PETSC_TRUE; 1207 /* Solve the algebraic equations */ 1208 CHKERRQ(SNESSolve(snes_alg,NULL,X)); 1209 1210 /* Just to set up the Jacobian structure */ 1211 CHKERRQ(VecDuplicate(X,&Xdot)); 1212 CHKERRQ(IJacobian(ts,0.0,X,Xdot,0.0,ctx->J,ctx->J,ctx)); 1213 CHKERRQ(VecDestroy(&Xdot)); 1214 1215 ctx->stepnum++; 1216 1217 /* 1218 Save trajectory of solution so that TSAdjointSolve() may be used 1219 */ 1220 CHKERRQ(TSSetSaveTrajectory(ts)); 1221 1222 CHKERRQ(TSSetTimeStep(ts,0.01)); 1223 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 1224 CHKERRQ(TSSetFromOptions(ts)); 1225 /* CHKERRQ(TSSetPostStep(ts,SaveSolution)); */ 1226 1227 /* Prefault period */ 1228 ctx->alg_flg = PETSC_FALSE; 1229 CHKERRQ(TSSetTime(ts,0.0)); 1230 CHKERRQ(TSSetMaxTime(ts,ctx->tfaulton)); 1231 CHKERRQ(TSSolve(ts,X)); 1232 CHKERRQ(TSGetStepNumber(ts,&steps1)); 1233 1234 /* Create the nonlinear solver for solving the algebraic system */ 1235 /* Note that although the algebraic system needs to be solved only for 1236 Idq and V, we reuse the entire system including xgen. The xgen 1237 variables are held constant by setting their residuals to 0 and 1238 putting a 1 on the Jacobian diagonal for xgen rows 1239 */ 1240 CHKERRQ(MatZeroEntries(ctx->J)); 1241 1242 /* Apply disturbance - resistive fault at ctx->faultbus */ 1243 /* This is done by adding shunt conductance to the diagonal location 1244 in the Ybus matrix */ 1245 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; /* Location for G */ 1246 val = 1/ctx->Rfault; 1247 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1248 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; /* Location for G */ 1249 val = 1/ctx->Rfault; 1250 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1251 1252 CHKERRQ(MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1253 CHKERRQ(MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1254 1255 ctx->alg_flg = PETSC_TRUE; 1256 /* Solve the algebraic equations */ 1257 CHKERRQ(SNESSolve(snes_alg,NULL,X)); 1258 1259 ctx->stepnum++; 1260 1261 /* Disturbance period */ 1262 ctx->alg_flg = PETSC_FALSE; 1263 CHKERRQ(TSSetTime(ts,ctx->tfaulton)); 1264 CHKERRQ(TSSetMaxTime(ts,ctx->tfaultoff)); 1265 CHKERRQ(TSSolve(ts,X)); 1266 CHKERRQ(TSGetStepNumber(ts,&steps2)); 1267 steps2 -= steps1; 1268 1269 /* Remove the fault */ 1270 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; 1271 val = -1/ctx->Rfault; 1272 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1273 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; 1274 val = -1/ctx->Rfault; 1275 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1276 1277 CHKERRQ(MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1278 CHKERRQ(MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1279 1280 CHKERRQ(MatZeroEntries(ctx->J)); 1281 1282 ctx->alg_flg = PETSC_TRUE; 1283 1284 /* Solve the algebraic equations */ 1285 CHKERRQ(SNESSolve(snes_alg,NULL,X)); 1286 1287 ctx->stepnum++; 1288 1289 /* Post-disturbance period */ 1290 ctx->alg_flg = PETSC_TRUE; 1291 CHKERRQ(TSSetTime(ts,ctx->tfaultoff)); 1292 CHKERRQ(TSSetMaxTime(ts,ctx->tmax)); 1293 CHKERRQ(TSSolve(ts,X)); 1294 CHKERRQ(TSGetStepNumber(ts,&steps3)); 1295 steps3 -= steps2; 1296 steps3 -= steps1; 1297 1298 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1299 Adjoint model starts here 1300 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1301 CHKERRQ(TSSetPostStep(ts,NULL)); 1302 CHKERRQ(MatCreateVecs(ctx->J,&lambda[0],NULL)); 1303 /* Set initial conditions for the adjoint integration */ 1304 CHKERRQ(VecZeroEntries(lambda[0])); 1305 1306 CHKERRQ(MatCreateVecs(ctx->Jacp,&mu[0],NULL)); 1307 CHKERRQ(VecZeroEntries(mu[0])); 1308 CHKERRQ(TSSetCostGradients(ts,1,lambda,mu)); 1309 1310 CHKERRQ(TSAdjointSetSteps(ts,steps3)); 1311 CHKERRQ(TSAdjointSolve(ts)); 1312 1313 CHKERRQ(MatZeroEntries(ctx->J)); 1314 /* Applying disturbance - resistive fault at ctx->faultbus */ 1315 /* This is done by deducting shunt conductance to the diagonal location 1316 in the Ybus matrix */ 1317 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; /* Location for G */ 1318 val = 1./ctx->Rfault; 1319 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1320 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; /* Location for G */ 1321 val = 1./ctx->Rfault; 1322 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1323 1324 CHKERRQ(MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1325 CHKERRQ(MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1326 1327 /* Set number of steps for the adjoint integration */ 1328 CHKERRQ(TSAdjointSetSteps(ts,steps2)); 1329 CHKERRQ(TSAdjointSolve(ts)); 1330 1331 CHKERRQ(MatZeroEntries(ctx->J)); 1332 /* remove the fault */ 1333 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; /* Location for G */ 1334 val = -1./ctx->Rfault; 1335 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1336 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; /* Location for G */ 1337 val = -1./ctx->Rfault; 1338 CHKERRQ(MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 1339 1340 CHKERRQ(MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1341 CHKERRQ(MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY)); 1342 1343 /* Set number of steps for the adjoint integration */ 1344 CHKERRQ(TSAdjointSetSteps(ts,steps1)); 1345 CHKERRQ(TSAdjointSolve(ts)); 1346 1347 CHKERRQ(ComputeSensiP(lambda[0],mu[0],DICDP,ctx)); 1348 CHKERRQ(VecCopy(mu[0],G)); 1349 1350 CHKERRQ(TSGetQuadratureTS(ts,NULL,&quadts)); 1351 CHKERRQ(TSGetSolution(quadts,&q)); 1352 CHKERRQ(VecGetArray(q,&x_ptr)); 1353 *f = x_ptr[0]; 1354 x_ptr[0] = 0; 1355 CHKERRQ(VecRestoreArray(q,&x_ptr)); 1356 1357 CHKERRQ(VecDestroy(&lambda[0])); 1358 CHKERRQ(VecDestroy(&mu[0])); 1359 1360 CHKERRQ(SNESDestroy(&snes_alg)); 1361 CHKERRQ(VecDestroy(&F_alg)); 1362 CHKERRQ(VecDestroy(&X)); 1363 CHKERRQ(TSDestroy(&ts)); 1364 for (i=0;i<3;i++) { 1365 CHKERRQ(VecDestroy(&DICDP[i])); 1366 } 1367 PetscFunctionReturn(0); 1368 } 1369 1370 /*TEST 1371 1372 build: 1373 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 1374 1375 test: 1376 args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2 1377 localrunfiles: petscoptions X.bin Ybus.bin 1378 1379 TEST*/ 1380