pylon/dyn.py

   1 #------------------------------------------------------------------------------
   2 # Copyright (C) 2007-2010 Richard Lincoln
   3 #
   4 # Licensed under the Apache License, Version 2.0 (the "License");
   5 # you may not use this file except in compliance with the License.
   6 # You may obtain a copy of the License at
   7 #
   8 #     http://www.apache.org/licenses/LICENSE-2.0
   9 #
  10 # Unless required by applicable law or agreed to in writing, software
  11 # distributed under the License is distributed on an "AS IS" BASIS,
  12 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13 # See the License for the specific language governing permissions and
  14 # limitations under the License.
  15 #------------------------------------------------------------------------------
  16
  17 """ Defines classes for dynamic simulation.
  18
  19 Based on MatDyn by Stijn Cole, developed at Katholieke Universiteit Leuven.
  20 See U{http://www.esat.kuleuven.be/electa/teaching/matdyn/} for more info.
  21 """
  22
  23 #------------------------------------------------------------------------------
  24 #  Imports:
  25 #------------------------------------------------------------------------------
  26
  27 import math
  28 import logging
  29
  30 from time import time
  31
  32 from numpy import \
  33     array, zeros, ones, exp, conj, pi, angle, abs, sin, cos, c_, r_, \
  34     flatnonzero, finfo
  35
  36 from scipy.sparse.linalg import spsolve, splu
  37
  38 from pylon import NewtonPF
  39
  40 from util import _Named, _Serializable
  41
  42 #------------------------------------------------------------------------------
  43 #  Logging:
  44 #------------------------------------------------------------------------------
  45
  46 logger = logging.getLogger(__name__)
  47
  48 #------------------------------------------------------------------------------
  49 #  Constants:
  50 #------------------------------------------------------------------------------
  51
  52 EPS = finfo(float).eps
  53
  54 CLASSICAL = "classical"
  55 FOURTH_ORDER = "fourth_order"
  56
  57 CONST_EXCITATION = "constant excitation"
  58 IEEE_DC1A = "IEEE DC1A"
  59
  60 CONST_POWER = "constant power"
  61 GENERAL_IEEE = "IEEE general speed-governing system"
  62
  63 BUS_CHANGE = "bus change"
  64 BRANCH_CHANGE = "branch change"
  65
  66 #------------------------------------------------------------------------------
  67 #  "DynamicCase" class:
  68 #------------------------------------------------------------------------------
  69
  70 class DynamicCase(_Named, _Serializable):
  71     """ Defines a dynamic simulation case.
  72     """
  73
  74     def __init__(self, case, frequency=50.0, stepsize=0.01, stoptime=1):
  75         """ Constructs a new DynamicCase instance.
  76         """
  77         #: Power flow case.
  78         self.case = case
  79
  80         #: Dynamic generators.
  81         self.dyn_generators = []
  82
  83         #: Generator exciters.
  84         self.exciters = []
  85
  86         #: Generator governors.
  87         self.governors = []
  88
  89         #: Network frequency (Hz).
  90         self.freq = frequency
  91
  92         #: Stepsize of the integration algorithm (s).
  93         self.stepsize = stepsize
  94
  95         #: Stoptime of the simulation (s).
  96         self.stoptime = stoptime
  97
  98
  99     def getAugYbus(self, U0, gbus):
 100         """ Based on AugYbus.m from MatDyn by Stijn Cole, developed at
 101         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 102         electa/teaching/matdyn/} for more information.
 103
 104         @rtype: csr_matrix
 105         @return: The augmented bus admittance matrix.
 106         """
 107         j = 0 + 1j
 108         buses = self.case.connected_buses
 109         nb = len(buses)
 110
 111         Ybus, _, _ = self.case.getYbus()
 112
 113         # Steady-state bus voltages.
 114
 115         # Calculate equivalent load admittance
 116         Sd = array([self.case.s_demand(bus) for bus in buses])
 117         Yd = conj(Sd) / abs(U0)**2
 118
 119         xd_tr = array([g.xd_tr for g in self.dyn_generators])
 120
 121         # Calculate equivalent generator admittance.
 122         Yg = zeros(nb)
 123         Yg[gbus] = 1 / (j * xd_tr)
 124
 125         # Add equivalent load and generator admittance to Ybus matrix
 126         for i in range(nb):
 127             Ybus[i, i] = Ybus[i, i] + Yg[i] + Yd[i]
 128
 129         return Ybus
 130
 131
 132     def generatorInit(self, U0):
 133         """ Based on GeneratorInit.m from MatDyn by Stijn Cole, developed at
 134         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 135         electa/teaching/matdyn/} for more information.
 136
 137         @rtype: tuple
 138         @return: Initial generator conditions.
 139         """
 140         j = 0 + 1j
 141         generators = self.dyn_generators
 142
 143         Efd0 = zeros(len(generators))
 144         Xgen0 = zeros((len(generators), 4))
 145
 146         typ1 = [g._i for g in generators if g.model == CLASSICAL]
 147         typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]
 148
 149         # Generator type 1: classical model
 150         x_tr = array([g.x_tr for g in generators])
 151
 152         omega0 = ones(len(typ1)) * 2 * pi * self.freq
 153
 154         # Initial machine armature currents.
 155         Sg = array([g.p + j * g.q for g in generators])
 156         Ia0 = conj(Sg[typ1]) / conj(U0) / self.base_mva
 157
 158         # Initial Steady-state internal EMF.
 159         Eq_tr0 = U0[typ1] + j * x_tr * Ia0
 160         delta0 = angle(Eq_tr0)
 161         Eq_tr0 = abs(Eq_tr0)
 162
 163         Xgen0[typ1, :] = c_[delta0, omega0, Eq_tr0]
 164
 165         # Generator type 2: 4th order model
 166         xd = array([g.xd for g in generators])
 167         xq = array([g.xq for g in generators])
 168         xd_tr = array([g.xd_tr for g in generators])
 169         xq_tr = array([g.xq_tr for g in generators])
 170
 171         omega0 = ones(len(typ2)) * 2 * pi * self.freq
 172
 173         # Initial machine armature currents.
 174         Ia0 = conj(Sg[typ2]) / conj(U0[typ2]) / self.base_mva
 175         phi0 = angle(Ia0)
 176
 177         # Initial Steady-state internal EMF.
 178         Eq0 = U0[typ2] + j * xq * Ia0
 179         delta0 = angle(Eq0)
 180
 181         # Machine currents in dq frame.
 182         Id0 = -abs(Ia0) * sin(delta0 - phi0)
 183         Iq0 =  abs(Ia0) * cos(delta0 - phi0)
 184
 185         # Field voltage.
 186         Efd0[typ2] = abs(Eq0) - (xd - xq) * Id0
 187
 188         # Initial Transient internal EMF.
 189         Eq_tr0 = Efd0[typ2] + (xd - xd_tr) * Id0
 190         Ed_tr0 = -(xq - xq_tr) * Iq0
 191
 192         Xgen0[typ2, :] = c_[delta0, omega0, Eq_tr0, Ed_tr0]
 193
 194         # Generator type 3:
 195
 196         # Generator type 4:
 197
 198         return Efd0, Xgen0
 199
 200
 201     def exciterInit(self, Xexc, Vexc):
 202         """ Based on ExciterInit.m from MatDyn by Stijn Cole, developed at
 203         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 204         electa/teaching/matdyn/} for more information.
 205
 206         @rtype: tuple
 207         @return: Exciter initial conditions.
 208         """
 209         exciters = self.exciters
 210
 211         Xexc0 = zeros(Xexc.shape)
 212         Pexc0 = zeros(len(exciters))
 213
 214         typ1 = [e.generator._i for e in exciters if e.model ==CONST_EXCITATION]
 215         typ2 = [e.generator._i for e in exciters if e.model == IEEE_DC1A]
 216
 217         # Exciter type 1: constant excitation
 218         Efd0 = Xexc[typ1, 0]
 219         Xexc0[typ1, 0] = Efd0
 220
 221         # Exciter type 2: IEEE DC1A
 222         Efd0 = Xexc[typ2, 0]
 223         Ka = array([e.ka for e in exciters])
 224         Ta = array([e.ta for e in exciters])
 225         Ke = array([e.ke for e in exciters])
 226         Te = array([e.te for e in exciters])
 227         Kf = array([e.kf for e in exciters])
 228         Tf = array([e.tf for e in exciters])
 229         Aex = array([e.aex for e in exciters])
 230         Bex = array([e.bex for e in exciters])
 231         Ur_min = array([e.ur_min for e in exciters])
 232         Ur_max = array([e.ur_max for e in exciters])
 233
 234         U = Vexc[typ2, 0]
 235
 236         Uf = zeros(len(typ2))
 237         Ux = Aex * exp(Bex * Efd0)
 238         Ur = Ux + Ke * Efd0
 239         Uref2 = U + (Ux + Ke * Efd0) / Ka - U
 240         Uref = U
 241
 242         Xexc0[typ2, :] = c_[Efd0, Uf, Ur]
 243         Pexc0[typ2, :] = c_[Ka, Ta, Ke, Te, Kf, Tf, Aex, Bex,
 244                             Ur_min, Ur_max, Uref, Uref2]
 245
 246         # Exciter type 3:
 247
 248         # Exciter type 4:
 249
 250         return Xexc0, Pexc0
 251
 252
 253     def governorInit(self, Xgov, Vgov):
 254         """ Based on GovernorInit.m from MatDyn by Stijn Cole, developed at
 255         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 256         electa/teaching/matdyn/} for more information.
 257
 258         @rtype: tuple
 259         @return: Initial governor conditions.
 260         """
 261         governors = self.governors
 262
 263         Xgov0 = zeros(Xgov.shape)
 264         Pgov0 = zeros(len(governors))
 265
 266         typ1 = [g.generator._i for g in governors if g.model == CONST_POWER]
 267         typ2 = [g.generator._i for g in governors if g.model == GENERAL_IEEE]
 268
 269         # Governor type 1: constant power
 270         Pm0 = Xgov[typ1, 0]
 271         Xgov0[typ1, 0] = Pm0
 272
 273         # Governor type 2: IEEE general speed-governing system
 274         Pm0 = Xgov[typ2, 0]
 275
 276         K = array([g.k for g in governors])
 277         T1 = array([g.t1 for g in governors])
 278         T2 = array([g.t2 for g in governors])
 279         T3 = array([g.t3 for g in governors])
 280         Pup = array([g.p_up for g in governors])
 281         Pdown = array([g.p_down for g in governors])
 282         Pmax = array([g.p_max for g in governors])
 283         Pmin = array([g.p_min for g in governors])
 284
 285         omega0 = Vgov[typ2, 0]
 286
 287         zz0 = Pm0
 288         PP0 = Pm0
 289
 290         P0 = K * (2 * pi * self.freq - omega0)
 291         xx0 = T1 * (1 - T2 / T1) * (2 * pi * self.freq - omega0)
 292
 293         Xgov0[typ2, :] = c_[Pm0, P0, xx0, zz0]
 294         Pgov0[typ2, :] = c_[K, T1, T2, T3, Pup, Pdown, Pmax, Pmin, PP0]
 295
 296         # Governor type 3:
 297
 298         # Governor type 4:
 299
 300         return Xgov0, Pgov0
 301
 302
 303     def machineCurrents(self, Xg, U):
 304         """ Based on MachineCurrents.m from MatDyn by Stijn Cole, developed at
 305         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 306         electa/teaching/matdyn/} for more information.
 307
 308         @param Xg: Generator state variables.
 309         @param U: Generator voltages.
 310
 311         @rtype: tuple
 312         @return: Currents and electric power of generators.
 313         """
 314         generators = self.dyn_generators
 315
 316         # Initialise.
 317         ng = len(generators)
 318         Id = zeros(ng)
 319         Iq = zeros(ng)
 320         Pe = zeros(ng)
 321
 322         typ1 = [g._i for g in generators if g.model == CLASSICAL]
 323         typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]
 324
 325         # Generator type 1: classical model
 326         delta = Xg[typ1, 0]
 327         Eq_tr = Xg[typ1, 2]
 328
 329         xd = array([g.xd for g in generators])
 330
 331         Pe[typ1] = \
 332             1 / xd * abs(U[typ1]) * abs(Eq_tr) * sin(delta - angle(U[typ1]))
 333
 334         # Generator type 2: 4th order model
 335         delta = Xg[typ1, 0]
 336         Eq_tr = Xg[typ1, 2]
 337         Ed_tr = Xg[typ1, 3]
 338
 339         xd_tr = array([g.xd_tr for g in generators])
 340         xq_tr = array([g.xq_tr for g in generators])
 341
 342         theta = angle(U)
 343
 344         # Transform U to rotor frame of reference.
 345         vd = -abs(U[typ2]) * sin(delta - theta[typ2])
 346         vq =  abs(U[typ2]) * cos(delta - theta[typ2])
 347
 348         Id[typ2] =  (vq - Eq_tr) / xd_tr
 349         Iq[typ2] = -(vd - Ed_tr) / xq_tr
 350
 351         Pe[typ2] = \
 352             Eq_tr * Iq[typ2] + Ed_tr * Id[typ2] + \
 353             (xd_tr - xq_tr) * Id[typ2] * Iq[typ2]
 354
 355         return Id, Iq, Pe
 356
 357
 358     def solveNetwork(self, Xgen, augYbus_solver, gbus):
 359         """ Based on SolveNetwork.m from MatDyn by Stijn Cole, developed at
 360         Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
 361         electa/teaching/matdyn/} for more information.
 362
 363         @rtype: array
 364         @return: Bus voltages.
 365         """
 366         generators = self.dyn_generators
 367         j = 0 + 1j
 368
 369         ng = len(gbus)
 370         Igen = zeros(ng)
 371
 372         s = len(augYbus_solver)
 373         Ig = zeros(s)
 374
 375         # Define generator types.
 376         typ1 = [g._i for g in generators if g.model == CLASSICAL]
 377         typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]
 378
 379         # Generator type 1: classical model
 380         delta = Xgen[typ1, 0]
 381         Eq_tr = Xgen[typ1, 2]
 382
 383         xd_tr = array([g.xd_tr for g in generators])[typ1]
 384
 385         # Calculate generator currents
 386         Igen[typ1] = (Eq_tr * exp(j * delta)) / (j * xd_tr)
 387
 388         # Generator type 2: 4th order model
 389         delta = Xgen[typ2, 0]
 390         Eq_tr = Xgen[typ2, 2]
 391         Ed_tr = Xgen[typ2, 3]
 392
 393         xd_tr = array([g.xd_tr for g in generators])[typ2] # Pgen(type2,8)
 394
 395         # Calculate generator currents. (Padiyar, p.417.)
 396         Igen[typ2] = (Eq_tr + j * Ed_tr) * exp(j * delta) / (j * xd_tr)
 397
 398         # Calculations --------------------------------------------------------
 399
 400         # Generator currents
 401         Ig[gbus] = Igen
 402
 403         # Calculate network voltages: U = Y/Ig
 404         U = augYbus_solver.solve(Ig)
 405
 406         return U
 407
 408
 409     def exciter(self, Xexc, Pexc, Vexc):
 410         """ Exciter model.
 411
 412         Based on Exciter.m from MatDyn by Stijn Cole, developed at Katholieke
 413         Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
 414         matdyn/} for more information.
 415         """
 416         exciters = self.exciters
 417
 418         F = zeros(Xexc.shape)
 419
 420         typ1 = [e.generator._i for e in exciters if e.model ==CONST_EXCITATION]
 421         typ2 = [e.generator._i for e in exciters if e.model == IEEE_DC1A]
 422
 423         # Exciter type 1: constant excitation
 424         F[typ1, :] = 0.0
 425
 426         # Exciter type 2: IEEE DC1A
 427         Efd = Xexc[typ2, 0]
 428         Uf = Xexc[typ2, 1]
 429         Ur = Xexc[typ2, 2]
 430
 431         Ka = Pexc[typ2, 0]
 432         Ta = Pexc[typ2, 1]
 433         Ke = Pexc[typ2, 2]
 434         Te = Pexc[typ2, 3]
 435         Kf = Pexc[typ2, 4]
 436         Tf = Pexc[typ2, 5]
 437         Aex = Pexc[typ2, 6]
 438         Bex = Pexc[typ2, 7]
 439         Ur_min = Pexc[typ2, 8]
 440         Ur_max = Pexc[typ2, 9]
 441         Uref = Pexc[typ2, 10]
 442         Uref2 = Pexc[typ2, 11]
 443
 444         U = Vexc[typ2, 1]
 445
 446         Ux = Aex * exp(Bex * Efd)
 447         dUr = 1 / Ta * (Ka * (Uref - U + Uref2 - Uf) - Ur)
 448         dUf = 1 / Tf * (Kf / Te * (Ur - Ux - Ke * Efd) - Uf)
 449
 450         if sum(flatnonzero(Ur > Ur_max)) >= 1:
 451             Ur2 = Ur_max
 452         elif sum(flatnonzero(Ur < Ur_max)) >= 1:
 453             Ur2 = Ur_min
 454         else:
 455             Ur2 = Ur
 456
 457         dEfd = 1 / Te * (Ur2 - Ux - Ke * Efd)
 458         F[typ2, :] = c_[dEfd, dUf, dUr]
 459
 460         # Exciter type 3:
 461
 462         # Exciter type 4:
 463
 464         return F
 465
 466
 467     def governor(self, Xgov, Pgov, Vgov):
 468         """ Governor model.
 469
 470         Based on Governor.m from MatDyn by Stijn Cole, developed at Katholieke
 471         Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
 472         matdyn/} for more information.
 473         """
 474         governors = self.governors
 475         omegas = 2 * pi * self.freq
 476
 477         F = zeros(Xgov.shape)
 478
 479         typ1 = [g.generator._i for g in governors if g.model == CONST_POWER]
 480         typ2 = [g.generator._i for g in governors if g.model == GENERAL_IEEE]
 481
 482         # Governor type 1: constant power
 483         F[typ1, 0] = 0
 484
 485         # Governor type 2: IEEE general speed-governing system
 486         Pm = Xgov[typ2, 0]
 487         P = Xgov[typ2, 1]
 488         x = Xgov[typ2, 2]
 489         z = Xgov[typ2, 3]
 490
 491         K = Pgov[typ2, 0]
 492         T1 = Pgov[typ2, 1]
 493         T2 = Pgov[typ2, 2]
 494         T3 = Pgov[typ2, 3]
 495         Pup = Pgov[typ2, 4]
 496         Pdown = Pgov[typ2, 5]
 497         Pmax = Pgov[typ2, 6]
 498         Pmin = Pgov[typ2, 7]
 499         P0 = Pgov[typ2, 8]
 500
 501         omega = Vgov[typ2, 0]
 502
 503         dx = K * (-1 / T1 * x + (1 - T2 / T1) * (omega - omegas))
 504         dP = 1 / T1 * x + T2 / T1 * (omega - omegas)
 505
 506         y = 1 / T3 * (P0 - P - Pm)
 507
 508         y2 = y
 509
 510         if sum(flatnonzero(y > Pup)) >= 1:
 511             y2 = (1 - flatnonzero(y > Pup)) * y2 + flatnonzero(y > Pup) * Pup
 512         if sum(flatnonzero(y < Pdown)) >= 1:
 513             y2 = (1 - flatnonzero(y<Pdown)) * y2 + flatnonzero(y<Pdown) * Pdown
 514
 515         dz = y2
 516
 517         dPm = y2
 518
 519         if sum(flatnonzero(z > Pmax)) >= 1:
 520             dPm = (1 - flatnonzero(z > Pmax)) * dPm + flatnonzero(z > Pmax) * 0
 521         if sum(flatnonzero(z < Pmin)) >= 1:
 522             dPm = (1 - flatnonzero(z < Pmin)) * dPm + flatnonzero(z < Pmin) * 0
 523
 524         F[typ2, :] = c_[dPm, dP, dx, dz]
 525
 526         # Governor type 3:
 527
 528         # Governor type 4:
 529
 530         return F
 531
 532
 533     def generator(self, Xgen, Xexc, Xgov, Vgen):
 534         """ Generator model.
 535
 536         Based on Generator.m from MatDyn by Stijn Cole, developed at Katholieke
 537         Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
 538         matdyn/} for more information.
 539         """
 540         generators = self.dyn_generators
 541         omegas = 2 * pi * self.freq
 542
 543         F = zeros(Xgen.shape)
 544
 545         typ1 = [g._i for g in generators if g.model == CLASSICAL]
 546         typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]
 547
 548         # Generator type 1: classical model
 549         omega = Xgen[typ1, 1]
 550         Pm0 = Xgov[typ1, 0]
 551
 552         H = array([g.h for g in generators])[typ1]
 553         D = array([g.d for g in generators])[typ1]
 554
 555         Pe = Vgen[typ1, 2]
 556
 557         ddelta = omega = omegas
 558         domega = pi * self.freq / H * (-D * (omega - omegas) + Pm0 - Pe)
 559         dEq = zeros(len(typ1))
 560
 561         F[typ1, :] = c_[ddelta, domega, dEq]
 562
 563         # Generator type 2: 4th order model
 564         omega = Xgen[typ2, 1]
 565         Eq_tr = Xgen[typ2, 2]
 566         Ed_tr = Xgen[typ2, 3]
 567
 568         H = array([g.h for g in generators])
 569         D = array([g.d for g in generators])
 570         xd = array([g.xd for g in generators])
 571         xq = array([g.xq for g in generators])
 572         xd_tr = array([g.xd_tr for g in generators])
 573         xq_tr = array([g.xq_tr for g in generators])
 574         Td0_tr = array([g.td for g in generators])
 575         Tq0_tr = array([g.tq for g in generators])
 576
 577         Id = Vgen[typ2, 0]
 578         Iq = Vgen[typ2, 1]
 579         Pe = Vgen[typ2, 2]
 580
 581         Efd = Xexc[typ2, 0]
 582         Pm = Xgov[typ2, 0]
 583
 584         ddelta = omega - omegas
 585         domega = pi * self.freq / H * (-D * (omega - omegas) + Pm - Pe)
 586         dEq = 1 / Td0_tr * (Efd - Eq_tr + (xd - xd_tr) * Id)
 587         dEd = 1 / Tq0_tr * (-Ed_tr - (xq - xq_tr) * Iq)
 588
 589         F[typ2, :] = c_[ddelta, domega, dEq, dEd]
 590
 591         # Generator type 3:
 592
 593         # Generator type 4:
 594
 595         return F
 596
 597 #------------------------------------------------------------------------------
 598 #  "DynamicSolver" class:
 599 #------------------------------------------------------------------------------
 600
 601 class DynamicSolver(object):
 602     """ Defines a solver for dynamic simulation.
 603
 604     The adaptive step size methods start with minimal step size. It is of
 605     interest to increase minimum step size as it speeds up the calculations.
 606     Generally, tolerance must be increased as well, or the integration will
 607     fail.
 608
 609     Based on rundyn.m from MatDyn by Stijn Cole, developed at Katholieke
 610     Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
 611     matdyn/} for more information.
 612     """
 613
 614     def __init__(self, dyn_case, method=None, tol=1e-04,
 615                  minstep=1e-03, maxstep=1e02, verbose=True, plot=True):
 616
 617         #: Dynamic case.
 618         self.dyn_case = dyn_case
 619
 620         #: Integration method.
 621         self.method = ModifiedEuler() if method is None else method
 622
 623         #: Specify the tolerance of the error. This argument is only used for
 624         #: the Runge-Kutta Fehlberg and Higham and Hall methods.
 625         self.tol = tol
 626
 627         #: Sets the minimum step size. Only used by the adaptive step size
 628         #: algorithms: Runge-Kutta Fehlberg and Higham and Hall methods.
 629         self.minstep = minstep
 630
 631         #: Sets the maximal step size. Only used by the adaptive step size
 632         #: algorithms: Runge-Kutta Fehlberg and Higham and Hall methods.
 633         self.maxstep = maxstep
 634
 635         #: Print progress output?
 636         self.verbose = verbose
 637
 638         #: Draw plot?
 639         self.plot = plot
 640
 641
 642     def solve(self):
 643         """ Runs dynamic simulation.
 644
 645         @rtype: dict
 646         @return: Solution dictionary with the following keys:
 647                    - C{angles} - generator angles
 648                    - C{speeds} - generator speeds
 649                    - C{eq_tr} - q component of transient voltage behind
 650                      reactance
 651                    - C{ed_tr} - d component of transient voltage behind
 652                      reactance
 653                    - C{efd} - Excitation voltage
 654                    - C{pm} - mechanical power
 655                    - C{voltages} - bus voltages
 656                    - C{stepsize} - step size integration method
 657                    - C{errest} - estimation of integration error
 658                    - C{failed} - failed steps
 659                    - C{time} - time points
 660         """
 661         t0 = time()
 662         buses = self.dyn_case.buses
 663
 664         solution = NewtonPF(self.case).solve()
 665
 666         if not solution["converged"]:
 667             logger.error("Power flow did not converge. Exiting...")
 668             return {}
 669         elif self.verbose:
 670             logger.info("Power flow converged.")
 671
 672         # Construct augmented Ybus.
 673         if self.verbose:
 674             logger.info("Constructing augmented admittance matrix...")
 675
 676         gbus = [g.bus._i for g in self.dyn_generators]
 677         ng = len(gbus)
 678
 679         Um = array([bus.v_magnitude for bus in buses])
 680         Ua = array([bus.v_angle * (pi / 180.0) for bus in buses])
 681         U0 = Um * exp(1j * Ua)
 682         U00 = U0
 683
 684         augYbus = self.dyn_case.getAugYbus(U0, gbus)
 685         augYbus_solver = splu(augYbus)
 686
 687         # Calculate initial machine state.
 688         if self.verbose:
 689             logger.info("Calculating initial state...")
 690
 691         Efd0, Xgen0 = self.dyn_case.generatorInit(U0)
 692         omega0 = Xgen0[:, 1]
 693
 694         Id0, Iq0, Pe0 = self.dyn_case.machineCurrents(Xgen0, U0)
 695         Vgen0 = r_[Id0, Iq0, Pe0]
 696
 697         # Exciter initial conditions.
 698         Vexc0 = abs(U0[gbus])
 699         Xexc0, Pexc0 = self.dyn_case.exciterInit(Efd0, Vexc0)
 700
 701         # Governor initial conditions.
 702         Pm0 = Pe0
 703         Xgov0, Pgov0 = self.dyn_case.governorInit(Pm0, omega0)
 704         Vgov0 = omega0
 705
 706         # Check steady-state.
 707         Fexc0 = self.dyn_case.exciter(Xexc0, Pexc0, Vexc0)
 708         Fgov0 = self.dyn_case.governor(Xgov0, Pgov0, Vgov0)
 709         Fgen0 = self.dyn_case.generator(Xgen0, Xexc0, Xgov0, Vgen0)
 710
 711         # Check Generator Steady-state
 712         if sum(abs(Fgen0)) > 1e-06:
 713             logger.error("Generator not in steady-state. Exiting...")
 714             return {}
 715         # Check Exciter Steady-state
 716         if sum(abs(Fexc0)) > 1e-06:
 717             logger.error("Exciter not in steady-state. Exiting...")
 718             return {}
 719         # Check Governor Steady-state
 720         if sum(abs(Fgov0)) > 1e-06:
 721             logger.error("Governor not in steady-state. Exiting...")
 722             return {}
 723
 724         if self.verbose:
 725             logger.info("System in steady-state.")
 726
 727         # Initialization of main stability loop.
 728         t = -0.02 # simulate 0.02s without applying events
 729         erst = False
 730         failed = False
 731         eulerfailed = False
 732
 733         stoptime = self.dyn_case.stoptime
 734
 735         if (isinstance(self.method, RungeKuttaFehlberg) or
 736             isinstance(self.method, RungeKuttaHighamHall)):
 737             stepsize = self.minstep
 738         else:
 739             stepsize = self.dyn_case.stepsize
 740
 741         ev = 0
 742         eventhappened = False
 743         i = 0
 744
 745         # Allocate memory for variables.
 746         if self.verbose:
 747             logger.info("Allocating memory...")
 748
 749         chunk = 5000
 750         time = zeros(chunk)
 751         time[0, :] = t
 752         errest = zeros(chunk)
 753         errest[0, :] = erst
 754         stepsizes = zeros(chunk)
 755         stepsizes[0, :] = stepsize
 756
 757         # System variables
 758         voltages = zeros(chunk)
 759         voltages[0, :] = U0.H
 760
 761         # Generator
 762         angles = zeros((chunk, ng))
 763         angles[0, :] = Xgen0[:, 0] * 180.0 / pi
 764         speeds = zeros((chunk, ng))
 765         speeds[0, :] = Xgen0[:, 0] / 2 * pi * self.dyn_case.freq
 766         Eq_tr = zeros((chunk, ng))
 767         Eq_tr[0, :] = Xgen0[:, 2]
 768         Ed_tr = zeros((chunk, ng))
 769         Ed_tr[0, :] = Xgen0[:, 3]
 770
 771         # Exciter and governor
 772         Efd = zeros((chunk, ng))
 773         Efd[0, :] = Efd0[:, 0]
 774         PM = zeros((chunk, ng))
 775         PM[0, :] = Pm0[:, 0]
 776
 777         # Main stability loop.
 778         while t < stoptime + stepsize:
 779             i += 1
 780             if i % 45 == 0 and self.verbose:
 781                 logger.info("%6.2f%% completed." % t / stoptime * 100)
 782
 783             # Numerical Method.
 784             Xgen0, self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0, \
 785                 Vgov0, U0, t, newstepsize = self.method.solve(t, Xgen0,
 786                     self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0,
 787                     Vgov0, augYbus_solver, gbus, stepsize)
 788
 789 #            if self.method == MODIFIED_EULER:
 790 #                solver = ModifiedEuler(t, Xgen0, self.Pgen0, Vgen0, Xexc0,
 791 #                                       Pexc0, Vexc0, Xgov0, Pgov0, Vgov0,
 792 #                                       augYbus_solver, gbus, stepsize)
 793 #
 794 #                Xgen0, self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0,
 795 #                Vgov0, U0, t, newstepsize = solver.solve()
 796 #            elif self.method == RUNGE_KUTTA:
 797 #                pass
 798 #            elif self.method == RUNGE_KUTTA_FEHLBERG:
 799 #                pass
 800 #            elif self.method == HIGHAM_HALL:
 801 #                pass
 802 #            elif self.method == MODIFIED_EULER2:
 803 #                pass
 804 #            else:
 805 #                raise ValueError
 806
 807             if eulerfailed:
 808                 logger.info("No solution found. Exiting... ")
 809                 return {}
 810
 811             if failed:
 812                 t = t - stepsize
 813
 814             # End exactly at stop time.
 815             if t + newstepsize > stoptime:
 816                 newstepsize = stoptime - t
 817             elif stepsize < self.minstep:
 818                 logger.info("No solution found with minimum step size. Exiting... ")
 819                 return {}
 820
 821             # Allocate new memory chunk if matrices are full.
 822             if i > time.shape[0]:
 823                 time = zeros(chunk)
 824                 errest = zeros(chunk)
 825                 stepsize = zeros(chunk)
 826                 voltages = zeros(chunk)
 827                 angles = zeros((chunk, ng))
 828                 speeds = zeros((chunk, ng))
 829                 Eq_tr = zeros((chunk, ng))
 830                 Ed_tr = zeros((chunk, ng))
 831                 Efd = zeros((chunk, ng))
 832                 PM = zeros((chunk, ng))
 833
 834             # Save values.
 835             stepsizes[i, :] = stepsize
 836             errest[i, :] = erst
 837             time[i, :] = t
 838
 839             voltages[i, :] = U0
 840
 841             # Exciters
 842             Efd[i, :] = Xexc0[:, 0]
 843             # TODO: Set Efd to zero when using classical generator model.
 844
 845             # Governors
 846             PM[i, :] = Xgov0[:, 0]
 847
 848             # Generators
 849             angles[i, :] = Xgen0[:, 0] * 180.0 / pi
 850             speeds[i, :] = Xgen0[:, 1] * (2 * pi * self.dyn_case.freq)
 851             Eq_tr[i, :] = Xgen0[:, 2]
 852             Ed_tr[i, :] = Xgen0[:, 3]
 853
 854             # Adapt step size if event will occur in next step.
 855             if (len(self.events) > 0 and ev <= len(self.events) and
 856                 isinstance(self.method, RungeKuttaFehlberg) and
 857                 isinstance(self.method, RungeKutta)):
 858
 859                 if t + newstepsize >= self.events[ev].t:
 860                     if self.events[ev] - t < newstepsize:
 861                         newstepsize = self.events[ev].t - t
 862
 863             # Check for events.
 864             if len(self.events) > 0 and ev <= len(self.events):
 865                 for event in self.events:
 866                     if (abs(t - self.events[ev].t) > 10 * EPS or
 867                         ev > len(self.events)):
 868                         break
 869                     else:
 870                         eventhappened = True
 871
 872                     event.obj.set_attr(event.param, event.newval)
 873
 874                     ev += 1
 875
 876                 if eventhappened:
 877                     # Refactorise.
 878                     self.dyn_case.getAugYbus(U00, gbus)
 879                     U0 = self.dyn_case.solveNetwork(Xgen0, self.Pgen0,
 880                                                     augYbus_solver, gbus)
 881
 882                     Id0, Iq0, Pe0 = self.dyn_case.machineCurrents(Xgen0,
 883                                                                   self.Pgen0,
 884                                                                   U0[gbus])
 885                     Vgen0 = r_[Id0, Iq0, Pe0]
 886                     Vexc0 = abs(U0[gbus])
 887
 888                     # Decrease stepsize after event occured.
 889                     if (isinstance(self.method, RungeKuttaFehlberg) or
 890                         isinstance(self.method, RungeKuttaHighamHall)):
 891
 892                         newstepsize = self.minstepsize
 893
 894                     # If event occurs, save values at t- and t+.
 895                     i += 1
 896
 897                     # Save values
 898                     stepsize[i, :] = stepsize.T
 899                     errest[i, :] = erst.T
 900                     time[i, :] = t
 901
 902                     voltages[i, :] = U0.T
 903
 904                     # Exciters.
 905                     # Set Efd to zero when using classical generator model.
 906 #                    Efd[i, :] = Xexc0[:, 1] * (flatnonzero(genmodel > 1))
 907
 908                     # Governors.
 909                     PM[i, :] = Xgov0[:, 1]
 910
 911                     # Generators.
 912                     angles[i, :] = Xgen0[:, 0] * 180.0 / pi
 913                     speeds[i, :] = Xgen0[:, 1] / (2.0 * pi * self.freq)
 914                     Eq_tr[i, :] = Xgen0[:, 2]
 915                     Ed_tr[i, :] = Xgen0[:, 3]
 916
 917                     eventhappened = False
 918
 919             # Advance time
 920             stepsize = newstepsize
 921             t += stepsize
 922
 923         # End of main stability loop ------------------------------------------
 924
 925         # Output --------------------------------------------------------------
 926
 927         if self.verbose:
 928             logger.info("100%% completed")
 929             elapsed = time() - t0
 930             logger.info("Simulation completed in %5.2f seconds." % elapsed)
 931
 932         # Save only the first i elements.
 933         angles = angles[0:i, :]
 934         speeds = speeds[0:i, :]
 935         Eq_tr = Eq_tr[0:i, :]
 936         Ed_tr = Ed_tr[0:i, :]
 937
 938         Efd = Efd[0:i, :]
 939         PM = PM[0:i, :]
 940
 941         voltages = voltages[0:i, :]
 942
 943         stepsize = stepsize[0:i, :]
 944         errest = errest[0:i, :]
 945         time = time[0:i, :]
 946
 947         if self.plot:
 948             raise NotImplementedError
 949
 950         return {}
 951
 952 #------------------------------------------------------------------------------
 953 #  "ModifiedEuler" class:
 954 #------------------------------------------------------------------------------
 955
 956 class ModifiedEuler(object):
 957     """ Modified Euler ODE solver.
 958
 959     Based on ModifiedEuler.m from MatDyn by Stijn Cole, developed at Katholieke
 960     Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
 961     matdyn/} for more information.
 962     """
 963
 964     def solve(self, t, Xgen0, Pgen, Vgen0, Xexc0, Pexc, Vexc0, Xgov0, Pgov,
 965             Vgov0, augYbus_solver, gbus, stepsize):
 966         case = self.dyn_case
 967
 968         # First Euler step ----------------------------------------------------
 969
 970         # Exciters.
 971         dFexc0 = case.exciter(Xexc0, Pexc, Vexc0)
 972         Xexc1 = Xexc0 + case.stepsize * dFexc0
 973
 974         # Governors.
 975         dFgov0 = case.governor(Xgov0, Pgov, Vgov0)
 976         Xgov1 = Xgov0 + case.stepsize * dFgov0
 977
 978         # Generators.
 979         dFgen0 = case.generator(Xgen0, Xexc1, Xgov1, Pgen,Vgen0)
 980         Xgen1 = Xgen0 + case.stepsize * dFgen0
 981
 982         # Calculate system voltages.
 983         U1 = case.solveNetwork(Xgen1, Pgen, augYbus_solver,gbus)
 984
 985         # Calculate machine currents and power.
 986         Id1, Iq1, Pe1 = case.machineCurrents(Xgen1, Pgen, U1[gbus])
 987
 988         # Update variables that have changed.
 989         Vexc1 = abs(U1[gbus])
 990         Vgen1 = r_[Id1, Iq1, Pe1]
 991         Vgov1 = Xgen1[:, 2]
 992
 993         # Second Euler step ---------------------------------------------------
 994
 995         # Exciters.
 996         dFexc1 = case.exciter(Xexc1, Pexc, Vexc1)
 997         Xexc2 = Xexc0 + case.stepsize / 2 * (dFexc0 + dFexc1)
 998
 999         # Governors.
