1 static char help[] = "Using finite difference for the problem in ex9busopt.c \n\n"; 2 3 /* 4 Use finite difference approximations to solve the same optimization problem as in ex9busopt.c. 5 */ 6 7 #include <petsctao.h> 8 #include <petscts.h> 9 #include <petscdm.h> 10 #include <petscdmda.h> 11 #include <petscdmcomposite.h> 12 13 PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*); 14 15 #define freq 60 16 #define w_s (2*PETSC_PI*freq) 17 18 /* Sizes and indices */ 19 const PetscInt nbus = 9; /* Number of network buses */ 20 const PetscInt ngen = 3; /* Number of generators */ 21 const PetscInt nload = 3; /* Number of loads */ 22 const PetscInt gbus[3] = {0,1,2}; /* Buses at which generators are incident */ 23 const PetscInt lbus[3] = {4,5,7}; /* Buses at which loads are incident */ 24 25 /* Generator real and reactive powers (found via loadflow) */ 26 PetscScalar PG[3] = { 0.69,1.59,0.69}; 27 /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/ 28 const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588}; 29 /* Generator constants */ 30 const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */ 31 const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */ 32 const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */ 33 const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */ 34 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 */ 35 const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */ 36 const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */ 37 const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */ 38 PetscScalar M[3]; /* M = 2*H/w_s */ 39 PetscScalar D[3]; /* D = 0.1*M */ 40 41 PetscScalar TM[3]; /* Mechanical Torque */ 42 /* Exciter system constants */ 43 const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */ 44 const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */ 45 const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */ 46 const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */ 47 const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */ 48 const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */ 49 const PetscScalar k1[3] = {0.0039,0.0039,0.0039}; 50 const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */ 51 52 PetscScalar Vref[3]; 53 /* Load constants 54 We use a composite load model that describes the load and reactive powers at each time instant as follows 55 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i 56 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i 57 where 58 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads 59 ld_alphap,ld_alphap - Percentage contribution (weights) or loads 60 P_D0 - Real power load 61 Q_D0 - Reactive power load 62 V_m(t) - Voltage magnitude at time t 63 V_m0 - Voltage magnitude at t = 0 64 ld_betap, ld_betaq - exponents describing the load model for real and reactive part 65 66 Note: All loads have the same characteristic currently. 67 */ 68 const PetscScalar PD0[3] = {1.25,0.9,1.0}; 69 const PetscScalar QD0[3] = {0.5,0.3,0.35}; 70 const PetscInt ld_nsegsp[3] = {3,3,3}; 71 const PetscScalar ld_alphap[3] = {1.0,0.0,0.0}; 72 const PetscScalar ld_betap[3] = {2.0,1.0,0.0}; 73 const PetscInt ld_nsegsq[3] = {3,3,3}; 74 const PetscScalar ld_alphaq[3] = {1.0,0.0,0.0}; 75 const PetscScalar ld_betaq[3] = {2.0,1.0,0.0}; 76 77 typedef struct { 78 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */ 79 DM dmpgrid; /* Composite DM to manage the entire power grid */ 80 Mat Ybus; /* Network admittance matrix */ 81 Vec V0; /* Initial voltage vector (Power flow solution) */ 82 PetscReal tfaulton,tfaultoff; /* Fault on and off times */ 83 PetscInt faultbus; /* Fault bus */ 84 PetscScalar Rfault; 85 PetscReal t0,tmax; 86 PetscInt neqs_gen,neqs_net,neqs_pgrid; 87 Mat Sol; /* Matrix to save solution at each time step */ 88 PetscInt stepnum; 89 PetscBool alg_flg; 90 PetscReal t; 91 IS is_diff; /* indices for differential equations */ 92 IS is_alg; /* indices for algebraic equations */ 93 PetscReal freq_u,freq_l; /* upper and lower frequency limit */ 94 PetscInt pow; /* power coefficient used in the cost function */ 95 Vec vec_q; 96 } Userctx; 97 98 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 99 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr, PetscScalar *Fi) 100 { 101 PetscFunctionBegin; 102 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 103 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 104 PetscFunctionReturn(0); 105 } 106 107 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 108 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd, PetscScalar *Fq) 109 { 110 PetscFunctionBegin; 111 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 112 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode SetInitialGuess(Vec X,Userctx *user) 117 { 118 PetscErrorCode ierr; 119 Vec Xgen,Xnet; 120 PetscScalar *xgen,*xnet; 121 PetscInt i,idx=0; 122 PetscScalar Vr,Vi,IGr,IGi,Vm,Vm2; 123 PetscScalar Eqp,Edp,delta; 124 PetscScalar Efd,RF,VR; /* Exciter variables */ 125 PetscScalar Id,Iq; /* Generator dq axis currents */ 126 PetscScalar theta,Vd,Vq,SE; 127 128 PetscFunctionBegin; 129 M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s; 130 D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2]; 131 132 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 133 134 /* Network subsystem initialization */ 135 ierr = VecCopy(user->V0,Xnet);CHKERRQ(ierr); 136 137 /* Generator subsystem initialization */ 138 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 139 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 140 141 for (i=0; i < ngen; i++) { 142 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 143 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 144 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 145 IGr = (Vr*PG[i] + Vi*QG[i])/Vm2; 146 IGi = (Vi*PG[i] - Vr*QG[i])/Vm2; 147 148 delta = PetscAtan2Real(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */ 149 150 theta = PETSC_PI/2.0 - delta; 151 152 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */ 153 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */ 154 155 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 156 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 157 158 Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */ 159 Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */ 160 161 TM[i] = PG[i]; 162 163 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */ 164 xgen[idx] = Eqp; 165 xgen[idx+1] = Edp; 166 xgen[idx+2] = delta; 167 xgen[idx+3] = w_s; 168 169 idx = idx + 4; 170 171 xgen[idx] = Id; 172 xgen[idx+1] = Iq; 173 174 idx = idx + 2; 175 176 /* Exciter */ 177 Efd = Eqp + (Xd[i] - Xdp[i])*Id; 178 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 179 VR = KE[i]*Efd + SE; 180 RF = KF[i]*Efd/TF[i]; 181 182 xgen[idx] = Efd; 183 xgen[idx+1] = RF; 184 xgen[idx+2] = VR; 185 186 Vref[i] = Vm + (VR/KA[i]); 187 188 idx = idx + 3; 189 } 190 191 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 192 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 193 194 /* ierr = VecView(Xgen,0);CHKERRQ(ierr); */ 195 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet);CHKERRQ(ierr); 196 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 /* Computes F = [-f(x,y);g(x,y)] */ 201 PetscErrorCode ResidualFunction(SNES snes,Vec X, Vec F, Userctx *user) 202 { 203 PetscErrorCode ierr; 204 Vec Xgen,Xnet,Fgen,Fnet; 205 PetscScalar *xgen,*xnet,*fgen,*fnet; 206 PetscInt i,idx=0; 207 PetscScalar Vr,Vi,Vm,Vm2; 208 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 209 PetscScalar Efd,RF,VR; /* Exciter variables */ 210 PetscScalar Id,Iq; /* Generator dq axis currents */ 211 PetscScalar Vd,Vq,SE; 212 PetscScalar IGr,IGi,IDr,IDi; 213 PetscScalar Zdq_inv[4],det; 214 PetscScalar PD,QD,Vm0,*v0; 215 PetscInt k; 216 217 PetscFunctionBegin; 218 ierr = VecZeroEntries(F);CHKERRQ(ierr); 219 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 220 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 221 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 222 ierr = DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet);CHKERRQ(ierr); 223 224 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 225 The generator current injection, IG, and load current injection, ID are added later 226 */ 227 /* Note that the values in Ybus are stored assuming the imaginary current balance 228 equation is ordered first followed by real current balance equation for each bus. 229 Thus imaginary current contribution goes in location 2*i, and 230 real current contribution in 2*i+1 231 */ 232 ierr = MatMult(user->Ybus,Xnet,Fnet);CHKERRQ(ierr); 233 234 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 235 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 236 ierr = VecGetArray(Fgen,&fgen);CHKERRQ(ierr); 237 ierr = VecGetArray(Fnet,&fnet);CHKERRQ(ierr); 238 239 /* Generator subsystem */ 240 for (i=0; i < ngen; i++) { 241 Eqp = xgen[idx]; 242 Edp = xgen[idx+1]; 243 delta = xgen[idx+2]; 244 w = xgen[idx+3]; 245 Id = xgen[idx+4]; 246 Iq = xgen[idx+5]; 247 Efd = xgen[idx+6]; 248 RF = xgen[idx+7]; 249 VR = xgen[idx+8]; 250 251 /* Generator differential equations */ 252 fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; 253 fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; 254 fgen[idx+2] = -w + w_s; 255 fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; 256 257 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 258 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 259 260 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 261 /* Algebraic equations for stator currents */ 262 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 263 264 Zdq_inv[0] = Rs[i]/det; 265 Zdq_inv[1] = Xqp[i]/det; 266 Zdq_inv[2] = -Xdp[i]/det; 267 Zdq_inv[3] = Rs[i]/det; 268 269 fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 270 fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 271 272 /* Add generator current injection to network */ 273 ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr); 274 275 fnet[2*gbus[i]] -= IGi; 276 fnet[2*gbus[i]+1] -= IGr; 277 278 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 279 280 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 281 282 /* Exciter differential equations */ 283 fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; 284 fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; 285 fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 286 287 idx = idx + 9; 288 } 289 290 ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr); 291 for (i=0; i < nload; i++) { 292 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 293 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 294 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 295 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 296 PD = QD = 0.0; 297 for (k=0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 298 for (k=0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 299 300 /* Load currents */ 301 IDr = (PD*Vr + QD*Vi)/Vm2; 302 IDi = (-QD*Vr + PD*Vi)/Vm2; 303 304 fnet[2*lbus[i]] += IDi; 305 fnet[2*lbus[i]+1] += IDr; 306 } 307 ierr = VecRestoreArray(user->V0,&v0);CHKERRQ(ierr); 308 309 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 310 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 311 ierr = VecRestoreArray(Fgen,&fgen);CHKERRQ(ierr); 312 ierr = VecRestoreArray(Fnet,&fnet);CHKERRQ(ierr); 313 314 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet);CHKERRQ(ierr); 315 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 316 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 317 PetscFunctionReturn(0); 318 } 319 320 /* \dot{x} - f(x,y) 321 g(x,y) = 0 322 */ 323 PetscErrorCode IFunction(TS ts,PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user) 324 { 325 PetscErrorCode ierr; 326 SNES snes; 327 PetscScalar *f; 328 const PetscScalar *xdot; 329 PetscInt i; 330 331 PetscFunctionBegin; 332 user->t = t; 333 334 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 335 ierr = ResidualFunction(snes,X,F,user);CHKERRQ(ierr); 336 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 337 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 338 for (i=0;i < ngen;i++) { 339 f[9*i] += xdot[9*i]; 340 f[9*i+1] += xdot[9*i+1]; 341 f[9*i+2] += xdot[9*i+2]; 342 f[9*i+3] += xdot[9*i+3]; 343 f[9*i+6] += xdot[9*i+6]; 344 f[9*i+7] += xdot[9*i+7]; 345 f[9*i+8] += xdot[9*i+8]; 346 } 347 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 348 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351 352 /* This function is used for solving the algebraic system only during fault on and 353 off times. It computes the entire F and then zeros out the part corresponding to 354 differential equations 355 F = [0;g(y)]; 356 */ 357 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx) 358 { 359 PetscErrorCode ierr; 360 Userctx *user=(Userctx*)ctx; 361 PetscScalar *f; 362 PetscInt i; 363 364 PetscFunctionBegin; 365 ierr = ResidualFunction(snes,X,F,user);CHKERRQ(ierr); 366 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 367 for (i=0; i < ngen; i++) { 368 f[9*i] = 0; 369 f[9*i+1] = 0; 370 f[9*i+2] = 0; 371 f[9*i+3] = 0; 372 f[9*i+6] = 0; 373 f[9*i+7] = 0; 374 f[9*i+8] = 0; 375 } 376 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user) 381 { 382 PetscErrorCode ierr; 383 PetscInt *d_nnz; 384 PetscInt i,idx=0,start=0; 385 PetscInt ncols; 386 387 PetscFunctionBegin; 388 ierr = PetscMalloc1(user->neqs_pgrid,&d_nnz);CHKERRQ(ierr); 389 for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0; 390 /* Generator subsystem */ 391 for (i=0; i < ngen; i++) { 392 393 d_nnz[idx] += 3; 394 d_nnz[idx+1] += 2; 395 d_nnz[idx+2] += 2; 396 d_nnz[idx+3] += 5; 397 d_nnz[idx+4] += 6; 398 d_nnz[idx+5] += 6; 399 400 d_nnz[user->neqs_gen+2*gbus[i]] += 3; 401 d_nnz[user->neqs_gen+2*gbus[i]+1] += 3; 402 403 d_nnz[idx+6] += 2; 404 d_nnz[idx+7] += 2; 405 d_nnz[idx+8] += 5; 406 407 idx = idx + 9; 408 } 409 410 start = user->neqs_gen; 411 412 for (i=0; i < nbus; i++) { 413 ierr = MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 414 d_nnz[start+2*i] += ncols; 415 d_nnz[start+2*i+1] += ncols; 416 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 417 } 418 419 ierr = MatSeqAIJSetPreallocation(J,0,d_nnz);CHKERRQ(ierr); 420 421 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 /* 426 J = [-df_dx, -df_dy 427 dg_dx, dg_dy] 428 */ 429 PetscErrorCode ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx) 430 { 431 PetscErrorCode ierr; 432 Userctx *user=(Userctx*)ctx; 433 Vec Xgen,Xnet; 434 PetscScalar *xgen,*xnet; 435 PetscInt i,idx=0; 436 PetscScalar Vr,Vi,Vm,Vm2; 437 PetscScalar Eqp,Edp,delta; /* Generator variables */ 438 PetscScalar Efd; /* Exciter variables */ 439 PetscScalar Id,Iq; /* Generator dq axis currents */ 440 PetscScalar Vd,Vq; 441 PetscScalar val[10]; 442 PetscInt row[2],col[10]; 443 PetscInt net_start=user->neqs_gen; 444 PetscScalar Zdq_inv[4],det; 445 PetscScalar dVd_dVr,dVd_dVi,dVq_dVr,dVq_dVi,dVd_ddelta,dVq_ddelta; 446 PetscScalar dIGr_ddelta,dIGi_ddelta,dIGr_dId,dIGr_dIq,dIGi_dId,dIGi_dIq; 447 PetscScalar dSE_dEfd; 448 PetscScalar dVm_dVd,dVm_dVq,dVm_dVr,dVm_dVi; 449 PetscInt ncols; 450 const PetscInt *cols; 451 const PetscScalar *yvals; 452 PetscInt k; 453 PetscScalar PD,QD,Vm0,*v0,Vm4; 454 PetscScalar dPD_dVr,dPD_dVi,dQD_dVr,dQD_dVi; 455 PetscScalar dIDr_dVr,dIDr_dVi,dIDi_dVr,dIDi_dVi; 456 457 PetscFunctionBegin; 458 ierr = MatZeroEntries(B);CHKERRQ(ierr); 459 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 460 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 461 462 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 463 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 464 465 /* Generator subsystem */ 466 for (i=0; i < ngen; i++) { 467 Eqp = xgen[idx]; 468 Edp = xgen[idx+1]; 469 delta = xgen[idx+2]; 470 Id = xgen[idx+4]; 471 Iq = xgen[idx+5]; 472 Efd = xgen[idx+6]; 473 474 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */ 475 row[0] = idx; 476 col[0] = idx; col[1] = idx+4; col[2] = idx+6; 477 val[0] = 1/ Td0p[i]; val[1] = (Xd[i] - Xdp[i])/ Td0p[i]; val[2] = -1/Td0p[i]; 478 479 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 480 481 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */ 482 row[0] = idx + 1; 483 col[0] = idx + 1; col[1] = idx+5; 484 val[0] = 1/Tq0p[i]; val[1] = -(Xq[i] - Xqp[i])/Tq0p[i]; 485 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 486 487 /* fgen[idx+2] = - w + w_s; */ 488 row[0] = idx + 2; 489 col[0] = idx + 2; col[1] = idx + 3; 490 val[0] = 0; val[1] = -1; 491 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 492 493 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */ 494 row[0] = idx + 3; 495 col[0] = idx; col[1] = idx + 1; col[2] = idx + 3; col[3] = idx + 4; col[4] = idx + 5; 496 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]; 497 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 498 499 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 500 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 501 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 502 503 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 504 505 Zdq_inv[0] = Rs[i]/det; 506 Zdq_inv[1] = Xqp[i]/det; 507 Zdq_inv[2] = -Xdp[i]/det; 508 Zdq_inv[3] = Rs[i]/det; 509 510 dVd_dVr = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta); 511 dVq_dVr = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta); 512 dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta); 513 dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta); 514 515 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */ 516 row[0] = idx+4; 517 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 518 val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0]; val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta; 519 col[3] = idx + 4; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 520 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; 521 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 522 523 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */ 524 row[0] = idx+5; 525 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 526 val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2]; val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta; 527 col[3] = idx + 5; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 528 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; 529 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 530 531 dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta); 532 dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta); 533 dIGr_dId = PetscSinScalar(delta); dIGr_dIq = PetscCosScalar(delta); 534 dIGi_dId = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta); 535 536 /* fnet[2*gbus[i]] -= IGi; */ 537 row[0] = net_start + 2*gbus[i]; 538 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 539 val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq; 540 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 541 542 /* fnet[2*gbus[i]+1] -= IGr; */ 543 row[0] = net_start + 2*gbus[i]+1; 544 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 545 val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq; 546 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 547 548 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 549 550 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */ 551 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */ 552 553 dSE_dEfd = k1[i]*k2[i]*PetscExpScalar(k2[i]*Efd); 554 555 row[0] = idx + 6; 556 col[0] = idx + 6; col[1] = idx + 8; 557 val[0] = (KE[i] + dSE_dEfd)/TE[i]; val[1] = -1/TE[i]; 558 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 559 560 /* Exciter differential equations */ 561 562 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */ 563 row[0] = idx + 7; 564 col[0] = idx + 6; col[1] = idx + 7; 565 val[0] = (-KF[i]/TF[i])/TF[i]; val[1] = 1/TF[i]; 566 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 567 568 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */ 569 /* Vm = (Vd^2 + Vq^2)^0.5; */ 570 dVm_dVd = Vd/Vm; dVm_dVq = Vq/Vm; 571 dVm_dVr = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr; 572 dVm_dVi = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi; 573 row[0] = idx + 8; 574 col[0] = idx + 6; col[1] = idx + 7; col[2] = idx + 8; 575 val[0] = (KA[i]*KF[i]/TF[i])/TA[i]; val[1] = -KA[i]/TA[i]; val[2] = 1/TA[i]; 576 col[3] = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1; 577 val[3] = KA[i]*dVm_dVr/TA[i]; val[4] = KA[i]*dVm_dVi/TA[i]; 578 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 579 idx = idx + 9; 580 } 581 582 for (i=0; i<nbus; i++) { 583 ierr = MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 584 row[0] = net_start + 2*i; 585 for (k=0; k<ncols; k++) { 586 col[k] = net_start + cols[k]; 587 val[k] = yvals[k]; 588 } 589 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 590 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 591 592 ierr = MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 593 row[0] = net_start + 2*i+1; 594 for (k=0; k<ncols; k++) { 595 col[k] = net_start + cols[k]; 596 val[k] = yvals[k]; 597 } 598 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 599 ierr = MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 600 } 601 602 ierr = MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 603 ierr = MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 604 605 ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr); 606 for (i=0; i < nload; i++) { 607 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 608 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 609 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; Vm4 = Vm2*Vm2; 610 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 611 PD = QD = 0.0; 612 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0; 613 for (k=0; k < ld_nsegsp[i]; k++) { 614 PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 615 dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2)); 616 dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2)); 617 } 618 for (k=0; k < ld_nsegsq[i]; k++) { 619 QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 620 dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2)); 621 dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2)); 622 } 623 624 /* IDr = (PD*Vr + QD*Vi)/Vm2; */ 625 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */ 626 627 dIDr_dVr = (dPD_dVr*Vr + dQD_dVr*Vi + PD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vr)/Vm4; 628 dIDr_dVi = (dPD_dVi*Vr + dQD_dVi*Vi + QD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vi)/Vm4; 629 630 dIDi_dVr = (-dQD_dVr*Vr + dPD_dVr*Vi - QD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vr)/Vm4; 631 dIDi_dVi = (-dQD_dVi*Vr + dPD_dVi*Vi + PD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vi)/Vm4; 632 633 /* fnet[2*lbus[i]] += IDi; */ 634 row[0] = net_start + 2*lbus[i]; 635 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 636 val[0] = dIDi_dVr; val[1] = dIDi_dVi; 637 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 638 /* fnet[2*lbus[i]+1] += IDr; */ 639 row[0] = net_start + 2*lbus[i]+1; 640 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 641 val[0] = dIDr_dVr; val[1] = dIDr_dVi; 642 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 643 } 644 ierr = VecRestoreArray(user->V0,&v0);CHKERRQ(ierr); 645 646 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 647 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 648 649 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 650 651 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 652 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 /* 657 J = [I, 0 658 dg_dx, dg_dy] 659 */ 660 PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx) 661 { 662 PetscErrorCode ierr; 663 Userctx *user=(Userctx*)ctx; 664 665 PetscFunctionBegin; 666 ierr = ResidualJacobian(snes,X,A,B,ctx);CHKERRQ(ierr); 667 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);CHKERRQ(ierr); 668 ierr = MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL);CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 /* 673 J = [a*I-df_dx, -df_dy 674 dg_dx, dg_dy] 675 */ 676 677 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user) 678 { 679 PetscErrorCode ierr; 680 SNES snes; 681 PetscScalar atmp = (PetscScalar) a; 682 PetscInt i,row; 683 684 PetscFunctionBegin; 685 user->t = t; 686 687 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 688 ierr = ResidualJacobian(snes,X,A,B,user);CHKERRQ(ierr); 689 for (i=0;i < ngen;i++) { 690 row = 9*i; 691 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 692 row = 9*i+1; 693 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 694 row = 9*i+2; 695 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 696 row = 9*i+3; 697 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 698 row = 9*i+6; 699 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 700 row = 9*i+7; 701 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 702 row = 9*i+8; 703 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 704 } 705 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 706 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 710 static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,Userctx *user) 711 { 712 PetscErrorCode ierr; 713 PetscScalar *r; 714 const PetscScalar *u; 715 PetscInt idx; 716 Vec Xgen,Xnet; 717 PetscScalar *xgen; 718 PetscInt i; 719 720 PetscFunctionBegin; 721 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 722 ierr = DMCompositeScatter(user->dmpgrid,U,Xgen,Xnet);CHKERRQ(ierr); 723 724 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 725 726 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 727 ierr = VecGetArray(R,&r);CHKERRQ(ierr); 728 r[0] = 0.; 729 730 idx = 0; 731 for (i=0;i<ngen;i++) { 732 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); 733 idx += 9; 734 } 735 ierr = VecRestoreArray(R,&r);CHKERRQ(ierr); 736 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 737 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 738 PetscFunctionReturn(0); 739 } 740 741 static PetscErrorCode MonitorUpdateQ(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx0) 742 { 743 PetscErrorCode ierr; 744 Vec C,*Y; 745 PetscInt Nr; 746 PetscReal h,theta; 747 Userctx *ctx=(Userctx*)ctx0; 748 749 PetscFunctionBegin; 750 theta = 0.5; 751 ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr); 752 ierr = TSGetTimeStep(ts,&h);CHKERRQ(ierr); 753 ierr = VecDuplicate(ctx->vec_q,&C);CHKERRQ(ierr); 754 /* compute integrals */ 755 if (stepnum>0) { 756 ierr = CostIntegrand(ts,time,X,C,ctx);CHKERRQ(ierr); 757 ierr = VecAXPY(ctx->vec_q,h*theta,C);CHKERRQ(ierr); 758 ierr = CostIntegrand(ts,time+h*theta,Y[0],C,ctx);CHKERRQ(ierr); 759 ierr = VecAXPY(ctx->vec_q,h*(1-theta),C);CHKERRQ(ierr); 760 } 761 ierr = VecDestroy(&C);CHKERRQ(ierr); 762 PetscFunctionReturn(0); 763 } 764 765 int main(int argc,char **argv) 766 { 767 Userctx user; 768 Vec p; 769 PetscScalar *x_ptr; 770 PetscErrorCode ierr; 771 PetscMPIInt size; 772 PetscInt i; 773 KSP ksp; 774 PC pc; 775 PetscInt *idx2; 776 Tao tao; 777 Vec lowerb,upperb; 778 779 PetscFunctionBeginUser; 780 ierr = PetscInitialize(&argc,&argv,"petscoptions",help);if (ierr) return ierr; 781 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 782 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 783 784 ierr = VecCreateSeq(PETSC_COMM_WORLD,1,&user.vec_q);CHKERRQ(ierr); 785 786 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 787 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 788 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 789 790 /* Create indices for differential and algebraic equations */ 791 ierr = PetscMalloc1(7*ngen,&idx2);CHKERRQ(ierr); 792 for (i=0; i<ngen; i++) { 793 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; 794 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 795 } 796 ierr = ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);CHKERRQ(ierr); 797 ierr = ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);CHKERRQ(ierr); 798 ierr = PetscFree(idx2);CHKERRQ(ierr); 799 800 /* Set run time options */ 801 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 802 { 803 user.tfaulton = 1.0; 804 user.tfaultoff = 1.2; 805 user.Rfault = 0.0001; 806 user.faultbus = 8; 807 ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr); 808 ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr); 809 ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr); 810 user.t0 = 0.0; 811 user.tmax = 1.5; 812 ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr); 813 ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr); 814 user.freq_u = 61.0; 815 user.freq_l = 59.0; 816 user.pow = 2; 817 ierr = PetscOptionsReal("-frequ","","",user.freq_u,&user.freq_u,NULL);CHKERRQ(ierr); 818 ierr = PetscOptionsReal("-freql","","",user.freq_l,&user.freq_l,NULL);CHKERRQ(ierr); 819 ierr = PetscOptionsInt("-pow","","",user.pow,&user.pow,NULL);CHKERRQ(ierr); 820 821 } 822 ierr = PetscOptionsEnd();CHKERRQ(ierr); 823 824 /* Create DMs for generator and network subsystems */ 825 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); 826 ierr = DMSetOptionsPrefix(user.dmgen,"dmgen_");CHKERRQ(ierr); 827 ierr = DMSetFromOptions(user.dmgen);CHKERRQ(ierr); 828 ierr = DMSetUp(user.dmgen);CHKERRQ(ierr); 829 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); 830 ierr = DMSetOptionsPrefix(user.dmnet,"dmnet_");CHKERRQ(ierr); 831 ierr = DMSetFromOptions(user.dmnet);CHKERRQ(ierr); 832 ierr = DMSetUp(user.dmnet);CHKERRQ(ierr); 833 /* Create a composite DM packer and add the two DMs */ 834 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);CHKERRQ(ierr); 835 ierr = DMSetOptionsPrefix(user.dmpgrid,"pgrid_");CHKERRQ(ierr); 836 ierr = DMCompositeAddDM(user.dmpgrid,user.dmgen);CHKERRQ(ierr); 837 ierr = DMCompositeAddDM(user.dmpgrid,user.dmnet);CHKERRQ(ierr); 838 839 /* Create TAO solver and set desired solution method */ 840 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 841 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr); 842 /* 843 Optimization starts 844 */ 845 /* Set initial solution guess */ 846 ierr = VecCreateSeq(PETSC_COMM_WORLD,3,&p);CHKERRQ(ierr); 847 ierr = VecGetArray(p,&x_ptr);CHKERRQ(ierr); 848 x_ptr[0] = PG[0]; x_ptr[1] = PG[1]; x_ptr[2] = PG[2]; 849 ierr = VecRestoreArray(p,&x_ptr);CHKERRQ(ierr); 850 851 ierr = TaoSetInitialVector(tao,p);CHKERRQ(ierr); 852 /* Set routine for function and gradient evaluation */ 853 ierr = TaoSetObjectiveRoutine(tao,FormFunction,(void *)&user);CHKERRQ(ierr); 854 ierr = TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&user);CHKERRQ(ierr); 855 856 /* Set bounds for the optimization */ 857 ierr = VecDuplicate(p,&lowerb);CHKERRQ(ierr); 858 ierr = VecDuplicate(p,&upperb);CHKERRQ(ierr); 859 ierr = VecGetArray(lowerb,&x_ptr);CHKERRQ(ierr); 860 x_ptr[0] = 0.5; x_ptr[1] = 0.5; x_ptr[2] = 0.5; 861 ierr = VecRestoreArray(lowerb,&x_ptr);CHKERRQ(ierr); 862 ierr = VecGetArray(upperb,&x_ptr);CHKERRQ(ierr); 863 x_ptr[0] = 2.0; x_ptr[1] = 2.0; x_ptr[2] = 2.0; 864 ierr = VecRestoreArray(upperb,&x_ptr);CHKERRQ(ierr); 865 ierr = TaoSetVariableBounds(tao,lowerb,upperb);CHKERRQ(ierr); 866 867 /* Check for any TAO command line options */ 868 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 869 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 870 if (ksp) { 871 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 872 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 873 } 874 875 /* SOLVE THE APPLICATION */ 876 ierr = TaoSolve(tao);CHKERRQ(ierr); 877 878 ierr = VecView(p,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 879 /* Free TAO data structures */ 880 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 881 ierr = VecDestroy(&user.vec_q);CHKERRQ(ierr); 882 ierr = VecDestroy(&lowerb);CHKERRQ(ierr); 883 ierr = VecDestroy(&upperb);CHKERRQ(ierr); 884 ierr = VecDestroy(&p);CHKERRQ(ierr); 885 ierr = DMDestroy(&user.dmgen);CHKERRQ(ierr); 886 ierr = DMDestroy(&user.dmnet);CHKERRQ(ierr); 887 ierr = DMDestroy(&user.dmpgrid);CHKERRQ(ierr); 888 ierr = ISDestroy(&user.is_diff);CHKERRQ(ierr); 889 ierr = ISDestroy(&user.is_alg);CHKERRQ(ierr); 890 ierr = PetscFinalize(); 891 return ierr; 892 } 893 894 /* ------------------------------------------------------------------ */ 895 /* 896 FormFunction - Evaluates the function and corresponding gradient. 897 898 Input Parameters: 899 tao - the Tao context 900 X - the input vector 901 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine() 902 903 Output Parameters: 904 f - the newly evaluated function 905 */ 906 PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0) 907 { 908 TS ts; 909 SNES snes_alg; 910 PetscErrorCode ierr; 911 Userctx *ctx = (Userctx*)ctx0; 912 Vec X; 913 Mat J; 914 /* sensitivity context */ 915 PetscScalar *x_ptr; 916 PetscViewer Xview,Ybusview; 917 Vec F_alg; 918 Vec Xdot; 919 PetscInt row_loc,col_loc; 920 PetscScalar val; 921 922 ierr = VecGetArrayRead(P,(const PetscScalar**)&x_ptr);CHKERRQ(ierr); 923 PG[0] = x_ptr[0]; 924 PG[1] = x_ptr[1]; 925 PG[2] = x_ptr[2]; 926 ierr = VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);CHKERRQ(ierr); 927 928 ctx->stepnum = 0; 929 930 ierr = VecZeroEntries(ctx->vec_q);CHKERRQ(ierr); 931 932 /* Read initial voltage vector and Ybus */ 933 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview);CHKERRQ(ierr); 934 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview);CHKERRQ(ierr); 935 936 ierr = VecCreate(PETSC_COMM_WORLD,&ctx->V0);CHKERRQ(ierr); 937 ierr = VecSetSizes(ctx->V0,PETSC_DECIDE,ctx->neqs_net);CHKERRQ(ierr); 938 ierr = VecLoad(ctx->V0,Xview);CHKERRQ(ierr); 939 940 ierr = MatCreate(PETSC_COMM_WORLD,&ctx->Ybus);CHKERRQ(ierr); 941 ierr = MatSetSizes(ctx->Ybus,PETSC_DECIDE,PETSC_DECIDE,ctx->neqs_net,ctx->neqs_net);CHKERRQ(ierr); 942 ierr = MatSetType(ctx->Ybus,MATBAIJ);CHKERRQ(ierr); 943 /* ierr = MatSetBlockSize(ctx->Ybus,2);CHKERRQ(ierr); */ 944 ierr = MatLoad(ctx->Ybus,Ybusview);CHKERRQ(ierr); 945 946 ierr = PetscViewerDestroy(&Xview);CHKERRQ(ierr); 947 ierr = PetscViewerDestroy(&Ybusview);CHKERRQ(ierr); 948 949 ierr = DMCreateGlobalVector(ctx->dmpgrid,&X);CHKERRQ(ierr); 950 951 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 952 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ctx->neqs_pgrid,ctx->neqs_pgrid);CHKERRQ(ierr); 953 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 954 ierr = PreallocateJacobian(J,ctx);CHKERRQ(ierr); 955 956 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 957 Create timestepping solver context 958 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 959 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 960 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 961 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 962 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);CHKERRQ(ierr); 963 ierr = TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,ctx);CHKERRQ(ierr); 964 ierr = TSSetApplicationContext(ts,ctx);CHKERRQ(ierr); 965 966 ierr = TSMonitorSet(ts,MonitorUpdateQ,ctx,NULL);CHKERRQ(ierr); 967 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 968 Set initial conditions 969 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 970 ierr = SetInitialGuess(X,ctx);CHKERRQ(ierr); 971 972 ierr = VecDuplicate(X,&F_alg);CHKERRQ(ierr); 973 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_alg);CHKERRQ(ierr); 974 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,ctx);CHKERRQ(ierr); 975 ierr = MatZeroEntries(J);CHKERRQ(ierr); 976 ierr = SNESSetJacobian(snes_alg,J,J,AlgJacobian,ctx);CHKERRQ(ierr); 977 ierr = SNESSetOptionsPrefix(snes_alg,"alg_");CHKERRQ(ierr); 978 ierr = SNESSetFromOptions(snes_alg);CHKERRQ(ierr); 979 ctx->alg_flg = PETSC_TRUE; 980 /* Solve the algebraic equations */ 981 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 982 983 /* Just to set up the Jacobian structure */ 984 ierr = VecDuplicate(X,&Xdot);CHKERRQ(ierr); 985 ierr = IJacobian(ts,0.0,X,Xdot,0.0,J,J,ctx);CHKERRQ(ierr); 986 ierr = VecDestroy(&Xdot);CHKERRQ(ierr); 987 988 ctx->stepnum++; 989 990 ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr); 991 ierr = TSSetMaxTime(ts,ctx->tmax);CHKERRQ(ierr); 992 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 993 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 994 995 /* Prefault period */ 996 ctx->alg_flg = PETSC_FALSE; 997 ierr = TSSetTime(ts,0.0);CHKERRQ(ierr); 998 ierr = TSSetMaxTime(ts,ctx->tfaulton);CHKERRQ(ierr); 999 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1000 1001 /* Create the nonlinear solver for solving the algebraic system */ 1002 /* Note that although the algebraic system needs to be solved only for 1003 Idq and V, we reuse the entire system including xgen. The xgen 1004 variables are held constant by setting their residuals to 0 and 1005 putting a 1 on the Jacobian diagonal for xgen rows 1006 */ 1007 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1008 1009 /* Apply disturbance - resistive fault at ctx->faultbus */ 1010 /* This is done by adding shunt conductance to the diagonal location 1011 in the Ybus matrix */ 1012 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; /* Location for G */ 1013 val = 1/ctx->Rfault; 1014 ierr = MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 1015 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; /* Location for G */ 1016 val = 1/ctx->Rfault; 1017 ierr = MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 1018 1019 ierr = MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1020 ierr = MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1021 1022 ctx->alg_flg = PETSC_TRUE; 1023 /* Solve the algebraic equations */ 1024 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 1025 1026 ctx->stepnum++; 1027 1028 /* Disturbance period */ 1029 ctx->alg_flg = PETSC_FALSE; 1030 ierr = TSSetTime(ts,ctx->tfaulton);CHKERRQ(ierr); 1031 ierr = TSSetMaxTime(ts,ctx->tfaultoff);CHKERRQ(ierr); 1032 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1033 1034 /* Remove the fault */ 1035 row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; 1036 val = -1/ctx->Rfault; 1037 ierr = MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 1038 row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; 1039 val = -1/ctx->Rfault; 1040 ierr = MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 1041 1042 ierr = MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1043 ierr = MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1044 1045 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1046 1047 ctx->alg_flg = PETSC_TRUE; 1048 1049 /* Solve the algebraic equations */ 1050 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 1051 1052 ctx->stepnum++; 1053 1054 /* Post-disturbance period */ 1055 ctx->alg_flg = PETSC_TRUE; 1056 ierr = TSSetTime(ts,ctx->tfaultoff);CHKERRQ(ierr); 1057 ierr = TSSetMaxTime(ts,ctx->tmax);CHKERRQ(ierr); 1058 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1059 1060 ierr = VecGetArray(ctx->vec_q,&x_ptr);CHKERRQ(ierr); 1061 *f = x_ptr[0]; 1062 ierr = VecRestoreArray(ctx->vec_q,&x_ptr);CHKERRQ(ierr); 1063 1064 ierr = MatDestroy(&ctx->Ybus);CHKERRQ(ierr); 1065 ierr = VecDestroy(&ctx->V0);CHKERRQ(ierr); 1066 ierr = SNESDestroy(&snes_alg);CHKERRQ(ierr); 1067 ierr = VecDestroy(&F_alg);CHKERRQ(ierr); 1068 ierr = MatDestroy(&J);CHKERRQ(ierr); 1069 ierr = VecDestroy(&X);CHKERRQ(ierr); 1070 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1071 1072 return 0; 1073 } 1074 1075 /*TEST 1076 1077 build: 1078 requires: double !complex !defined(USE_64BIT_INDICES) 1079 1080 test: 1081 args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2 1082 localrunfiles: petscoptions X.bin Ybus.bin 1083 1084 TEST*/ 1085