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