1000         dFgov1 = case.governor(Xgov1, Pgov, Vgov1)
1001         Xgov2 = Xgov0 + case.stepsize / 2 * (dFgov0 + dFgov1)
1002
1003         # Generators.
1004         dFgen1 = case.generator(Xgen1, Xexc2, Xgov2, Pgen, Vgen1)
1005         Xgen2 = Xgen0 + case.stepsize / 2 * (dFgen0 + dFgen1)
1006
1007         # Calculate system voltages.
1008         U2 = case.solveNetwork(Xgen2, Pgen, augYbus_solver, gbus)
1009
1010         # Calculate machine currents and power.
1011         Id2, Iq2, Pe2 = case.machineCurrents(Xgen2, Pgen, U2[gbus])
1012
1013         # Update variables that have changed.
1014         Vgen2 = r_[Id2, Iq2, Pe2]
1015         Vexc2 = abs(U2[gbus])
1016         Vgov2 = Xgen2[:, 2]
1017
1018         return Xgen2, Pgen, Vgen2, Xexc2, Pexc, Vexc2, \
1019             Xgov2, Pgov, Vgov2, U2, t, stepsize
1020
1021 #------------------------------------------------------------------------------
1022 #  "RungeKutta" class:
1023 #------------------------------------------------------------------------------
1024
1025 class RungeKutta(object):
1026     """ Standard 4th order Runge-Kutta ODE solver.
1027
1028     Based on RundeKutta.m from MatDyn by Stijn Cole, developed at Katholieke
1029     Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/teaching/
1030     matdyn/} for more information.
1031     """
1032
1033     def __init__(self):
1034         #: Runge-Kutta coefficients.
1035         self._a = array([0.0, 0.0, 0.0, 0.0, 1.0/2.0, 0.0, 0.0, 0.0, 0.0,
1036                          1.0/2.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0])
1037         #: Runge-Kutta coefficients.
1038         self._b = array([1.0/6.0, 2.0/6.0, 2.0/6.0, 1.0/6.0])
1039         # self._c = array([0.0, 1.0/2.0, 1.0/2.0, 1.0]) # not used
1040
1041
1042     def solve(self, t, Xgen0, Pgen, Vgen0, Xexc0, Pexc, Vexc0, Xgov0, Pgov,
1043             Vgov0, augYbus_solver, gbus, stepsize):
1044
1045         Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1, Vgen1, Kgen1, _ = \
1046             self._k1()
1047
1048         Xexc2, Vexc2, Kexc2, Xgov2, Vgov2, Kgov2, Xgen2, Vgen2, Kgen2, _ = \
1049             self._k2(Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1, Vgen1,
1050                 Kgen1)
1051
1052         Xexc3, Vexc3, Kexc3, Xgov3, Vgov3, Kgov3, Xgen3, Vgen3, Kgen3, _ = \
1053             self._k3(Xexc2, Vexc2, Kexc2, Kexc1, Xgov2, Vgov2, Kgov2, Kgov1,
1054                 Xgen2, Vgen2, Kgen2, Kgen1)
1055
1056         Xexc4, Vexc4, _, Xgov4, Vgov4, _, Xgen4, Vgen4, _, U4 = \
1057             self._k4(Xexc3, Vexc3, Kexc3, Kexc2, Kexc1, Xgov3, Vgov3, Kgov3,
1058                 Kgov2, Kgov1, Xgen3, Vgen3, Kgen3, Kgen2, Kgen1)
1059
1060         return Xgen4, Vgen4, Xexc4, Vexc4, Xgov4, Vgov4, U4
1061
1062
1063     def _k1(self, Xexc0, Pexc, Vexc0, Xgov0, Pgov, Vgov0, Xgen0, Pgen, Vgen0,
1064             augYbus_solver, gbus):
1065         case = self.dyn_case
1066         a = self._a
1067
1068         # Exciters.
1069         Kexc1 = case.exciter(Xexc0, Pexc, Vexc0)
1070         Xexc1 = Xexc0 + case.stepsize * a[1, 0] * Kexc1
1071
1072         # Governors
1073         Kgov1 = case.governor(Xgov0, Pgov, Vgov0)
1074         Xgov1 = Xgov0 + case.stepsize * a[1, 0] * Kgov1
1075
1076         # Generators.
1077         Kgen1 = case.generator(Xgen0, Xexc1, Xgov1, Pgen, Vgen0)
1078         Xgen1 = Xgen0 + case.stepsize * a[1, 0] * Kgen1
1079
1080         # Calculate system voltages.
1081         U1 = case.solveNetwork(Xgen1, Pgen, augYbus_solver, gbus)
1082
1083         # Calculate machine currents and power.
1084         Id1, Iq1, Pe1 = case.machineCurrents(Xgen1, Pgen, U1[gbus])
1085
1086         # Update variables that have changed
1087         Vexc1 = abs(U1[gbus])
1088         Vgen1 = r_[Id1, Iq1, Pe1];
1089         Vgov1 = Xgen1[:, 1]
1090
1091         return Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1, Vgen1, Kgen1,U1
1092
1093
1094     def _k2(self, Xexc1, Xexc0, Pexc, Vexc1, Kexc1, Xgov1, Xgov0, Pgov, Vgov1,
1095             Kgov1, Xgen1, Xgen0, Pgen, Vgen1, Kgen1, augYbus_solver, gbus):
1096         case = self.dyn_case
1097         a = self._a
1098
1099         # Exciters.
1100         Kexc2 = case.exciter(Xexc1, Pexc, Vexc1)
1101         Xexc2 = Xexc0 + case.stepsize * (a[2, 0] * Kexc1 + a[2, 1] *Kexc2)
1102
1103         # Governors.
1104         Kgov2 = case.governor(Xgov1, Pgov, Vgov1)
1105         Xgov2 = Xgov0 + case.stepsize * (a[2, 0] * Kgov1 + a[2, 1] *Kgov2)
1106
1107         # Generators.
1108         Kgen2 = case.generator(Xgen1, Xexc2, Xgov2, Pgen, Vgen1)
1109         Xgen2 = Xgen0 + case.stepsize * (a[2, 0] * Kgen1 + a[2, 1] *Kgen2)
1110
1111         # Calculate system voltages.
1112         U2 = case.solveNetwork(Xgen2, Pgen, augYbus_solver,gbus)
1113
1114         # Calculate machine currents and power.
1115         Id2, Iq2, Pe2 = case.machineCurrents(Xgen2, Pgen, U2[gbus])
1116
1117         # Update variables that have changed
1118         Vexc2 = abs(U2[gbus])
1119         Vgen2 = r_[Id2, Iq2, Pe2]
1120         Vgov2 = Xgen2[:, 1]
1121
1122         return Xexc2, Vexc2, Kexc2, Xgov2, Vgov2, Kgov2, Xgen2, Vgen2, Kgen2,U2
1123
1124
1125     def _k3(self, Xexc2, Xexc0, Pexc, Vexc2, Kexc2, Kexc1, Xgov2, Xgov0, Pgov,
1126             Vgov2, Kgov2, Kgov1, Xgen2, Xgen0, Pgen, Vgen2, Kgen2, Kgen1,
1127             augYbus_solver, gbus):
1128         case = self.dyn_case
1129         a = self._a
1130
1131         # Exciters.
1132         Kexc3 = case.exciter(Xexc2, Pexc, Vexc2)
1133         Xexc3 = Xexc0 + case.stepsize * \
1134             (a[3, 0] * Kexc1 + a[3, 1] * Kexc2 + a[3, 2] * Kexc3)
1135
1136         # Governors.
1137         Kgov3 = case.governor(Xgov2, Pgov, Vgov2)
1138         Xgov3 = Xgov0 + case.stepsize * \
1139             (a[3,0] * Kgov1 + a[3, 1] * Kgov2 + a[3, 2] * Kgov3)
1140
1141         # Generators.
1142         Kgen3 = case.generator(Xgen2, Xexc3, Xgov3, Pgen, Vgen2)
1143         Xgen3 = Xgen0 + case.stepsize * \
1144             (a[3, 0] * Kgen1 + a[3, 1] * Kgen2 + a[3, 2] * Kgen3)
1145
1146         # Calculate system voltages.
1147         U3 = case.solveNetwork(Xgen3, Pgen, augYbus_solver,gbus)
1148
1149         # Calculate machine currents and power
1150         Id3, Iq3, Pe3 = case.machineCurrents(Xgen3, Pgen, U3[gbus])
1151
1152         # Update variables that have changed.
1153         Vexc3 = abs(U3[gbus])
1154         Vgen3 = r_[Id3, Iq3, Pe3]
1155         Vgov3 = Xgen3[:, 1]
1156
1157         return Xexc3, Vexc3, Kexc3, Xgov3, Vgov3, Kgov3, Xgen3, Vgen3, Kgen3,U3
1158
1159
1160     def _k4(self, Xexc3, Xexc0, Pexc, Vexc3, Kexc3, Kexc2, Kexc1, Xgov3, Xgov0,
1161             Pgov, Vgov3, Kgov3, Kgov2, Kgov1, Xgen3, Xgen0, Pgen, Vgen3, Kgen3,
1162             Kgen2, Kgen1, augYbus_solver, gbus):
1163         case = self.dyn_case
1164         b = self._b
1165
1166         # Exciters.
1167         Kexc4 = case.exciter(Xexc3, Pexc, Vexc3)
1168         Xexc4 = Xexc0 + case.stepsize * \
1169             (b[0] * Kexc1 + b[1] * Kexc2 + b[2] * Kexc3 + b[3] * Kexc4)
1170
1171         # Governors.
1172         Kgov4 = case.governor(Xgov3, Pgov, Vgov3)
1173         Xgov4 = Xgov0 + case.stepsize * \
1174             (b[0] * Kgov1 + b[1] * Kgov2 + b[2] * Kgov3 + b[3] * Kgov4)
1175
1176         # Generators.
1177         Kgen4 = case.generator(Xgen3, Xexc4, Xgov4, Pgen, Vgen3)
1178         Xgen4 = Xgen0 + case.stepsize * \
1179             (b[0] * Kgen1 + b[1] * Kgen2 + b[2] * Kgen3 + b[3] * Kgen4)
1180
1181         # Calculate system voltages.
1182         U4 = case.solveNetwork(Xgen4, Pgen, augYbus_solver,gbus)
1183
1184         # Calculate machine currents and power.
1185         Id4, Iq4, Pe4 = case.machineCurrents(Xgen4, Pgen, U4[gbus])
1186
1187         # Update variables that have changed
1188         Vexc4 = abs(U4[gbus])
1189         Vgen4 = r_[Id4, Iq4, Pe4]
1190         Vgov4 = Xgen4[:, 1]
1191
1192         return Xexc4, Vexc4, Kexc4, Xgov4, Vgov4, Kgov4, Xgen4, Vgen4, Kgen4,U4
1193
1194 #------------------------------------------------------------------------------
1195 #  "RungeKuttaFehlberg" class:
1196 #------------------------------------------------------------------------------
1197
1198 class RungeKuttaFehlberg(RungeKutta):
1199     """ Runge-Kutta Fehlberg ODE solver.
1200
1201     Based on RungeKuttaFehlberg.m from MatDyn by Stijn Cole, developed at
1202     Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/
1203     teaching/matdyn/} for more information.
1204     """
1205
1206     def __init__(self, t, Xgen0, Pgen, Vgen0, Xexc0, Pexc, Vexc0, Xgov0, Pgov,
1207             Vgov0, augYbus_solver, gbus, stepsize):
1208
1209         super(self, RungeKuttaFehlberg).__init__(t, Xgen0, Pgen, Vgen0, Xexc0,
1210             Pexc, Vexc0, Xgov0, Pgov, Vgov0, augYbus_solver, gbus, stepsize)
1211
1212         #: Runge-Kutta coefficients
1213         self._a = array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0/4.0, 0.0, 0.0, 0.0, 0.0,
1214                          3.0/32.0, 9.0/32.0, 0.0, 0.0, 0.0, 1932.0/2197.0,
1215                          -7200.0/2197.0, 7296.0/2197.0, 0.0, 0.0, 439.0/216.0,
1216                          -8.0, 3680.0/513.0, -845.0/4104.0, 0.0, -8.0/27.0,
1217                          2.0, -3544.0/2565.0, 1859.0/4104.0, -11.0/40.0])
1218         #: Runge-Kutta coefficients.
1219         self._b1 = array([25.0/216.0, 0.0, 1408.0/2565.0, 2197.0/4104.0,
1220                           -1.0/5.0, 0.0])
1221         #: Runge-Kutta coefficients.
1222         self._b2 = array([16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0,
1223                           -9.0/50.0, 2.0/55.0,])
1224         # c = array([0.0, 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0,])#not used
1225
1226
1227     def solve(self):
1228         case = self.dyn_case
1229         b2 = self._b2
1230
1231         accept = False
1232         facmax = 4
1233         failed = False
1234
1235         i = 0
1236         while accept == False:
1237             i += 1
1238
1239             Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1, Vgen1, Kgen1, _ =\
1240                 self._k1()
1241
1242             Xexc2, Vexc2, Kexc2, Xgov2, Vgov2, Kgov2, Xgen2, Vgen2, Kgen2, _ =\
1243                 self._k2(Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1,
1244                     Vgen1, Kgen1)
1245
1246             Xexc3, Vexc3, Kexc3, Xgov3, Vgov3, Kgov3, Xgen3, Vgen3, Kgen3, _ =\
1247                 self._k3(Xexc2, Vexc2, Kexc2, Kexc1, Xgov2, Vgov2, Kgov2,
1248                     Kgov1, Xgen2, Vgen2, Kgen2, Kgen1)
1249
1250             Xexc4, Vexc4, Kexc4, Xgov4, Vgov4, Kgov4, Xgen4, Vgen4, Kgen4, _ =\
1251                 self._k4(Xexc3, Vexc3, Kexc3, Kexc2, Kexc1, Xgov3, Vgov3,
1252                     Kgov3, Kgov2, Kgov1, Xgen3, Vgen3, Kgen3, Kgen2, Kgen1)
1253
1254             Xexc5, Vexc5, Kexc5, Xgov5, Vgov5, Kgov5, Xgen5, Vgen5, Kgen5, _ =\
1255                 self._k5(Xexc4, Vexc4, Kexc4, Kexc3, Kexc2, Kexc1, Xgov4,
1256                     Vgov4, Kgov4, Kgov3, Kgov2, Kgov1, Xgen4, Vgen4, Kgen4,
1257                     Kgen3, Kgen2, Kgen1)
1258
1259             Xexc6, Vexc6, Kexc6, Xgov6, Vgov6, Kgov6, Xgen6, Vgen6, Kgen6, U6=\
1260                 self._k6(Xexc5, Vexc5, Kexc5, Kexc4, Kexc3, Kexc2, Kexc1,
1261                     Xgov5, Vgov5, Kgov5, Kgov4, Kgov3, Kgov2, Kgov1, Xgen5,
1262                     Vgen5, Kgen5, Kgen4, Kgen3, Kgen2, Kgen1)
1263
1264             # Second, higher order solution.
1265             Xexc62 = self.Xexc0 + case.stepsize * \
1266                 (b2[0] * Kexc1 + b2[1] * Kexc2 + b2[2] * Kexc3 + \
1267                 b2[3] * Kexc4 + b2[4] * Kexc5 + b2[5] * Kexc6)
1268             Xgov62 = self.Xgov0 + case.stepsize * \
1269                 (b2[0] * Kgov1 + b2[1] * Kgov2 + b2[2] * Kgov3 + \
1270                 b2[3] * Kgov4 + b2[4] * Kgov5 + b2[5] * Kgov6)
1271             Xgen62 = self.Xgen0 + case.stepsize * \
1272                 (b2[0] * Kgen1 + b2[1] * Kgen2 + b2[2] * Kgen3 + \
1273                 b2[3] * Kgen4 + b2[4] * Kgen5 + b2[5] * Kgen6)
1274
1275             # Error estimate.
1276             Xexc = abs((Xexc62 - Xexc6).T)
1277             Xgov = abs((Xgov62 - Xgov6).T)
1278             Xgen = abs((Xgen62 - Xgen6).T)
1279             errest = max( r_[max(max(Xexc)), max(max(Xgov)), max(max(Xgen))] )
1280
1281             if errest < EPS:
1282                 errest = EPS
1283
1284             q = 0.84 * (self.tol / errest)^(1.0/4.0)
1285
1286             if errest < self.tol:
1287                 accept = True
1288                 U0 = U6
1289                 Vgen0 = Vgen6;
1290                 Vgov0 = Vgov6
1291                 Vexc0 = Vexc6
1292
1293                 Xgen0 = Xgen6
1294                 Xexc0 = Xexc6
1295                 Xgov0 = Xgov6
1296
1297                 Pgen0 = self.Pgen
1298                 Pexc0 = self.Pexc
1299                 Pgov0 = self.Pgov
1300
1301 #                t = t0
1302             else:
1303                 failed += 1
1304                 facmax = 1
1305
1306 #                t = t0
1307                 Pgen0 = self.Pgen
1308                 Pexc0 = self.Pexc
1309                 Pgov0 = self.Pgov
1310
1311                 stepsize = min(max(q, 0.1), facmax) * case.stepsize
1312                 return
1313
1314             stepsize = min(max(q, 0.1), facmax) * case.stepsize
1315
1316             if stepsize > self.maxstepsize:
1317                 stepsize = self.maxstepsize
1318
1319         return Xgen0, Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0, Vgov0, \
1320             U0, errest, failed, self.t, stepsize
1321
1322
1323     def _k4(self, Xexc3, Vexc3, Kexc3, Kexc2, Kexc1, Xgov3, Vgov3, Kgov3,
1324             Kgov2, Kgov1, Xgen3, Vgen3, Kgen3, Kgen2, Kgen1):
1325         # Overrides the standard Runge-Kutta K4 method.
1326         case = self.dyn_case
1327         a = self._a
1328
1329         # Exciters.
1330         Kexc4 = case.exciter(Xexc3, self.Pexc, Vexc3)
1331         Xexc4 = self.Xexc0 + case.stepsize * \
1332             (a[4, 0] * Kexc1 + a[4, 1] * Kexc2 + a[4, 2] * Kexc3 +
1333              a[4, 3] * Kexc4)
1334
1335         # Governors.
1336         Kgov4 = case.governor(Xgov3, self.Pgov, Vgov3)
1337         Xgov4 = self.Xgov0 + case.stepsize * \
1338             (a[4, 0] * Kgov1 + a[4, 1] * Kgov2 + a[4, 2] * Kgov3 +
1339              a[4, 3] * Kgov4)
1340
1341         # Generators.
1342         Kgen4 = case.generator(Xgen3, Xexc4, Xgov4, self.Pgen, Vgen3)
1343         Xgen4 = self.Xgen0 + case.stepsize * \
1344             (a[4, 0] * Kgen1 + a[4, 1] * Kgen2 + a[4, 2] * Kgen3 +
1345              a[4, 3] * Kgen4)
1346
1347         # Calculate system voltages.
1348         U4 = case.solveNetwork(Xgen4, self.Pgen, self.augYbus_solver,self.gbus)
1349
1350         # Calculate machine currents and power.
1351         Id4, Iq4, Pe4 = case.machineCurrents(Xgen4, self.Pgen, U4[self.gbus])
1352
1353         # Update variables that have changed
1354         Vexc4 = abs(U4[self.gbus])
1355         Vgen4 = r_[Id4, Iq4, Pe4]
1356         Vgov4 = Xgen4[:, 1]
1357
1358         return Xexc4, Vexc4, Kexc4, Xgov4, Vgov4, Kgov4, Xgen4, Vgen4, Kgen4,U4
1359
1360
1361     def _k5(self, Xexc4, Vexc4, Kexc4, Kexc3, Kexc2, Kexc1, Xgov4, Vgov4,Kgov4,
1362             Kgov3, Kgov2, Kgov1, Xgen4, Vgen4, Kgen4, Kgen3, Kgen2, Kgen1):
1363         case = self.dyn_case
1364         a = self._a
1365
1366         # Exciters.
1367         Kexc5 = case.exciter(Xexc4, self.Pexc, Vexc4)
1368         Xexc5 = self.Xexc0 + case.stepsize * \
1369             (a[5, 0] * Kexc1 + a[5, 1] * Kexc2 + a[5, 2] * Kexc3 + \
1370             a[5, 3] * Kexc4 + a[5, 4] * Kexc5)
1371
1372         # Governors.
1373         Kgov5 = case.governor(Xgov4, self.Pgov, Vgov4)
1374         Xgov5 = self.Xgov0 + case.stepsize * \
1375             (a[5, 0] * Kgov1 + a[5, 1] * Kgov2 + a[5, 2] * Kgov3 + \
1376             a[5, 3] * Kgov4 + a[5, 4] * Kgov5)
1377
1378         # Generators.
1379         Kgen5 = case.generator(Xgen4, Xexc5, Xgov5, self.Pgen, Vgen4)
1380         Xgen5 = self.Xgen0 + case.stepsize * \
1381             (a[5, 0] * Kgen1 + a[5, 1] * Kgen2 + a[5, 2] * Kgen3 + \
1382             a[5, 3] * Kgen4 + a[5, 4] * Kgen5)
1383
1384         # Calculate system voltages.
1385         U5 = case.solveNetwork(Xgen5, self.Pgen, self.augYbus_solver,
1386                                self.gbus)
1387
1388         # Calculate machine currents and power.
1389         Id5, Iq5, Pe5 = case.machineCurrents(Xgen5, self.Pgen,
1390                                              U5[self.gbus])
1391
1392         # Update variables that have changed.
1393         Vexc5 = abs(U5[self.gbus])
1394         Vgen5 = r_[Id5, Iq5, Pe5]
1395         Vgov5 = Xgen5[:, 1]
1396
1397         return Xexc5, Vexc5, Kexc5, Xgov5, Vgov5, Kgov5, Xgen5, Vgen5, Kgen5,U5
1398
1399
1400     def _k6(self, Xexc5, Vexc5, Kexc5, Kexc4, Kexc3, Kexc2, Kexc1, Xgov5,
1401             Vgov5, Kgov5, Kgov4, Kgov3, Kgov2, Kgov1, Xgen5, Vgen5, Kgen5,
1402             Kgen4, Kgen3, Kgen2, Kgen1):
1403         case = self.dyn_case
1404         b1 = self._b1
1405
1406         # Exciters.
1407         Kexc6 = case.exciter(Xexc5, self.Pexc, Vexc5)
1408         Xexc6 = self.Xexc0 + case.stepsize * \
1409             (b1[0] * Kexc1 + b1[1] * Kexc2 + b1[2] * Kexc3 + \
1410             b1[3] * Kexc4 + b1[4] * Kexc5 + b1[5] * Kexc6)
1411
1412         # Governors.
1413         Kgov6 = case.governor(Xgov5, self.Pgov, Vgov5)
1414         Xgov6 = self.Xgov0 + case.stepsize * \
1415             (b1[0] * Kgov1 + b1[1] * Kgov2 + b1[2] * Kgov3 + \
1416             b1[3] * Kgov4 + b1[4] * Kgov5 + b1[5] * Kgov6)
1417
1418         # Generators.
1419         Kgen6 = case.generator(Xgen5, Xexc6, Xgov6, self.Pgen, Vgen5)
1420         Xgen6 = self.Xgen0 + case.stepsize * \
1421             (b1[0] * Kgen1 + b1[1] * Kgen2 + b1[2] * Kgen3 + \
1422             b1[3] * Kgen4 + b1[4] * Kgen5 + b1[5] * Kgen6)
1423
1424         # Calculate system voltages.
1425         U6 = case.solveNetwork(Xgen6, self.Pgen, self.augYbus_solver,
1426                                self.gbus)
1427
1428         # Calculate machine currents and power.
1429         Id6, Iq6, Pe6 = case.machineCurrents(Xgen6, self.Pgen,
1430                                              U6[self.gbus])
1431
1432         # Update variables that have changed.
1433         Vexc6 = abs(U6[self.gbus])
1434         Vgen6 = r_[Id6, Iq6, Pe6]
1435         Vgov6 = Xgen6[:, 1]
1436
1437         return Xexc6, Vexc6, Kexc6, Xgov6, Vgov6, Kgov6, Xgen6, Vgen6, Kgen6,U6
1438
1439 #            # K1 --------------------------------------------------------------
1440 #
1441 #            # Exciters.
1442 #            Kexc1 = case.exciter(self.Xexc0, self.Pexc, self.Vexc0)
1443 #            Xexc1 = self.Xexc0 + case.stepsize * a[1, 0] * Kexc1
1444 #
1445 #            # Governors.
1446 #            Kgov1 = case.governor(self.Xgov0, self.Pgov, self.Vgov0)
1447 #            Xgov1 = self.Xgov0 + case.stepsize * a[1, 0] * Kgov1
1448 #
1449 #            # Generators.
1450 #            Kgen1 = case.generator(self.Xgen0, Xexc1, Xgov1, self.Pgen,
1451 #                                   self.Vgen0)
1452 #            Xgen1 = self.Xgen0 + case.stepsize * a[1, 0] * Kgen1
1453 #
1454 #            # Calculate system voltages.
1455 #            U1 = case.solveNetwork(Xgen1, self.Pgen, self.augYbus_solver,
1456 #                                   self.gbus)
1457 #
1458 #            # Calculate machine currents and power.
1459 #            Id1, Iq1, Pe1 = case.machineCurrents(Xgen1, self.Pgen,
1460 #                                                 U1[self.gbus])
1461 #
1462 #            # Update variables that have changed.
1463 #            Vexc1 = abs(U1[self.gbus])
1464 #            Vgen1 = r_[Id1, Iq1, Pe1]
1465 #            Vgov1 = Xgen1[:, 1]
1466 #
1467 #            # K2 --------------------------------------------------------------
1468 #
1469 #            # Exciters.
1470 #            Kexc2 = case.exciter(Xexc1, self.Pexc, Vexc1)
1471 #            Xexc2 = self.Xexc0 + case.stepsize * \
1472 #                (a[2, 0] * Kexc1 + a[2, 1] * Kexc2)
1473 #
1474 #            # Governors.
1475 #            Kgov2 = case.governor(Xgov1, self.Pgov, Vgov1)
1476 #            Xgov2 = self.Xgov0 + case.stepsize * \
1477 #                (a[2, 0] * Kgov1 + a[2, 1] * Kgov2)
1478 #
1479 #            # Generators.
1480 #            Kgen2 = case.generator(Xgen1, Xexc2, Xgov2, self.Pgen, Vgen1)
1481 #            Xgen2 = self.Xgen0 + case.stepsize * \
1482 #                (a[2, 0] * Kgen1 + a[2, 1] * Kgen2)
1483 #
1484 #            # Calculate system voltages.
1485 #            U2 = case.solveNetwork(Xgen2, self.Pgen, self.augYbus_solver,
1486 #                                   self.gbus)
1487 #
1488 #            # Calculate machine currents and power.
1489 #            Id2, Iq2, Pe2  = case.machineCurrents(Xgen2, self.Pgen,
1490 #                                                  U2[self.gbus])
1491 #
1492 #            # Update variables that have changed
1493 #            Vexc2 = abs(U2[self.gbus])
1494 #            Vgen2 = r_[Id2, Iq2, Pe2]
1495 #            Vgov2 = Xgen2[:, 1]
1496 #
1497 #            # K3 --------------------------------------------------------------
1498 #
1499 #            # Exciters.
1500 #            Kexc3 = case.exciter(Xexc2, self.Pexc, Vexc2)
1501 #            Xexc3 = self.Xexc0 + case.stepsize * \
1502 #                (a[3, 0] * Kexc1 + a[3, 1] * Kexc2 + a[3, 2] * Kexc3)
1503 #
1504 #            # Governors.
1505 #            Kgov3 = case.governor(Xgov2, self.Pgov, Vgov2)
1506 #            Xgov3 = self.Xgov0 + case.stepsize * \
1507 #                (a[3, 0] * Kgov1 + a[3, 1] * Kgov2 + a[3, 2] * Kgov3)
1508 #
1509 #            # Generators.
1510 #            Kgen3 = case.generator(Xgen2, Xexc3, Xgov3, self.Pgen, Vgen2)
1511 #            Xgen3 = self.Xgen0 + case.stepsize * \
1512 #                (a[3, 0] * Kgen1 + a[3, 1] * Kgen2 + a[3, 2] * Kgen3)
1513 #
1514 #            # Calculate system voltages
1515 #            U3 = case.solveNetwork(Xgen3, self.Pgen, self.augYbus_solver,
1516 #                                   self.gbus)
1517 #
1518 #            # Calculate machine currents and power.
1519 #            Id3, Iq3, Pe3 = case.machineCurrents(Xgen3, self.Pgen,
1520 #                                                 U3[self.gbus])
1521 #
1522 #            # Update variables that have changed
1523 #            Vexc3 = abs(U3[self.gbus])
1524 #            Vgen3 = r_[Id3, Iq3, Pe3]
1525 #            Vgov3 = Xgen3[:, 1]
1526 #
1527 #            # K4 --------------------------------------------------------------
1528 #
1529 #            # Exciters.
1530 #            Kexc4 = case.exciter(Xexc3, self.Pexc, Vexc3)
1531 #            Xexc4 = self.Xexc0 + case.stepsize * \
1532 #                (a[4, 0] * Kexc1 + a[4, 1] * Kexc2 + a[4, 2] * Kexc3 +
1533 #                 a[4, 3] * Kexc4)
1534 #
1535 #            # Governors.
1536 #            Kgov4 = case.governor(Xgov3, self.Pgov, Vgov3)
1537 #            Xgov4 = self.Xgov0 + case.stepsize * \
1538 #                (a[4, 0] * Kgov1 + a[4, 1] * Kgov2 + a[4, 2] * Kgov3 +
1539 #                 a[4, 3] * Kgov4)
1540 #
1541 #            # Generators.
1542 #            Kgen4 = case.generator(Xgen3, Xexc4, Xgov4, self.Pgen, Vgen3)
1543 #            Xgen4 = self.Xgen0 + case.stepsize * \
1544 #                (a[4, 0] * Kgen1 + a[4, 1] * Kgen2 + a[4, 2] * Kgen3 +
1545 #                 a[4, 3] * Kgen4)
1546 #
1547 #            # Calculate system voltages.
1548 #            U4 = case.solveNetwork(Xgen4, self.Pgen, self.augYbus_solver,
1549 #                                   self.gbus)
1550 #
1551 #            # Calculate machine currents and power.
1552 #            Id4, Iq4, Pe4 = case.machineCurrents(Xgen4, self.Pgen,
1553 #                                                 U4[self.gbus])
1554 #
1555 #            # Update variables that have changed.
1556 #            Vexc4 = abs(U4[self.gbus])
1557 #            Vgen4 = r_[Id4, Iq4, Pe4]
1558 #            Vgov4 = Xgen4[:, 1]
1559
1560 #------------------------------------------------------------------------------
1561 #  "RungeKuttaHighamHall" class:
1562 #------------------------------------------------------------------------------
1563
1564 class RungeKuttaHighamHall(RungeKuttaFehlberg):
1565     """ Runge-Kutta Higham and Hall ODE solver.
1566
1567     Based on RungeKuttaHighamHall.m from MatDyn by Stijn Cole, developed at
1568     Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/
1569     teaching/matdyn/} for more information.
1570     """
1571
1572     def __init__(self, t, Xgen0, Pgen, Vgen0, Xexc0, Pexc, Vexc0, Xgov0, Pgov,
1573             Vgov0, augYbus_solver, gbus, stepsize):
1574
1575         super(self, RungeKuttaHighamHall).__init__(t, Xgen0, Pgen, Vgen0,Xexc0,
1576             Pexc, Vexc0, Xgov0, Pgov, Vgov0, augYbus_solver, gbus, stepsize)
1577
1578         #: Runge-Kutta coefficients.
1579         self._a = array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0/9.0, 0.0, 0.0, 0.0,
1580                          0.0, 0.0, 1.0/12.0, 1.0/4.0, 0.0, 0.0, 0.0, 0.0,
1581                          1.0/8.0, 0.0, 3.0/8.0, 0.0, 0.0, 0.0, 91.0/500.0,
1582                          -27.0/100.0, 78.0/125.0, 8.0/125.0, 0.0, 0.0,
1583                          -11.0/20.0, 27.0/20.0, 12.0/5.0, -36.0/5.0, 5.0, 0.0,
1584                          1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0,
1585                          5.0/48.0])
1586         #: Runge-Kutta coefficients.
1587         self._b1 = array([1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0,
1588                           5.0/48.0, 0.0])
1589         #: Runge-Kutta coefficients.
1590         self._b2 = array([2.0/15.0, 0.0, 27.0/80.0, -2.0/15.0, 25.0/48.0,
1591                           1.0/24.0, 1.0/10.0,])
1592         # c = array([0.0, 2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0,])
1593
1594
1595     def solve(self):
1596         case = self.dyn_case
1597         b2 = self._b2
1598
1599         accept = False
1600         facmax = 4
1601         failed = False
1602
1603         i = 0
1604         while accept == False:
1605             i += 1
1606
1607             Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1, Vgen1, Kgen1, _ =\
1608                 self._k1()
1609
1610             Xexc2, Vexc2, Kexc2, Xgov2, Vgov2, Kgov2, Xgen2, Vgen2, Kgen2, _ =\
1611                 self._k2(Xexc1, Vexc1, Kexc1, Xgov1, Vgov1, Kgov1, Xgen1,
1612                     Vgen1, Kgen1)
1613
1614             Xexc3, Vexc3, Kexc3, Xgov3, Vgov3, Kgov3, Xgen3, Vgen3, Kgen3, _ =\
1615                 self._k3(Xexc2, Vexc2, Kexc2, Kexc1, Xgov2, Vgov2, Kgov2,
1616                     Kgov1, Xgen2, Vgen2, Kgen2, Kgen1)
1617
1618             Xexc4, Vexc4, Kexc4, Xgov4, Vgov4, Kgov4, Xgen4, Vgen4, Kgen4, _ =\
1619                 self._k4(Xexc3, Vexc3, Kexc3, Kexc2, Kexc1, Xgov3, Vgov3,
1620                     Kgov3, Kgov2, Kgov1, Xgen3, Vgen3, Kgen3, Kgen2, Kgen1)
1621
1622             Xexc5, Vexc5, Kexc5, Xgov5, Vgov5, Kgov5, Xgen5, Vgen5, Kgen5, _ =\
1623                 self._k5(Xexc4, Vexc4, Kexc4, Kexc3, Kexc2, Kexc1, Xgov4,
1624                     Vgov4, Kgov4, Kgov3, Kgov2, Kgov1, Xgen4, Vgen4, Kgen4,
1625                     Kgen3, Kgen2, Kgen1)
1626
1627             Xexc6, Vexc6, Kexc6, Xgov6, Vgov6, Kgov6, Xgen6, Vgen6, Kgen6, _ =\
1628                 self._k6(Xexc5, Vexc5, Kexc5, Kexc4, Kexc3, Kexc2, Kexc1,
1629                     Xgov5, Vgov5, Kgov5, Kgov4, Kgov3, Kgov2, Kgov1, Xgen5,
1630                     Vgen5, Kgen5, Kgen4, Kgen3, Kgen2, Kgen1)
1631
1632             Xexc7, _, Kexc7, Xgov7, _, Kgov7, Xgen7, _, Kgen7, U7 = \
1633                 self._k6(Xexc6, Vexc6, Kexc6, Kexc5, Kexc4, Kexc3, Kexc2,
1634                     Kexc1, Xgov6, Vgov6, Kgov6, Kgov5, Kgov4, Kgov3, Kgov2,
1635                     Kgov1, Xgen6, Vgen6, Kgen6, Kgen5, Kgen4, Kgen3, Kgen2,
1636                     Kgen1)
1637
1638             # Second, higher order solution.
1639             Xexc72 = self.Xexc0 + case.stepsize * \
1640                 (b2[0] * Kexc1 + b2[1] * Kexc2 + b2[2] * Kexc3 +
1641                  b2[3] * Kexc4 + b2[4] * Kexc5 + b2[5] * Kexc6 + b2[6] * Kexc7)
1642             Xgov72 = self.Xgov0 + case.stepsize * \
1643                 (b2[0] * Kgov1 + b2[1] * Kgov2 + b2[2] * Kgov3 +
1644                  b2[3] * Kgov4 + b2[4] * Kgov5 + b2[5] * Kgov6 + b2[6] * Kgov7)
1645             Xgen72 = self.Xgen0 + case.stepsize * \
1646                 (b2[0] * Kgen1 + b2[1] * Kgen2 + b2[2] * Kgen3 +
1647                  b2[3] * Kgen4 + b2[4] * Kgen5 + b2[5] * Kgen6 + b2[6] * Kgen7)
1648
1649             # Error estimate
1650             Xexc = abs((Xexc72 - Xexc7).T)
1651             Xgov = abs((Xgov72 - Xgov7).T)
1652             Xgen = abs((Xgen72 - Xgen7).T)
1653             errest = max( r_[max(max(Xexc)), max(max(Xgov)), max(max(Xgen))] )
1654
1655             if errest < EPS:
1656                 errest = EPS
1657
1658             q = 0.84 * (self.tol / errest)^(1.0 / 4.0)
1659
1660             if errest < self.tol:
1661                 accept = True
1662                 U0 = U7
1663                 Vgen0 = Vgen6;
1664                 Vgov0 = Vgov6
1665                 Vexc0 = Vexc6
1666
1667                 Xgen0 = Xgen6
1668                 Xexc0 = Xexc6
1669                 Xgov0 = Xgov6
1670
1671                 Pgen0 = self.Pgen
1672                 Pexc0 = self.Pexc
1673                 Pgov0 = self.Pgov
1674
1675 #                t = t0
1676             else:
1677                 failed += 1
1678                 facmax = 1
1679
1680 #                t = t0
1681                 Pgen0 = self.Pgen
1682                 Pexc0 = self.Pexc
1683                 Pgov0 = self.Pgov
1684
1685                 stepsize = min(max(q, 0.1), facmax) * case.stepsize
1686                 return
1687
1688             stepsize = min(max(q, 0.1), facmax) * case.stepsize
1689
1690             if stepsize > self.maxstepsize:
1691                 stepsize = self.maxstepsize
1692
1693         return Xgen0, Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0, Vgov0, \
1694             U0, errest, failed, self.t, stepsize
1695
1696
1697     def _k6(self, Xexc5, Vexc5, Kexc5, Kexc4, Kexc3, Kexc2, Kexc1, Xgov5,
1698             Vgov5, Kgov5, Kgov4, Kgov3, Kgov2, Kgov1, Xgen5, Vgen5, Kgen5,
1699             Kgen4, Kgen3, Kgen2, Kgen1):
1700         # Overrides the K6 method from the Runge-Kutta Fehlberg solver.
1701         case = self.dyn_case
1702         a = self._a
1703
1704         # Exciters.
1705         Kexc6 = case.exciter(Xexc5, self.Pexc, Vexc5)
1706         Xexc6 = self.Xexc0 + case.stepsize * \
1707             (a[6, 0] * Kexc1 + a[6, 1] * Kexc2 + a[6, 2] * Kexc3 + \
1708             a[6, 3] * Kexc4 + a[6, 4] * Kexc5 + a[6, 5] * Kexc6)
1709
1710         # Governors.
1711         Kgov6 = case.governor(Xgov5, self.Pgov, Vgov5)
1712         Xgov6 = self.Xgov0 + case.stepsize * \
1713             (a[6, 0] * Kgov1 + a[6, 1] * Kgov2 + a[6, 2] * Kgov3 + \
1714             a[6, 3] * Kgov4 + a[6, 4] * Kgov5 + a[6, 5] * Kgov6)
1715
1716         # Generators.
1717         Kgen6 = case.generator(Xgen5, Xexc6, Xgov6, self.Pgen, Vgen5)
1718         Xgen6 = self.Xgen0 + case.stepsize * \
1719             (a[6, 0] * Kgen1 + a[6, 1] * Kgen2 + a[6, 2] * Kgen3 + \
1720             a[6, 3] * Kgen4 + a[6, 4] * Kgen5 + a[6, 5] * Kgen6)
1721
1722         # Calculate system voltages.
1723         U6 = case.solveNetwork(Xgen6, self.Pgen, self.augYbus_solver,
1724                                self.gbus)
1725
1726         # Calculate machine currents and power.
1727         Id6, Iq6, Pe6 = case.machineCurrents(Xgen6, self.Pgen,
1728                                              U6[self.gbus])
1729
1730         # Update variables that have changed.
1731         Vexc6 = abs(U6[self.gbus])
1732         Vgen6 = r_[Id6, Iq6, Pe6]
1733         Vgov6 = Xgen6[:, 1]
1734
1735         return Xexc6, Vexc6, Kexc6, Xgov6, Vgov6, Kgov6, Xgen6, Vgen6, Kgen6,U6
1736
1737
1738     def _k7(self, Xexc6, Vexc6, Kexc6, Kexc5, Kexc4, Kexc3, Kexc2, Kexc1,
1739             Xgov6, Vgov6, Kgov6, Kgov5, Kgov4, Kgov3, Kgov2, Kgov1, Xgen6,
1740             Vgen6, Kgen6, Kgen5, Kgen4, Kgen3, Kgen2, Kgen1):
1741         case = self.dyn_case
1742         b1 = self.b1
1743
1744         # Exciters.
1745         Kexc7 = case.exciter(Xexc6, self.Pexc, Vexc6)
1746         Xexc7 = self.Xexc0 + case.stepsize * \
1747             (b1[0] * Kexc1 + b1[1] * Kexc2 + b1[2] * Kexc3 + b1[3] * Kexc4 +
1748              b1[4] * Kexc5 + b1[5] * Kexc6 + b1[6] * Kexc7)
1749
1750         # Governors.
1751         Kgov7 = case.governor(Xgov6, self.Pgov, Vgov6)
1752         Xgov7 = self.Xgov0 + case.stepsize * \
1753             (b1[0]* Kgov1 + b1[1] * Kgov2 + b1[2] * Kgov3 + b1[3] * Kgov4 +
1754              b1[4] * Kgov5 + b1[5] * Kgov6 + b1[6] * Kgov7)
1755
1756         # Generators.
1757         Kgen7 = case.generator(Xgen6, Xexc7, Xgov7, self.Pgen, Vgen6)
1758         Xgen7 = self.Xgen0 + case.stepsize * \
1759             (b1[0] * Kgen1 + b1[1] * Kgen2 + b1[2] * Kgen3 + b1[3] * Kgen4 +
1760              b1[4] * Kgen5 + b1[5] * Kgen6 + b1[6] * Kgen7)
1761
1762         # Calculate system voltages.
1763         U7 = case.solveNetwork(Xgen7, self.Pgen, self.augYbus_solver,
1764                                self.gbus)
1765
1766         # Calculate machine currents and power.
1767         Id7, Iq7, Pe7 = case.machineCurrents(Xgen7, self.Pgen, U7[self.gbus])
1768
1769         # Update variables that have changed
1770         Vexc7 = abs(U7[self.gbus])
1771         Vgen7 = r_[Id7, Iq7, Pe7]
1772         Vgov7 = Xgen7[:, 1]
1773
1774         return Xexc7, Vexc7, Kexc7, Xgov7, Vgov7, Kgov7, Xgen7, Vgen7, Kgen7,U7
1775
1776 #------------------------------------------------------------------------------
1777 #  "ModifiedEuler2" class:
1778 #------------------------------------------------------------------------------
1779
1780 class ModifiedEuler2(object):
1781     """ Modified Euler ODE solver with check on interface errors.
1782
1783     Based on ModifiedEuler2.m from MatDyn by Stijn Cole, developed at
1784     Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/electa/
1785     teaching/matdyn/} for more information.
1786     """
1787
1788     def __init__(self, t, Xgen0, Pgen, Vgen0, Xexc0, Pexc, Vexc0, Xgov0, Pgov,
1789             Vgov0, augYbus_solver, gbus, stepsize):
1790
1791         self.dyn_case = None
1792
1793         self.t = t
1794         self.Xgen0
1795         self.Pgen
1796         self.Vgen0
1797         self.Xexc0
1798         self.Pexc
1799         self.Vexc0
1800         self.Xgov0
1801         self.Pgov
1802         self.Vgov0
1803         self.augYbus_solver
1804         self.gbus
1805         self.stepsize
1806
1807
1808     def solve(self):
1809         case = self.dyn_case
1810
1811         eulerfailed = False
1812
1813         # First Prediction Step -----------------------------------------------
1814
1815         # Exciters.
1816         dFexc0 = case.exciter(self.Xexc0, self.Pexc, self.Vexc0)
1817         Xexc_new = self.Xexc0 + case.stepsize * dFexc0
1818
1819         # Governor.
1820         dFgov0 = case.governor(self.Xgov0, self.Pgov, self.Vgov0)
1821         Xgov_new = self.Xgov0 + case.stepsize * dFgov0
1822
1823         # Generators.
1824         dFgen0 = case.generator(self.Xgen0, Xexc_new, Xgov_new, self.Pgen,
1825                                 self.Vgen0)
1826         Xgen_new = self.Xgen0 + case.stepsize * dFgen0
1827
1828         Vexc_new = self.Vexc0
1829         Vgov_new = self.Vgov0
1830         Vgen_new = self.Vgen0
1831
1832         for i in range(self.maxit):
1833             Xexc_old = Xexc_new
1834             Xgov_old = Xgov_new
1835             Xgen_old = Xgen_new
1836
1837             Vexc_old = Vexc_new
1838             Vgov_old = Vgov_new
1839             Vgen_old = Vgen_new
1840
1841             # Calculate system voltages
1842             U_new = case.solveNetwork(Xgen_new, self.Pgen, self.augYbus_solver,
1843                                       self.gbus)
1844
1845             # Calculate machine currents and power.
1846             Id_new, Iq_new, Pe_new = case.machineCurrents(Xgen_new, self.Pgen,
1847                                                           U_new[self.gbus])
1848
1849             # Update variables that have changed.
1850             Vgen_new = r_[Id_new, Iq_new, Pe_new]
1851             Vexc_new = abs(U_new[self.gbus])
1852             Vgov_new = Xgen_new[:,1]
1853
1854             # Correct the prediction, and find new values of x ----------------
1855
1856             # Exciters.
1857             dFexc1 = case.exciter(Xexc_old, self.Pexc, Vexc_new)
1858             Xexc_new = self.Xexc0 + case.stepsize / 2.0 * (dFexc0 + dFexc1)
1859
1860             # Governors.
1861             dFgov1 = case.governor(Xgov_old, self.Pgov, Vgov_new)
1862             Xgov_new = self.Xgov0 + case.stepsize / 2.0 * (dFgov0 + dFgov1)
1863
1864             # Generators.
1865             dFgen1 = case.generator(Xgen_old, Xexc_new, Xgov_new, self.Pgen,
1866                                     Vgen_new)
1867             Xgen_new = self.Xgen0 + case.stepsize / 2.0 * (dFgen0 + dFgen1)
1868
1869             # Calculate error.
1870             Xexc_d = abs((Xexc_new - Xexc_old).T)
1871             Xgov_d = abs((Xgov_new - Xgov_old).T)
1872             Xgen_d = abs((Xgen_new - Xgen_old).T)
1873
1874             Vexc_d = abs((Vexc_new - Vexc_old).T)
1875             Vgov_d = abs((Vgov_new - Vgov_old).T)
1876             Vgen_d = abs((Vgen_new - Vgen_old).T)
1877
1878             errest = max( r_[max(max(Vexc_d)), max(max(Vgov_d)), max(max(Vgen_d)), max(max(Xexc_d)), max(max(Xgov_d)), max(max(Xgen_d))] )
1879
1880             if errest < self.tol:
1881                 break    # solution found
1882             else:
1883                 if i == self.maxit:
1884                     U0 = U_new
1885                     Vexc0 = Vexc_new
1886                     Vgov0 = Vgov_new
1887                     Vgen0 = Vgen_new
1888                     Xgen0 = Xgen_new
1889                     Xexc0 = Xexc_new
1890                     Xgov0 = Xgov_new
1891                     Pgen0 = self.Pgen
1892                     Pexc0 = self.Pexc
1893                     Pgov0 = self.Pgov
1894                     eulerfailed = True
1895                     return Xgen0, Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, \
1896                         Pgov0, Vgov0, U0, self.t, eulerfailed, case.stepsize
1897         # Update
1898         U0 = U_new
1899
1900         Vexc0 = Vexc_new
1901         Vgov0 = Vgov_new
1902         Vgen0 = Vgen_new
1903
1904         Xgen0 = Xgen_new
1905         Xexc0 = Xexc_new
1906         Xgov0 = Xgov_new
1907
1908         Pgen0 = self.Pgen
1909         Pexc0 = self.Pexc
1910         Pgov0 = self.Pgov
1911
1912         return Xgen0, Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0, Vgov0, \
1913             U0, self.t, eulerfailed, case.stepsize
1914
1915 #------------------------------------------------------------------------------
1916 #  "DynamicGenerator" class:
1917 #------------------------------------------------------------------------------
1918
1919 class DynamicGenerator(object):
1920     """ Defines a classical generator model and a fourth order generator model
1921     for dynamic simulation.
1922     """
1923
1924     def __init__(self, generator, exciter, governor, model=CLASSICAL, h=1.0,
1925                  d=0.01, x=1.0, x_tr=1.0, xd=1.0, xq=1.0, xd_tr=1.0, xq_tr=1.0,
1926                  td=1.0, tq=1.0):
1927         #: Power flow generator.
1928         self.generator = generator
1929
1930         #: Exciter model.
1931         self.exciter = exciter
1932
1933         #: Governor model.
1934         self.governor = governor
1935
1936         #: Classical or 4th order model.
1937         self.model = model
1938
1939         #: Inertia constant.
1940         self.h = h
1941
1942         #: Damping constant.
1943         self.d = d
1944
1945         # Classical model -----------------------------------------------------
1946
1947         #: Reactance (p.u.).
1948         self.x = x
1949
1950         #: Transient reactance (p.u.).
1951         self.x_tr = x_tr
1952
1953         # 4th order model -----------------------------------------------------
1954
1955         #: d-axis reactance (p.u.).
1956         self.xd = xd
1957
1958         #: q-axis reactance (p.u.).
1959         self.xq = xq
1960
1961         #: d-axis transient reactance (p.u.).
1962         self.xd_tr = xd_tr
1963
1964         #: q-axis transient reactance (p.u.).
1965         self.xq_tr = xq_tr
1966
1967         #: d-axis time constant (s).
1968         self.td = td
1969
1970         #: q-axis time constant (s).
1971         self.tq = tq
1972
1973 #------------------------------------------------------------------------------
1974 #  "Exciter" class:
1975 #------------------------------------------------------------------------------
1976
1977 class Exciter(object):
1978     """ Defines an IEEE DC1A excitation system.
1979     """
1980
1981     def __init__(self, generator, model=CONST_EXCITATION, ka=0.0, ta=0.0,
1982                  ke=0.0, te=0.0, kf=0.0, tf=0.0, aex=0.0, bex=0.0,
1983                  ur_min=-1.5, ur_max=1.5):
1984
1985         #: Power flow generator.
1986         self.generator = generator
1987
1988         #: Exciter model.
1989         self.model = model
1990
1991         #: Amplifier gain.
1992         self.ka = ka
1993
1994         #: Amplifier time constant.
1995         self.ta = ta
1996
1997         #: Exciter gain.
1998         self.ke = ke
1999
2000         #: Exciter time constant.
2001         self.te = te
2002
2003         #: Stabiliser gain.
2004         self.kf = kf
2005
2006         #: Stabiliser time constant.
2007         self.tf = tf
2008
2009         #: Parameter saturation function.
2010         self.aex = aex
2011
2012         #: Parameter saturation function.
2013         self.bex = bex
2014
2015         #: Lower voltage limit.
2016         self.ur_min = ur_min
2017
2018         #: Upper voltage limit.
2019         self.ur_max = ur_max
2020
2021 #------------------------------------------------------------------------------
2022 #  "Governor" class:
2023 #------------------------------------------------------------------------------
2024
2025 class Governor(object):
2026     """ Defines an IEEE model of a general speed-governing system for steam
2027     turbines. It can represent a mechanical-hydraulic or electro-hydraulic
2028     system by the appropriate selection of parameters.
2029     """
2030
2031     def __init__(self, generator, model=CONST_POWER, k=0.0, t1=0.0, t2=0.0,
2032                  t3=0.0, p_up=0.0, p_down=0.0, p_max=0.0, p_min=0.0):
2033
2034         #: Power flow generator.
2035         self.generator = generator
2036
2037         #: Governor model.
2038         self.model = model
2039
2040         #: Droop.
2041         self.k = k
2042
2043         #: Time constant.
2044         self.t1 = t1
2045
2046         #: Time constant.
2047         self.t2 = t2
2048
2049         #: Servo motor time constant.
2050         self.t3 = t3
2051
2052         #: Upper ramp limit.
2053         self.p_up = p_up
2054
2055         #: Lower ramp limit.
2056         self.p_down = p_down
2057
2058         #: Maxmimum turbine output.
2059         self.p_max = p_max
2060
2061         #: Minimum turbine output.
2062         self.p_min = p_min
2063
2064 #------------------------------------------------------------------------------
2065 #  "Event" class:
2066 #------------------------------------------------------------------------------
2067
2068 class Event(object):
2069     """ Defines an event.
2070     """
2071
2072     def __init__(self, time, type):
2073         """ Constructs a new Event instance.
2074         """
2075         #: Instant of change (s).
2076         self.time = time
2077
2078         #: Bus or branch event.
2079         self.type = type
2080
2081 #------------------------------------------------------------------------------
2082 #  "BusChange" class:
2083 #------------------------------------------------------------------------------
2084
2085 class BusChange(object):
2086     """ Defines a bus parameter change event.
2087
2088     Three-phase bus faults can be simulated by changing the shunt
2089     susceptance of the bus in a bus change event.
2090     """
2091
2092     def __init__(self, bus, time, param, newval):
2093
2094         #: Bus with new parameter value.
2095         self.bus = bus
2096
2097         #: Instant of change (s).
2098         self.time = time
2099
2100         #: Bus parameter to change.
2101         self.param = param
2102
2103         #: New parameter value.
2104         self.newval = newval
2105
2106 #------------------------------------------------------------------------------
2107 #  "BranchChange" class:
2108 #------------------------------------------------------------------------------
2109
2110 class BranchChange(object):
2111     """ Defines a branch parameter change event.
2112     """
2113
2114     def __init__(self, branch, time, param, newval):
2115
2116         #: Branch with new parameter value.
2117         self.branch = branch
2118
2119         #: Instant of change (s).
2120         self.time = time
2121
2122         #: Bus parameter to change.
2123         self.param = param
2124
2125         #: New parameter value.
2126         self.newval = newval
2127
2128 # EOF -------------------------------------------------------------------------