+ 1import numpy as np
+ 2from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
+ 3from pySDC.implementations.datatype_classes.mesh import mesh
+ 4from pySDC.core.Errors import ParameterError
+ 5
+ 6
+ 7def WSCC9Bus():
+ 8 """
+ 9 Returns the Ybus for the power system.
+ 10
+ 11 Returns
+ 12 -------
+ 13 ppc_res : dict
+ 14 The data with buses, branches, generators (with power flow result) and the Ybus to define the power system.
+ 15 """
+ 16 ppc_res = {}
+ 17
+ 18 ppc_res['baseMVA'] = 100
+ 19 ppc_res['Ybus'] = get_initial_Ybus()
+ 20 ppc_res['bus'] = np.array(
+ 21 [
+ 22 [
+ 23 1.0,
+ 24 3.0,
+ 25 0.0,
+ 26 0.0,
+ 27 0.0,
+ 28 0.0,
+ 29 1.0,
+ 30 1.0,
+ 31 0.0,
+ 32 345.0,
+ 33 1.0,
+ 34 1.100000000000000089e00,
+ 35 9.000000000000000222e-01,
+ 36 ],
+ 37 [
+ 38 2.0,
+ 39 2.0,
+ 40 0.0,
+ 41 0.0,
+ 42 0.0,
+ 43 0.0,
+ 44 1.0,
+ 45 9.999999999999998890e-01,
+ 46 9.668741126628123794e00,
+ 47 345.0,
+ 48 1.0,
+ 49 1.100000000000000089e00,
+ 50 9.000000000000000222e-01,
+ 51 ],
+ 52 [
+ 53 3.0,
+ 54 2.0,
+ 55 0.0,
+ 56 0.0,
+ 57 0.0,
+ 58 0.0,
+ 59 1.0,
+ 60 9.999999999999998890e-01,
+ 61 4.771073237177319015e00,
+ 62 345.0,
+ 63 1.0,
+ 64 1.100000000000000089e00,
+ 65 9.000000000000000222e-01,
+ 66 ],
+ 67 [
+ 68 4.0,
+ 69 1.0,
+ 70 0.0,
+ 71 0.0,
+ 72 0.0,
+ 73 0.0,
+ 74 1.0,
+ 75 9.870068523919054426e-01,
+ 76 -2.406643919519410257e00,
+ 77 345.0,
+ 78 1.0,
+ 79 1.100000000000000089e00,
+ 80 9.000000000000000222e-01,
+ 81 ],
+ 82 [
+ 83 5.0,
+ 84 1.0,
+ 85 90.0,
+ 86 30.0,
+ 87 0.0,
+ 88 0.0,
+ 89 1.0,
+ 90 9.754721770850530715e-01,
+ 91 -4.017264326707549849e00,
+ 92 345.0,
+ 93 1.0,
+ 94 1.100000000000000089e00,
+ 95 9.000000000000000222e-01,
+ 96 ],
+ 97 [
+ 98 6.0,
+ 99 1.0,
+ 100 0.0,
+ 101 0.0,
+ 102 0.0,
+ 103 0.0,
+ 104 1.0,
+ 105 1.003375436452800251e00,
+ 106 1.925601686828564363e00,
+ 107 345.0,
+ 108 1.0,
+ 109 1.100000000000000089e00,
+ 110 9.000000000000000222e-01,
+ 111 ],
+ 112 [
+ 113 7.0,
+ 114 1.0,
+ 115 100.0,
+ 116 35.0,
+ 117 0.0,
+ 118 0.0,
+ 119 1.0,
+ 120 9.856448817249467975e-01,
+ 121 6.215445553889322738e-01,
+ 122 345.0,
+ 123 1.0,
+ 124 1.100000000000000089e00,
+ 125 9.000000000000000222e-01,
+ 126 ],
+ 127 [
+ 128 8.0,
+ 129 1.0,
+ 130 0.0,
+ 131 0.0,
+ 132 0.0,
+ 133 0.0,
+ 134 1.0,
+ 135 9.961852458090698637e-01,
+ 136 3.799120192692319264e00,
+ 137 345.0,
+ 138 1.0,
+ 139 1.100000000000000089e00,
+ 140 9.000000000000000222e-01,
+ 141 ],
+ 142 [
+ 143 9.0,
+ 144 1.0,
+ 145 125.0,
+ 146 50.0,
+ 147 0.0,
+ 148 0.0,
+ 149 1.0,
+ 150 9.576210404299042578e-01,
+ 151 -4.349933576561006987e00,
+ 152 345.0,
+ 153 1.0,
+ 154 1.100000000000000089e00,
+ 155 9.000000000000000222e-01,
+ 156 ],
+ 157 ]
+ 158 )
+ 159 ppc_res['branch'] = np.array(
+ 160 [
+ 161 [
+ 162 1.0,
+ 163 4.0,
+ 164 0.0,
+ 165 5.759999999999999842e-02,
+ 166 0.0,
+ 167 250.0,
+ 168 250.0,
+ 169 250.0,
+ 170 0.0,
+ 171 0.0,
+ 172 1.0,
+ 173 -360.0,
+ 174 360.0,
+ 175 7.195470158922189796e01,
+ 176 2.406895777275899206e01,
+ 177 -7.195470158922189796e01,
+ 178 -2.075304453873995314e01,
+ 179 ],
+ 180 [
+ 181 4.0,
+ 182 5.0,
+ 183 1.700000000000000122e-02,
+ 184 9.199999999999999845e-02,
+ 185 1.580000000000000016e-01,
+ 186 250.0,
+ 187 250.0,
+ 188 250.0,
+ 189 0.0,
+ 190 0.0,
+ 191 1.0,
+ 192 -360.0,
+ 193 360.0,
+ 194 3.072828027973678999e01,
+ 195 -5.858508226424568033e-01,
+ 196 -3.055468555805444453e01,
+ 197 -1.368795049942141873e01,
+ 198 ],
+ 199 [
+ 200 5.0,
+ 201 6.0,
+ 202 3.899999999999999994e-02,
+ 203 1.700000000000000122e-01,
+ 204 3.579999999999999849e-01,
+ 205 150.0,
+ 206 150.0,
+ 207 150.0,
+ 208 0.0,
+ 209 0.0,
+ 210 1.0,
+ 211 -360.0,
+ 212 360.0,
+ 213 -5.944531444194475966e01,
+ 214 -1.631204950057851022e01,
+ 215 6.089386583276659337e01,
+ 216 -1.242746953108854591e01,
+ 217 ],
+ 218 [
+ 219 3.0,
+ 220 6.0,
+ 221 0.0,
+ 222 5.859999999999999931e-02,
+ 223 0.0,
+ 224 300.0,
+ 225 300.0,
+ 226 300.0,
+ 227 0.0,
+ 228 0.0,
+ 229 1.0,
+ 230 -360.0,
+ 231 360.0,
+ 232 8.499999999999997158e01,
+ 233 -3.649025534209526800e00,
+ 234 -8.499999999999997158e01,
+ 235 7.890678351196221740e00,
+ 236 ],
+ 237 [
+ 238 6.0,
+ 239 7.0,
+ 240 1.190000000000000085e-02,
+ 241 1.008000000000000007e-01,
+ 242 2.089999999999999913e-01,
+ 243 150.0,
+ 244 150.0,
+ 245 150.0,
+ 246 0.0,
+ 247 0.0,
+ 248 1.0,
+ 249 -360.0,
+ 250 360.0,
+ 251 2.410613416723325741e01,
+ 252 4.536791179891427994e00,
+ 253 -2.401064777894146474e01,
+ 254 -2.440076244077697609e01,
+ 255 ],
+ 256 [
+ 257 7.0,
+ 258 8.0,
+ 259 8.500000000000000611e-03,
+ 260 7.199999999999999456e-02,
+ 261 1.489999999999999936e-01,
+ 262 250.0,
+ 263 250.0,
+ 264 250.0,
+ 265 0.0,
+ 266 0.0,
+ 267 1.0,
+ 268 -360.0,
+ 269 360.0,
+ 270 -7.598935222105758669e01,
+ 271 -1.059923755922268285e01,
+ 272 7.649556434279409700e01,
+ 273 2.562394697223899231e-01,
+ 274 ],
+ 275 [
+ 276 8.0,
+ 277 2.0,
+ 278 0.0,
+ 279 6.250000000000000000e-02,
+ 280 0.0,
+ 281 250.0,
+ 282 250.0,
+ 283 250.0,
+ 284 0.0,
+ 285 0.0,
+ 286 1.0,
+ 287 -360.0,
+ 288 360.0,
+ 289 -1.629999999999997158e02,
+ 290 2.276189879408803574e00,
+ 291 1.629999999999997442e02,
+ 292 1.446011953112515869e01,
+ 293 ],
+ 294 [
+ 295 8.0,
+ 296 9.0,
+ 297 3.200000000000000067e-02,
+ 298 1.610000000000000042e-01,
+ 299 3.059999999999999942e-01,
+ 300 250.0,
+ 301 250.0,
+ 302 250.0,
+ 303 0.0,
+ 304 0.0,
+ 305 1.0,
+ 306 -360.0,
+ 307 360.0,
+ 308 8.650443565720313188e01,
+ 309 -2.532429349130207452e00,
+ 310 -8.403988686535042518e01,
+ 311 -1.428198298779915731e01,
+ 312 ],
+ 313 [
+ 314 9.0,
+ 315 4.0,
+ 316 1.000000000000000021e-02,
+ 317 8.500000000000000611e-02,
+ 318 1.759999999999999898e-01,
+ 319 250.0,
+ 320 250.0,
+ 321 250.0,
+ 322 0.0,
+ 323 0.0,
+ 324 1.0,
+ 325 -360.0,
+ 326 360.0,
+ 327 -4.096011313464404680e01,
+ 328 -3.571801701219811775e01,
+ 329 4.122642130948177197e01,
+ 330 2.133889536138378062e01,
+ 331 ],
+ 332 ]
+ 333 )
+ 334 ppc_res['gen'] = np.array(
+ 335 [
+ 336 [
+ 337 1.0,
+ 338 71.0,
+ 339 24.0,
+ 340 300.0,
+ 341 -300.0,
+ 342 1.0,
+ 343 100.0,
+ 344 1.0,
+ 345 250.0,
+ 346 10.0,
+ 347 0.0,
+ 348 0.0,
+ 349 0.0,
+ 350 0.0,
+ 351 0.0,
+ 352 0.0,
+ 353 0.0,
+ 354 0.0,
+ 355 0.0,
+ 356 0.0,
+ 357 0.0,
+ 358 ],
+ 359 [
+ 360 2.0,
+ 361 163.0,
+ 362 14.0,
+ 363 300.0,
+ 364 -300.0,
+ 365 1.0,
+ 366 100.0,
+ 367 1.0,
+ 368 300.0,
+ 369 10.0,
+ 370 0.0,
+ 371 0.0,
+ 372 0.0,
+ 373 0.0,
+ 374 0.0,
+ 375 0.0,
+ 376 0.0,
+ 377 0.0,
+ 378 0.0,
+ 379 0.0,
+ 380 0.0,
+ 381 ],
+ 382 [
+ 383 3.0,
+ 384 85.0,
+ 385 -3.0,
+ 386 300.0,
+ 387 -300.0,
+ 388 1.0,
+ 389 100.0,
+ 390 1.0,
+ 391 270.0,
+ 392 10.0,
+ 393 0.0,
+ 394 0.0,
+ 395 0.0,
+ 396 0.0,
+ 397 0.0,
+ 398 0.0,
+ 399 0.0,
+ 400 0.0,
+ 401 0.0,
+ 402 0.0,
+ 403 0.0,
+ 404 ],
+ 405 ]
+ 406 )
+ 407
+ 408 return ppc_res
+ 409
+ 410
+ 411def get_initial_Ybus():
+ 412 ybus = np.array(
+ 413 [
+ 414 [0 - 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j],
+ 415 [0 + 0j, 0 - 16j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 16j, 0 + 0j],
+ 416 [0 + 0j, 0 + 0j, 0 - 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 0j],
+ 417 [
+ 418 0 + 17.36111111111111j,
+ 419 0 + 0j,
+ 420 0 + 0j,
+ 421 3.307378962025306 - 39.30888872611897j,
+ 422 -1.942191248714727 + 10.51068205186793j,
+ 423 0 + 0j,
+ 424 0 + 0j,
+ 425 0 + 0j,
+ 426 -1.36518771331058 + 11.60409556313993j,
+ 427 ],
+ 428 [
+ 429 0 + 0j,
+ 430 0 + 0j,
+ 431 0 + 0j,
+ 432 -1.942191248714727 + 10.51068205186793j,
+ 433 3.224200387138842 - 15.84092701422946j,
+ 434 -1.282009138424115 + 5.588244962361526j,
+ 435 0 + 0j,
+ 436 0 + 0j,
+ 437 0 + 0j,
+ 438 ],
+ 439 [
+ 440 0 + 0j,
+ 441 0 + 0j,
+ 442 0 + 17.06484641638225j,
+ 443 0 + 0j,
+ 444 -1.282009138424115 + 5.588244962361526j,
+ 445 2.437096619314212 - 32.15386180510696j,
+ 446 -1.155087480890097 + 9.784270426363173j,
+ 447 0 + 0j,
+ 448 0 + 0j,
+ 449 ],
+ 450 [
+ 451 0 + 0j,
+ 452 0 + 0j,
+ 453 0 + 0j,
+ 454 0 + 0j,
+ 455 0 + 0j,
+ 456 -1.155087480890097 + 9.784270426363173j,
+ 457 2.772209954136233 - 23.30324902327162j,
+ 458 -1.617122473246136 + 13.69797859690844j,
+ 459 0 + 0j,
+ 460 ],
+ 461 [
+ 462 0 + 0j,
+ 463 0 + 16j,
+ 464 0 + 0j,
+ 465 0 + 0j,
+ 466 0 + 0j,
+ 467 0 + 0j,
+ 468 -1.617122473246136 + 13.69797859690844j,
+ 469 2.804726852537284 - 35.44561313021703j,
+ 470 -1.187604379291148 + 5.975134533308591j,
+ 471 ],
+ 472 [
+ 473 0 + 0j,
+ 474 0 + 0j,
+ 475 0 + 0j,
+ 476 -1.36518771331058 + 11.60409556313993j,
+ 477 0 + 0j,
+ 478 0 + 0j,
+ 479 0 + 0j,
+ 480 -1.187604379291148 + 5.975134533308591j,
+ 481 2.552792092601728 - 17.33823009644852j,
+ 482 ],
+ 483 ],
+ 484 dtype=complex,
+ 485 )
+ 486
+ 487 return ybus
+ 488
+ 489
+ 490def get_event_Ybus():
+ 491 ybus = np.array(
+ 492 [
+ 493 [0 - 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j],
+ 494 [0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j],
+ 495 [0 + 0j, 0 + 0j, 0 - 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 17.06484641638225j],
+ 496 [
+ 497 0 + 17.36111111111111j,
+ 498 0 + 0j,
+ 499 0 + 0j,
+ 500 3.307378962025306 - 39.30888872611897j,
+ 501 -1.36518771331058 + 11.60409556313993j,
+ 502 -1.942191248714727 + 10.51068205186793j,
+ 503 0 + 0j,
+ 504 0 + 0j,
+ 505 0 + 0j,
+ 506 ],
+ 507 [
+ 508 0 + 0j,
+ 509 0 + 0j,
+ 510 0 + 0j,
+ 511 -1.36518771331058 + 11.60409556313993j,
+ 512 2.552792092601728 - 17.33823009644852j,
+ 513 0 + 0j,
+ 514 -1.187604379291148 + 5.975134533308591j,
+ 515 0 + 0j,
+ 516 0 + 0j,
+ 517 ],
+ 518 [
+ 519 0 + 0j,
+ 520 0 + 0j,
+ 521 0 + 0j,
+ 522 -1.942191248714727 + 10.51068205186793j,
+ 523 0 + 0j,
+ 524 3.224200387138842 - 15.84092701422946j,
+ 525 0 + 0j,
+ 526 0 + 0j,
+ 527 -1.282009138424115 + 5.588244962361526j,
+ 528 ],
+ 529 [
+ 530 0 + 0j,
+ 531 0 + 0j,
+ 532 0 + 0j,
+ 533 0 + 0j,
+ 534 -1.187604379291148 + 5.975134533308591j,
+ 535 0 + 0j,
+ 536 2.804726852537284 - 19.44561313021703j,
+ 537 -1.617122473246136 + 13.69797859690844j,
+ 538 0 + 0j,
+ 539 ],
+ 540 [
+ 541 0 + 0j,
+ 542 0 + 0j,
+ 543 0 + 0j,
+ 544 0 + 0j,
+ 545 0 + 0j,
+ 546 0 + 0j,
+ 547 -1.617122473246136 + 13.69797859690844j,
+ 548 2.772209954136233 - 23.30324902327162j,
+ 549 -1.155087480890097 + 9.784270426363173j,
+ 550 ],
+ 551 [
+ 552 0 + 0j,
+ 553 0 + 0j,
+ 554 0 + 17.06484641638225j,
+ 555 0 + 0j,
+ 556 0 + 0j,
+ 557 -1.282009138424115 + 5.588244962361526j,
+ 558 0 + 0j,
+ 559 -1.155087480890097 + 9.784270426363173j,
+ 560 2.437096619314212 - 32.15386180510696j,
+ 561 ],
+ 562 ],
+ 563 dtype=complex,
+ 564 )
+ 565
+ 566 return ybus
+ 567
+ 568
+ 569# def get_YBus(ppc):
+ 570
+ 571# ppci = ext2int(ppc)
+ 572# Ybus, yf, yt = makeYbus(ppci['baseMVA'],ppci['bus'],ppci['branch'])
+ 573
+ 574# return ppc['Ybus']()
+ 575
+ 576
+ 577class WSCC9BusSystem(ptype_dae):
+ 578 r"""
+ 579 Example implementing the WSCC 9 Bus system [1]_. For this complex model, the equations can be found in [2]_, and
+ 580 sub-transient and turbine parameters are taken from [3]_. The network data of the system are taken from MATPOWER [4]_.
+ 581
+ 582 Parameters
+ 583 ----------
+ 584 nvars : int
+ 585 Number of unknowns of the system of DAEs (not used here, since it is set up inside this class).
+ 586 newton_tol : float
+ 587 Tolerance for Newton solver.
+ 588
+ 589 Attributes
+ 590 ----------
+ 591 mpc : dict
+ 592 Contains the data for the buses, branches, generators, and the Ybus.
+ 593 m : int
+ 594 Number of machines used in the network.
+ 595 n : int
+ 596 Number of buses used in the network.
+ 597 baseMVA : float
+ 598 Base value of the apparent power.
+ 599 ws : float
+ 600 Generator synchronous speed in rad per second.
+ 601 ws_vector : np.1darray
+ 602 Vector containing ``ws``.
+ 603 MD : np.2darray
+ 604 Machine data.
+ 605 ED : np.2darray
+ 606 Excitation data.
+ 607 TD : np.2darray
+ 608 Turbine data.
+ 609 bus : np.2darray
+ 610 Data for the buses.
+ 611 branch : np.2darray
+ 612 Data for the branches in the power system.
+ 613 gen : np.2darray
+ 614 Data for generators in the system.
+ 615 Ybus : np.2darray
+ 616 Ybus.
+ 617 YBus_line6_8_outage : np.2darray
+ 618 Contains the data for the line outage in the power system, where line at bus6 is outaged.
+ 619 psv_max : float
+ 620 Maximum valve position.
+ 621 IC1 : list
+ 622 Contains the :math:`8`-th row of the ``bus`` matrix.
+ 623 IC2 : list
+ 624 Contains the :math:`9`-th row of the ``bus`` matrix.
+ 625 IC3 : list
+ 626 Generator values divided by ``baseMVA``.
+ 627 IC4 : list
+ 628 Generator values divided by ``baseMVA``.
+ 629 IC5 : list
+ 630 Loads divided by ``baseMVA``.
+ 631 IC6 : list
+ 632 Loads divided by ``baseMVA``.
+ 633 genP : list
+ 634 Contains data for one generator from the ``mpc``.
+ 635 IC : list
+ 636 Contains all the values of `IC1, IC2, IC3, IC4, IC5, IC6`.
+ 637 PL : list
+ 638 Contains the :math:`5`-th column of ``IC``.
+ 639 QL : list
+ 640 Contains the :math:`6`-th column of ``IC``.
+ 641 PG : np.1darray
+ 642 Contains the :math:`3`-rd column of ``IC``.
+ 643 QG : np.1darray
+ 644 Contains the :math:`4`-th column of ``IC``.
+ 645 TH0 : np.1darray
+ 646 Initial condition for angle of bus voltage in rad.
+ 647 V0 : np.1darray
+ 648 Contains the :math:`1`-st column of ``IC``, initial condition for magnitude of bus voltage in per unit.
+ 649 VG0 : np.1darray
+ 650 Initial condition for complex voltage phasor.
+ 651 THG0 : np.1darray
+ 652 Initial condition for angle of the bus voltage in rad.
+ 653 H : np.1darray
+ 654 Shaft inertia constant in second.
+ 655 Xd : np.1darray
+ 656 d-axis reactance in per unit.
+ 657 Xdp : np.1darray
+ 658 Transient d-axis reactance in per unit.
+ 659 Xdpp : np.1darray
+ 660 Sub-transient d-axis reactance in per unit.
+ 661 Xq : np.1darray
+ 662 q-axis reactance in per unit.
+ 663 Xqp : np.1darray
+ 664 Transient q-axis reactance in per unit.
+ 665 Xqpp : np.1darray
+ 666 Sub-transient q-axis reactance in per unit.
+ 667 Td0p : np.1darray
+ 668 d-axis time constant associated with :math:`E_q'` in second.
+ 669 Td0pp : np.1darray
+ 670 d-axis time constant associated with :math:`\psi_{1d}` in second.
+ 671 Tq0p : np.1darray
+ 672 q-axis time constant associated with :math:`E_d'` in second.
+ 673 Tq0pp : np.1darray
+ 674 q-axis time constant associated with :math:`\psi_{2q}` in second.
+ 675 Rs : np.1darray
+ 676 Stator resistance in per unit.
+ 677 Xls : np.1darray
+ 678 Parameter :math:`X_{\ell s}`.
+ 679 Dm : np.1darray
+ 680 Rotor angle in rad.
+ 681 KA : np.1darray
+ 682 Amplifier gain.
+ 683 TA : np.1darray
+ 684 Amplifier time constant in second.
+ 685 KE : np.1darray
+ 686 Separate or self-excited constant.
+ 687 TE : np.1darray
+ 688 Parameter :math:`T_E`.
+ 689 KF : np.1darray
+ 690 Parameter _math:`K_F`.
+ 691 TF : np.1darray
+ 692 Parameter :math:`T_F`.
+ 693 Ax : np.1darray
+ 694 Constant :math:`A_x` of the saturation function :math:`S_{E_i}`.
+ 695 Bx : np.1darray
+ 696 Constant :math:`B_x` of the saturation function :math:`S_{E_i}`.
+ 697 TCH : np.1darray
+ 698 Incremental steam chest time constant in second.
+ 699 TSV : np.1darray
+ 700 Steam valve time constant in second.
+ 701 RD : np.1darray
+ 702 Speed regulation quantity in Hz/per unit.
+ 703 MH : float
+ 704 Factor :math:`\frac{2 H_i}{w_s}`.
+ 705 QG : np.1darray
+ 706 Used to compute :math:`I_{phasor}`.
+ 707 Vphasor : np.1darray
+ 708 Complex voltage phasor.
+ 709 Iphasor : np.1darray
+ 710 Complex current phasor.
+ 711 E0 : np.1darray
+ 712 Initial internal voltage of the synchronous generator.
+ 713 Em : np.1darray
+ 714 Absolute values of ``E0``.
+ 715 D0 : np.1darray
+ 716 Initial condition for rotor angle in rad.
+ 717 Id0 : np.1darray
+ 718 Initial condition for d-axis current in per unit.
+ 719 Iq0 : np.1darray
+ 720 Initial condition for q-axis current in per unit.
+ 721 Edp0 : np.1darray
+ 722 Initial condition for d-axis transient internal voltages in per unit.
+ 723 Si2q0 : np.1darray
+ 724 Initial condition for damper winding 2q flux linkages in per unit.
+ 725 Eqp0 : np.1darray
+ 726 Initial condition for q-axis transient internal voltages in per unit.
+ 727 Si1d0 : np.1darray
+ 728 Initial condition for damper winding 1d flux linkages in per unit.
+ 729 Efd0 : np.1darray
+ 730 Initial condition for field winding fd flux linkages in per unit.
+ 731 TM0 : np.1darray
+ 732 Initial condition for mechanical input torque in per unit.
+ 733 VR0 : np.1darray
+ 734 Initial condition for exciter input in per unit.
+ 735 RF0 : np.1darray
+ 736 Initial condition for exciter feedback in per unit.
+ 737 Vref : np.1darray
+ 738 Reference voltage input in per unit.
+ 739 PSV0 : np.1darray
+ 740 Initial condition for steam valve position in per unit.
+ 741 PC : np.1darray
+ 742 Initial condition for control power input in per unit.
+ 743 alpha : int
+ 744 Active load parameter.
+ 745 beta : int
+ 746 Reactive load parameter.
+ 747 bb1, aa1 : list of ndarrays
+ 748 Used to access on specific values of ``TH``.
+ 749 bb2, aa2 : list of ndarrays
+ 750 Used to access on specific values of ``TH``.
+ 751 t_switch : float
+ 752 Time the event found by detection.
+ 753 nswitches : int
+ 754 Number of events found by detection.
+ 755
+ 756 References
+ 757 ----------
+ 758 .. [1] WSCC 9-Bus System - Illinois Center for a Smarter Electric Grid. https://icseg.iti.illinois.edu/wscc-9-bus-system/
+ 759 .. [2] P. W. Sauer, M. A. Pai. Power System Dynamics and Analysis. John Wiley & Sons (2008).
+ 760 .. [3] I. Abdulrahman. MATLAB-Based Programs for Power System Dynamics Analysis. IEEE Open Access Journal of Power and Energy.
+ 761 Vol. 7, No. 1, pp. 59–69 (2020).
+ 762 .. [4] R. D. Zimmerman, C. E. Murillo-Sánchez, R. J. Thomas. MATPOWER: Steady-State Operations, Planning, and Analysis Tools
+ 763 for Power Systems Research and Education. IEEE Transactions on Power Systems. Vol. 26, No. 1, pp. 12–19 (2011).
+ 764 """
+ 765
+ 766 dtype_u = mesh
+ 767 dtype_f = mesh
+ 768
+ 769 def __init__(self, nvars=None, newton_tol=1e-10, m=3, n=9):
+ 770 """Initialization routine"""
+ 771
+ 772 nvars = 11 * m + 2 * m + 2 * n
+ 773 # invoke super init, passing number of dofs
+ 774 super().__init__(nvars, newton_tol)
+ 775 self._makeAttributeAndRegister('nvars', 'newton_tol', localVars=locals(), readOnly=True)
+ 776 self._makeAttributeAndRegister('m', 'n', localVars=locals())
+ 777 self.mpc = WSCC9Bus()
+ 778
+ 779 self.baseMVA = self.mpc['baseMVA']
+ 780 self.ws = 2 * np.pi * 60
+ 781 self.ws_vector = self.ws * np.ones(self.m)
+ 782
+ 783 # Machine data (MD) as a 2D NumPy array
+ 784 self.MD = np.array(
+ 785 [
+ 786 [23.640, 6.4000, 3.0100], # 1 - H
+ 787 [0.1460, 0.8958, 1.3125], # 2 - Xd
+ 788 [0.0608, 0.1198, 0.1813], # 3 - Xdp
+ 789 [0.0489, 0.0881, 0.1133], # 4 - Xdpp
+ 790 [0.0969, 0.8645, 1.2578], # 5 - Xq
+ 791 [0.0969, 0.1969, 0.2500], # 6 - Xqp
+ 792 [0.0396, 0.0887, 0.0833], # 7 - Xqpp
+ 793 [8.960000000000001, 6.0000, 5.8900], # 8 - Tdop
+ 794 [0.1150, 0.0337, 0.0420], # 9 - Td0pp
+ 795 [0.3100, 0.5350, 0.6000], # 10 - Tqop
+ 796 [0.0330, 0.0780, 0.1875], # 11 - Tq0pp
+ 797 [0.0041, 0.0026, 0.0035], # 12 - RS
+ 798 [0.1200, 0.1020, 0.0750], # 13 - Xls
+ 799 [
+ 800 0.1 * (2 * 23.64) / self.ws,
+ 801 0.2 * (2 * 6.4) / self.ws,
+ 802 0.3 * (2 * 3.01) / self.ws,
+ 803 ], # 14 - Dm (ws should be defined)
+ 804 ]
+ 805 )
+ 806
+ 807 # Excitation data (ED) as a 2D NumPy array
+ 808 self.ED = np.array(
+ 809 [
+ 810 20.000 * np.ones(self.m), # 1 - KA
+ 811 0.2000 * np.ones(self.m), # 2 - TA
+ 812 1.0000 * np.ones(self.m), # 3 - KE
+ 813 0.3140 * np.ones(self.m), # 4 - TE
+ 814 0.0630 * np.ones(self.m), # 5 - KF
+ 815 0.3500 * np.ones(self.m), # 6 - TF
+ 816 0.0039 * np.ones(self.m), # 7 - Ax
+ 817 1.5550 * np.ones(self.m), # 8 - Bx
+ 818 ]
+ 819 )
+ 820
+ 821 # Turbine data (TD) as a 2D NumPy array
+ 822 self.TD = np.array(
+ 823 [
+ 824 0.10 * np.ones(self.m), # 1 - TCH
+ 825 0.05 * np.ones(self.m), # 2 - TSV
+ 826 0.05 * np.ones(self.m), # 3 - RD
+ 827 ]
+ 828 )
+ 829
+ 830 self.bus = self.mpc['bus']
+ 831 self.branch = self.mpc['branch']
+ 832 self.gen = self.mpc['gen']
+ 833 self.YBus = get_initial_Ybus()
+ 834
+ 835 temp_mpc = self.mpc
+ 836 temp_mpc['branch'] = np.delete(
+ 837 temp_mpc['branch'], 6, 0
+ 838 ) # note that this is correct but not necessary, because we have the event Ybus already
+ 839 self.YBus_line6_8_outage = get_event_Ybus()
+ 840
+ 841 # excitation limiter vmax
+ 842 # self.vmax = 2.1
+ 843 self.psv_max = 1.0
+ 844
+ 845 self.IC1 = [row[7] for row in self.bus] # Column 8 in MATLAB is indexed as 7 in Python (0-based index)
+ 846 self.IC2 = [row[8] for row in self.bus] # Column 9 in MATLAB is indexed as 8 in Python
+ 847
+ 848 n_prev, m_prev = self.n, self.m
+ 849 self.n = len(self.bus) # Number of rows in 'bus' list; self.n already defined above?!
+ 850 self.m = len(self.gen) # Number of rows in 'gen' list; self.m already defined above?!
+ 851 if n_prev != 9 or m_prev != 3:
+ 852 raise ParameterError("Number of rows in bus or gen not equal to initialised n or m!")
+ 853
+ 854 gen0 = [0] * self.n
+ 855 for i in range(self.m):
+ 856 gen0[i] = self.gen[i][1]
+ 857 self.genP = gen0
+ 858 self.IC3 = [val / self.baseMVA for val in self.genP]
+ 859
+ 860 gen0 = [0] * self.n
+ 861 for i in range(self.m):
+ 862 gen0[i] = self.gen[i][2]
+ 863 genQ = gen0
+ 864 for i in range(self.n):
+ 865 genQ[i] += self.bus[i][5] # Column 6 in MATLAB is indexed as 5 in Python
+ 866 self.IC4 = [val / self.baseMVA for val in genQ]
+ 867
+ 868 self.IC5 = [row[2] for row in self.bus] # Column 3 in MATLAB is indexed as 2 in Python
+ 869 self.IC5 = [val / self.baseMVA for val in self.IC5]
+ 870
+ 871 self.IC6 = [row[3] for row in self.bus] # Column 4 in MATLAB is indexed as 3 in Python
+ 872 self.IC6 = [val / self.baseMVA for val in self.IC6]
+ 873
+ 874 self.IC = list(zip(self.IC1, self.IC2, self.IC3, self.IC4, self.IC5, self.IC6))
+ 875
+ 876 self.PL = [row[4] for row in self.IC] # Column 5 in MATLAB is indexed as 4 in Python
+ 877 self.QL = [row[5] for row in self.IC] # Column 6 in MATLAB is indexed as 5 in Python
+ 878
+ 879 self.PG = np.array([row[2] for row in self.IC]) # Column 3 in MATLAB is indexed as 2 in Python
+ 880 self.QG = np.array([row[3] for row in self.IC]) # Column 4 in MATLAB is indexed as 3 in Python
+ 881
+ 882 self.TH0 = np.array([row[1] * np.pi / 180 for row in self.IC])
+ 883 self.V0 = np.array([row[0] for row in self.IC])
+ 884 self.VG0 = self.V0[: self.m]
+ 885 self.THG0 = self.TH0[: self.m]
+ 886
+ 887 # Extracting values from the MD array
+ 888 self.H = self.MD[0, :]
+ 889 self.Xd = self.MD[1, :]
+ 890 self.Xdp = self.MD[2, :]
+ 891 self.Xdpp = self.MD[3, :]
+ 892 self.Xq = self.MD[4, :]
+ 893 self.Xqp = self.MD[5, :]
+ 894 self.Xqpp = self.MD[6, :]
+ 895 self.Td0p = self.MD[7, :]
+ 896 self.Td0pp = self.MD[8, :]
+ 897 self.Tq0p = self.MD[9, :]
+ 898 self.Tq0pp = self.MD[10, :]
+ 899 self.Rs = self.MD[11, :]
+ 900 self.Xls = self.MD[12, :]
+ 901 self.Dm = self.MD[13, :]
+ 902
+ 903 # Extracting values from the ED array
+ 904 self.KA = self.ED[0, :]
+ 905 self.TA = self.ED[1, :]
+ 906 self.KE = self.ED[2, :]
+ 907 self.TE = self.ED[3, :]
+ 908 self.KF = self.ED[4, :]
+ 909 self.TF = self.ED[5, :]
+ 910 self.Ax = self.ED[6, :]
+ 911 self.Bx = self.ED[7, :]
+ 912
+ 913 # Extracting values from the TD array
+ 914 self.TCH = self.TD[0, :]
+ 915 self.TSV = self.TD[1, :]
+ 916 self.RD = self.TD[2, :]
+ 917
+ 918 # Calculate MH
+ 919 self.MH = 2 * self.H / self.ws
+ 920
+ 921 # Represent QG as complex numbers
+ 922 self.QG = self.QG.astype(complex)
+ 923
+ 924 # Convert VG0 and THG0 to complex phasors
+ 925 self.Vphasor = self.VG0 * np.exp(1j * self.THG0)
+ 926
+ 927 # Calculate Iphasor
+ 928 self.Iphasor = np.conj(np.divide(self.PG[:m] + 1j * self.QG[:m], self.Vphasor))
+ 929
+ 930 # Calculate E0
+ 931 self.E0 = self.Vphasor + (self.Rs + 1j * self.Xq) * self.Iphasor
+ 932
+ 933 # Calculate Em, D0, Id0, and Iq0
+ 934 self.Em = np.abs(self.E0)
+ 935 self.D0 = np.angle(self.E0)
+ 936 self.Id0 = np.real(self.Iphasor * np.exp(-1j * (self.D0 - np.pi / 2)))
+ 937 self.Iq0 = np.imag(self.Iphasor * np.exp(-1j * (self.D0 - np.pi / 2)))
+ 938
+ 939 # Calculate Edp0, Si2q0, Eqp0, and Si1d0
+ 940 self.Edp0 = (self.Xq - self.Xqp) * self.Iq0
+ 941 self.Si2q0 = (self.Xls - self.Xq) * self.Iq0
+ 942 self.Eqp0 = self.Rs * self.Iq0 + self.Xdp * self.Id0 + self.V0[: self.m] * np.cos(self.D0 - self.TH0[: self.m])
+ 943 self.Si1d0 = self.Eqp0 - (self.Xdp - self.Xls) * self.Id0
+ 944
+ 945 # Calculate Efd0 and TM0
+ 946 self.Efd0 = self.Eqp0 + (self.Xd - self.Xdp) * self.Id0
+ 947 self.TM0 = (
+ 948 ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * self.Eqp0 * self.Iq0
+ 949 + ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * self.Si1d0 * self.Iq0
+ 950 + ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * self.Edp0 * self.Id0
+ 951 - ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * self.Si2q0 * self.Id0
+ 952 + (self.Xqpp - self.Xdpp) * self.Id0 * self.Iq0
+ 953 )
+ 954
+ 955 # Calculate VR0 and RF0
+ 956 self.VR0 = (self.KE + self.Ax * np.exp(self.Bx * self.Efd0)) * self.Efd0
+ 957 self.RF0 = (self.KF / self.TF) * self.Efd0
+ 958
+ 959 # Calculate Vref and PSV0
+ 960 self.Vref = self.V0[: self.m] + self.VR0 / self.KA
+ 961 self.PSV0 = self.TM0
+ 962 self.PC = self.PSV0
+ 963
+ 964 self.alpha = 2
+ 965 self.beta = 2
+ 966
+ 967 self.bb1, self.aa1 = np.meshgrid(np.arange(0, self.m), np.arange(0, self.n))
+ 968 self.bb1, self.aa1 = self.bb1.astype(int), self.aa1.astype(int)
+ 969
+ 970 # Create matrix grid to eliminate for-loops (load buses)
+ 971 self.bb2, self.aa2 = np.meshgrid(np.arange(self.m, self.n), np.arange(0, self.n))
+ 972 self.bb2, self.aa2 = self.bb2.astype(int), self.aa2.astype(int)
+ 973
+ 974 self.t_switch = None
+ 975 self.nswitches = 0
+ 976
+ 977 def eval_f(self, u, du, t):
+ 978 r"""
+ 979 Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`.
+ 980
+ 981 Parameters
+ 982 ----------
+ 983 u : dtype_u
+ 984 Current values of the numerical solution at time t.
+ 985 du : dtype_u
+ 986 Current values of the derivative of the numerical solution at time t.
+ 987 t : float
+ 988 Current time of the numerical solution.
+ 989
+ 990 Returns
+ 991 -------
+ 992 f : dtype_f
+ 993 The right-hand side of f (contains two components).
+ 994 """
+ 995
+ 996 dEqp, dSi1d, dEdp = du[0 : self.m], du[self.m : 2 * self.m], du[2 * self.m : 3 * self.m]
+ 997 dSi2q, dDelta = du[3 * self.m : 4 * self.m], du[4 * self.m : 5 * self.m]
+ 998 dw, dEfd, dRF = du[5 * self.m : 6 * self.m], du[6 * self.m : 7 * self.m], du[7 * self.m : 8 * self.m]
+ 999 dVR, dTM, dPSV = du[8 * self.m : 9 * self.m], du[9 * self.m : 10 * self.m], du[10 * self.m : 11 * self.m]
+ 1000
+ 1001 Eqp, Si1d, Edp = u[0 : self.m], u[self.m : 2 * self.m], u[2 * self.m : 3 * self.m]
+ 1002 Si2q, Delta = u[3 * self.m : 4 * self.m], u[4 * self.m : 5 * self.m]
+ 1003 w, Efd, RF = u[5 * self.m : 6 * self.m], u[6 * self.m : 7 * self.m], u[7 * self.m : 8 * self.m]
+ 1004 VR, TM, PSV = u[8 * self.m : 9 * self.m], u[9 * self.m : 10 * self.m], u[10 * self.m : 11 * self.m]
+ 1005
+ 1006 Id, Iq = u[11 * self.m : 11 * self.m + self.m], u[11 * self.m + self.m : 11 * self.m + 2 * self.m]
+ 1007 V = u[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n]
+ 1008 TH = u[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n]
+ 1009
+ 1010 # line outage disturbance:
+ 1011 if t >= 0.05:
+ 1012 self.YBus = self.YBus_line6_8_outage
+ 1013
+ 1014 self.Yang = np.angle(self.YBus)
+ 1015 self.Yabs = np.abs(self.YBus)
+ 1016
+ 1017 COI = np.sum(w * self.MH) / np.sum(self.MH)
+ 1018
+ 1019 # Voltage-dependent active loads PL2, and voltage-dependent reactive loads QL2
+ 1020 PL2 = np.array(self.PL)
+ 1021 QL2 = np.array(self.QL)
+ 1022
+ 1023 V = V.T
+ 1024
+ 1025 # Vectorized calculations
+ 1026 Vectorized_angle1 = (
+ 1027 np.array([TH.take(indices) for indices in self.bb1.T])
+ 1028 - np.array([TH.take(indices) for indices in self.aa1.T])
+ 1029 - self.Yang[: self.m, : self.n]
+ 1030 )
+ 1031 Vectorized_mag1 = (V[: self.m] * V[: self.n].reshape(-1, 1)).T * self.Yabs[: self.m, : self.n]
+ 1032
+ 1033 sum1 = np.sum(Vectorized_mag1 * np.cos(Vectorized_angle1), axis=1)
+ 1034 sum2 = np.sum(Vectorized_mag1 * np.sin(Vectorized_angle1), axis=1)
+ 1035
+ 1036 VG = V[: self.m]
+ 1037 THG = TH[: self.m]
+ 1038 Angle_diff = Delta - THG
+ 1039
+ 1040 Vectorized_angle2 = (
+ 1041 np.array([TH.take(indices) for indices in self.bb2.T])
+ 1042 - np.array([TH.take(indices) for indices in self.aa2.T])
+ 1043 - self.Yang[self.m : self.n, : self.n]
+ 1044 )
+ 1045 Vectorized_mag2 = (V[self.m : self.n] * V[: self.n].reshape(-1, 1)).T * self.Yabs[self.m : self.n, : self.n]
+ 1046
+ 1047 sum3 = np.sum(Vectorized_mag2 * np.cos(Vectorized_angle2), axis=1)
+ 1048 sum4 = np.sum(Vectorized_mag2 * np.sin(Vectorized_angle2), axis=1)
+ 1049
+ 1050 # Initialise f
+ 1051 f = self.dtype_f(self.init)
+ 1052
+ 1053 t_switch = np.inf if self.t_switch is None else self.t_switch
+ 1054
+ 1055 # Equations as list
+ 1056 eqs = []
+ 1057 eqs.append(
+ 1058 (1.0 / self.Td0p)
+ 1059 * (
+ 1060 -Eqp
+ 1061 - (self.Xd - self.Xdp)
+ 1062 * (
+ 1063 Id
+ 1064 - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls) ** 2) * (Si1d + (self.Xdp - self.Xls) * Id - Eqp)
+ 1065 )
+ 1066 + Efd
+ 1067 )
+ 1068 - dEqp
+ 1069 ) # (1)
+ 1070 eqs.append((1.0 / self.Td0pp) * (-Si1d + Eqp - (self.Xdp - self.Xls) * Id) - dSi1d) # (2)
+ 1071 eqs.append(
+ 1072 (1.0 / self.Tq0p)
+ 1073 * (
+ 1074 -Edp
+ 1075 + (self.Xq - self.Xqp)
+ 1076 * (
+ 1077 Iq
+ 1078 - ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls) ** 2) * (Si2q + (self.Xqp - self.Xls) * Iq + Edp)
+ 1079 )
+ 1080 )
+ 1081 - dEdp
+ 1082 ) # (3)
+ 1083 eqs.append((1.0 / self.Tq0pp) * (-Si2q - Edp - (self.Xqp - self.Xls) * Iq) - dSi2q) # (4)
+ 1084 eqs.append(w - COI - dDelta) # (5)
+ 1085 eqs.append(
+ 1086 (self.ws / (2.0 * self.H))
+ 1087 * (
+ 1088 TM
+ 1089 - ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * Eqp * Iq
+ 1090 - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * Si1d * Iq
+ 1091 - ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * Edp * Id
+ 1092 + ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * Si2q * Id
+ 1093 - (self.Xqpp - self.Xdpp) * Id * Iq
+ 1094 - self.Dm * (w - self.ws)
+ 1095 )
+ 1096 - dw
+ 1097 ) # (6)
+ 1098 eqs.append((1.0 / self.TE) * ((-(self.KE + self.Ax * np.exp(self.Bx * Efd))) * Efd + VR) - dEfd) # (7)
+ 1099 eqs.append((1.0 / self.TF) * (-RF + (self.KF / self.TF) * Efd) - dRF) # (8)
+ 1100 eqs.append(
+ 1101 (1.0 / self.TA)
+ 1102 * (-VR + self.KA * RF - ((self.KA * self.KF) / self.TF) * Efd + self.KA * (self.Vref - V[: self.m]))
+ 1103 - dVR
+ 1104 ) # (9)
+ 1105
+ 1106 # Limitation of valve position Psv with limiter start
+ 1107 if PSV[0] >= self.psv_max or t >= t_switch:
+ 1108 _temp_dPSV_g1 = (1.0 / self.TSV[1]) * (
+ 1109 -PSV[1] + self.PSV0[1] - (1.0 / self.RD[1]) * (w[1] / self.ws - 1)
+ 1110 ) - dPSV[1]
+ 1111 _temp_dPSV_g2 = (1.0 / self.TSV[2]) * (
+ 1112 -PSV[2] + self.PSV0[2] - (1.0 / self.RD[2]) * (w[2] / self.ws - 1)
+ 1113 ) - dPSV[2]
+ 1114 eqs.append(np.array([dPSV[0], _temp_dPSV_g1, _temp_dPSV_g2]))
+ 1115 else:
+ 1116 eqs.append((1.0 / self.TSV) * (-PSV + self.PSV0 - (1.0 / self.RD) * (w / self.ws - 1)) - dPSV)
+ 1117 # Limitation of valve position Psv with limiter end
+ 1118
+ 1119 eqs.append((1.0 / self.TCH) * (-TM + PSV) - dTM) # (10)
+ 1120 eqs.append(
+ 1121 self.Rs * Id
+ 1122 - self.Xqpp * Iq
+ 1123 - ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * Edp
+ 1124 + ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * Si2q
+ 1125 + VG * np.sin(Angle_diff)
+ 1126 ) # (12)
+ 1127 eqs.append(
+ 1128 self.Rs * Iq
+ 1129 + self.Xdpp * Id
+ 1130 - ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * Eqp
+ 1131 - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * Si1d
+ 1132 + VG * np.cos(Angle_diff)
+ 1133 ) # (13)
+ 1134 eqs.append((Id * VG.T * np.sin(Angle_diff) + Iq * VG.T * np.cos(Angle_diff)) - PL2[0 : self.m] - sum1) # (14)
+ 1135 eqs.append((Id * VG.T * np.cos(Angle_diff) - Iq * VG.T * np.sin(Angle_diff)) - QL2[0 : self.m] - sum2) # (15)
+ 1136 eqs.append(-PL2[self.m : self.n] - sum3) # (16)
+ 1137 eqs.append(-QL2[self.m : self.n] - sum4) # (17)
+ 1138 eqs_flatten = [item for sublist in eqs for item in sublist]
+ 1139
+ 1140 f[:] = eqs_flatten
+ 1141 return f
+ 1142
+ 1143 def u_exact(self, t):
+ 1144 r"""
+ 1145 Returns the initial conditions at time :math:`t=0`.
+ 1146
+ 1147 Parameters
+ 1148 ----------
+ 1149 t : float
+ 1150 Time of the initial conditions.
+ 1151
+ 1152 Returns
+ 1153 -------
+ 1154 me : dtype_u
+ 1155 Initial conditions.
+ 1156 """
+ 1157 assert t == 0, 'ERROR: u_exact only valid for t=0'
+ 1158
+ 1159 me = self.dtype_u(self.init)
+ 1160 me[0 : self.m] = self.Eqp0
+ 1161 me[self.m : 2 * self.m] = self.Si1d0
+ 1162 me[2 * self.m : 3 * self.m] = self.Edp0
+ 1163 me[3 * self.m : 4 * self.m] = self.Si2q0
+ 1164 me[4 * self.m : 5 * self.m] = self.D0
+ 1165 me[5 * self.m : 6 * self.m] = self.ws_vector
+ 1166 me[6 * self.m : 7 * self.m] = self.Efd0
+ 1167 me[7 * self.m : 8 * self.m] = self.RF0
+ 1168 me[8 * self.m : 9 * self.m] = self.VR0
+ 1169 me[9 * self.m : 10 * self.m] = self.TM0
+ 1170 me[10 * self.m : 11 * self.m] = self.PSV0
+ 1171 me[11 * self.m : 11 * self.m + self.m] = self.Id0
+ 1172 me[11 * self.m + self.m : 11 * self.m + 2 * self.m] = self.Iq0
+ 1173 me[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n] = self.V0
+ 1174 me[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n] = self.TH0
+ 1175 return me
+ 1176
+ 1177 def get_switching_info(self, u, t, du=None):
+ 1178 r"""
+ 1179 Provides information about the state function of the problem. When the state function changes its sign,
+ 1180 typically an event occurs. So the check for an event should be done in the way that the state function
+ 1181 is checked for a sign change. If this is the case, the intermediate value theorem states a root in this
+ 1182 step.
+ 1183
+ 1184 The state function for this problem is given by
+ 1185
+ 1186 .. math::
+ 1187 h(P_{SV,1}(t)) = P_{SV,1}(t) - P_{SV,1, max}
+ 1188
+ 1189 for :math:`P_{SV,1,max}=1.0`.
+ 1190
+ 1191 Parameters
+ 1192 ----------
+ 1193 u : dtype_u
+ 1194 Current values of the numerical solution at time :math:`t`.
+ 1195 t : float
+ 1196 Current time of the numerical solution.
+ 1197
+ 1198 Returns
+ 1199 -------
+ 1200 switch_detected : bool
+ 1201 Indicates whether a discrete event is found or not.
+ 1202 m_guess : int
+ 1203 The index before the sign changes.
+ 1204 state_function : list
+ 1205 Defines the values of the state function at collocation nodes where it changes the sign.
+ 1206 """
+ 1207
+ 1208 switch_detected = False
+ 1209 m_guess = -100
+ 1210 for m in range(1, len(u)):
+ 1211 h_prev_node = u[m - 1][10 * self.m] - self.psv_max
+ 1212 h_curr_node = u[m][10 * self.m] - self.psv_max
+ 1213 if h_prev_node < 0 and h_curr_node >= 0:
+ 1214 switch_detected = True
+ 1215 m_guess = m - 1
+ 1216 break
+ 1217
+ 1218 state_function = [u[m][10 * self.m] - self.psv_max for m in range(len(u))]
+ 1219 return switch_detected, m_guess, state_function
+ 1220
+ 1221 def count_switches(self):
+ 1222 """
+ 1223 Setter to update the number of switches if one is found.
+ 1224 """
+ 1225 self.nswitches += 1
+
+